Cointegrating Polynomial Regressions with Power Law Trends: Environmental Kuznets Curve or Omitted Time Effects?

spacing{1.2} \begin{abstract} The environmental Kuznets curve predicts an inverted U-shaped relationship between environmental pollution and economic growth. Current analyses frequently employ models which restrict nonlinearities in the data to be explained by the economic growth variable only. We propose a Generalized Cointegrating Polynomial Regression (GCPR) to allow for an alternative source of nonlinearity. More specifically, the GCPR is a seemingly unrelated regression with (1) integer powers of deterministic and stochastic trends for the individual units, and (2) a common flexible global trend. We estimate this GCPR by nonlinear least squares and derive its asymptotic distribution. Endogeneity of the regressors will introduce nuisance parameters into the limiting distribution but a simulation-based approach nevertheless enables us to conduct valid inference. A multivariate subsampling KPSS test is proposed to verify the correct specification of the cointegrating relation. Our simulation study shows good performance of the simulated inference approach and subsampling KPSS test. We illustrate the GCPR approach using data for Austria, Belgium, Finland, the Netherlands, Switzerland, and the UK. A single global trend accurately captures all nonlinearities leading to a linear cointegrating relation between GDP and CO\textsubscript{2} for all countries. This suggests that the environmental improvement of the last years is due to economic factors different from GDP. JEL Classification: C12, C13, C32, O44, Q20 Keywords: Cointegration Testing, Environmental Kuznets Curve, Generalized Cointegrating Polynomial Regression, Power Law Trends \end{abstract}

spacing{2} \section{Introduction} On page 370 of their seminal paper, grossmankrueger1995 conclude: \begin{displayquote} “Contrary to the alarmist cries of some environmental groups, we find no evidence that economic growth does unavoidable harm to the natural habitat. Instead we find that while increases in GDP may be associated with worsening environmental conditions in very poor countries, air and water quality appear to benefit from economic growth once some critical level of income has been reached.” \end{displayquote} The quote above suggests an inverted U-shaped relationship between environmental degradation and economic growth. This relationship is currently known as the Environmental Kuznets Curve (EKC) and it forms an active research area. Its relevance becomes clear if we look at some forecasts of long-run economic growth. The projected GDP per capita growth of the world is about 2.1% per year for the next decades (chapter 3 in nordhaus2013; christensengillinghamnordhaus2018) and this growth is partially powered by carbon-based energy resources, water usage, and material consumption. In absence of an EKC, economic growth will place more and more stress on the environment. Alternatively, if the EKC exists, then the inverted U-shape eventually implies a turning point after which economic growth and environmental improvement go hand in hand. Due to such considerations, there is now, some 25 years after its first conception, a rich literature that (1) reports on the experimental evidence on the existence/nonexistence of the EKC, (2) provides economic theory to explain the EKC, and/or (3) refines the econometric tools that are used to analyse the EKC.\footnote{Further references to these specific areas of research can be found in the review articles by dasguptaetal2002, stern2004, and carson2009 among others.} To quantify the volume of the literature we have entered the search query “Environmental Kuznets Curve” into the Web of Science: more than 4,200 references are found.\footnote{Web of Science, accessed on December 6, 2021, http://www.webofknowledge.com.} Driven by contradictory empirical results as well as the variability in estimated turning points, the EKC has been criticised on two main points. First, the income variable was initially treated as a stationary variable whereas later research shows that the unit root hypothesis often cannot be rejected (see galeottimaneralanza2009, p. 553; stern2017, p. 14--15). Nonstationarity has further implications because EKC regressions include higher integer powers of GDP as well. This combination of nonstationarity and nonlinearity places the EKC in the nonlinear cointegration literature and appropriate econometric techniques should be employed. Such techniques have been developed in wagner2015 and wagnerhong2016 under the name of Cointegrating Polynomial Regressions (CPRs). CPRs contain deterministic variables, integrated processes, and their integer powers.\footnote{stypkawagnergrabarczykkawka2017 reiterate the need to model the income variable as nonstationary. It is less important to use an estimation procedure acknowledging the fact that several integer powers of the same integrated process appear as regressors. That is, stypkawagnergrabarczykkawka2017 find that the “standard estimator” which treats higher order powers of the integrated regressor as additional I(1) variables has the same limiting distribution as the CPR estimator.} Multivariate generalizations of CPRs, Seemingly Unrelated Cointegrating Polynomial Regressions (SUCPRs), are discussed in wagnergrabarczykhong2019 and linreuvers2019. As a second point of critique, there is an ongoing debate on the model specification. Various functional forms can describe the relationship between national income and the pollution variable. The quadratic specification is widespread but cubic relationships (harbaughlevinsonwilson2002; wagner2015) and double-nonlinear transformation (lintuyao2020) are also in use. Various specification tests are helpful while deciding on the right parametric specification (hongphillips2010; wangphillips2012; wangwuzhu2018). Alternatively, one could resort to nonparametric estimation procedures altogether (wangphillips2009; lintonwang2016). Whereas such modelling approaches do allow for a more flexible relationship, they also implicitly assume that nonlinear environmental effects are solely attributable to economic growth. Relevant variables are thus potentially missing from the model specification. Such omitted variables are a valid concern because advances in green technology, pollution policy, and environmental awareness, may all influence pollution levels. However, such data is typically available for short time spans only (and for that reason often excluded from the model). Time effects can control for time-variation in unobserved effects (volleberghmelenbergdijkgraaf2009). Current developments on nonlinear cointegration emphasize the role of the nonstationary regressor yet pay less attention to time effects. Time effects are important. The small simulation exercise in Table (ref) illustrates the point. Foreshadowing our proposed model, we consider a multivariate setting with a global nonlinear, smooth time trend. The global trend is omitted by the researcher and a quadratic EKC-specification is estimated: $y_{i,t} = \tau_{1,i} + \tau_{2,i} t + \phi_{1,i} x_{i,t} + \phi_{2,i} x_{i,t}^2 + u_{i,t}$ ($i=1,\ldots,3$), where $x_{i,t}$ and $y_{i,t}$ are unit-specific variables measuring income and environmental pollution, respectively. We test $H_0: \phi_{2,1}=\phi_{2,2}=\phi_{2,3} =0$ because a significantly negative coefficient in front of $x_{i,t}^2$ is typically interpreted as evidence of an EKC.\footnote{For the moment, we focus on the curvature parameter. Clearly, for an inverted U-shaped relationship the coefficient in front of the linear term should be positive.} Panel (A) reveals exacerbated rejection frequency as curvature caused by the global deterministic trend is mistakenly interpreted as curvature caused by the income variable. In other words, negative and significant coefficients in front of squared GDP are possibly caused by omitted nonlinear deterministic trends rather than being indicative of an EKC. To be on the safe side, we recommend researchers to include a nonlinear trend component in their model specification. If unnecessary, then this is rather innocuous. Indeed, Panel (B) of Table (ref) shows that significant results for nonlinear economic growth effects continue to be found with modest losses in statistical power. \begin{table}[t] \caption{The rejection rate (in $\%$) when testing $H_0: \phi_{2,1}=\phi_{2,2}=\phi_{2,3}=0$. (A) Falsely inflated rejections of $H_0: \phi_{2,1}=\phi_{2,2}=\phi_{2,3}=0$ when time effects are omitted. (B) Adding an additional global deterministic trend to the model specification hardly influences the power of the test $H_0: \phi_{2,1}=\phi_{2,2}=\phi_{2,3}=0$. That is, significant coefficients in front of $x_{i,t}^2$ remain significant after adding a redundant flexible global trend.} \resizebox{\textwidth}{!}{ \begin{threeparttable} \begin{tabular}{ccccccccc} \toprule & \multicolumn{3}{c}{Panel (A): Omitted Global Trend} & & & \multicolumn{3}{c}{Panel (B): Redundant Global Trend} \\ \midrule \addlinespace[0.1cm] DGP & \multicolumn{3}{c}{$y_{i,t}= \tau_g t^{\theta} + \tau_{1,i} + \tau_{2,i} t + \phi_{1,i} x_{i,t} + u_{i,t}$} & & \multicolumn{1}{l} & \multicolumn{3}{c}{$\text{\small{$y_{i,t}$}}=\text{\small{$\tau_{1,i} + \tau_{2,i} t + \phi_{1,i} x_{i,t}$}} + \phi_2 x_{i,t}^2 + \text{\small{$u_{i,t}$}}$} \\ \addlinespace[0.1cm] \cmidrule{2-4}\cmidrule{7-9} \addlinespace[0.1cm] Model & \multicolumn{3}{c}{$y_{i,t} = \tau_{1,i} + \tau_{2,i} t + \phi_{1,i} x_{i,t} + \phi_{2,i} x_{i,t}^2 + u_{i,t}$} & & \multicolumn{1}{l} & \multicolumn{1}{l}{ Correct Specification} & \multicolumn{1}{l} & \multicolumn{1}{l}{$\text{ \small{$y_{i,t}$}} = \tau_g t^{\theta} + \text{ \small{$\tau_{1,i} + \tau_{2,i} t + \phi_{1,i} x_{i,t} + \phi_{2,i} x_{i,t}^2 + u_{i,t}$}}$} \\ \addlinespace[0.1cm] \midrule $\tau_g\,(\times 10^{-5})$ & FM-SOLS & & FM-SUR & & $\phi_2$ & SimNLS & & SimNLS \\ \midrule 0 & 6.30 & & 6.63 & & 0 & 6.93 & & 5.70 \\ \addlinespace[0.1cm] -0.5 & 13.27 & & 12.77 & & -0.5 & 9.20 & & 6.97 \\ -1 & 30.07 & & 27.50 & & -1 & 14.97 & & 9.97 \\ -1.5 & 46.23 & & 41.57 & & -1.5 & 30.47 & & 20.70 \\ \addlinespace[0.1cm] -2 & 56.60 & & 50.30 & & -2 & 55.33 & & 40.53 \\ -2.5 & 64.50 & & 56.00 & & -2.5 & 81.73 & & 69.20 \\ -3 & 68.60 & & 57.57 & & -3 & 93.97 & & 89.50 \\ \bottomrule \end{tabular} \begin{tablenotes} • Note 1: For illustrative purpose, we consider a stylised example in this introduction. The exact parametrisation is available in Section of the Supplementary Material. More elaborate simulation results based on the empirical application are reported as simulation DGP2 in Section (ref). • Note 2: FM-SOLS and FM-SUR are documented in wagnergrabarczykhong2019. The results in Panel (B) are based on simulation-based inference, see Section (ref). \end{tablenotes} \end{threeparttable} } \end{table} The contributions of this paper are fourfold. First, we propose the Generalized Cointegrating Polynomial Regression (GCPR). This multivariate model features a global power law trend to capture time effects. Power law trends have been employed to model non-constant growth rates in technology indices (duggalsaltzmanklein1999; duggal2007infrastructure) and production functions (klein2004). Within the EKC context, this flexible trend can capture common time effects that are implicit in omitted variables. Alternatively, as in lilinton2020, the reader can view the global flexible trend as an outside option (next to the income variable) to describe nonlinearities in the data. Limiting distributions for estimators in models with purely deterministic power law trends have been reported in phillips2007, robinson2012 and gaolintonpeng2020. The presence of integrated variables requires an alternative asymptotic framework. Moreover, due to endogeneity, approaches assuming pre-determined integrated regressors (parkphillips1999; parkphillips2001; changparkphillips2001) are inappropriate and we instead opt for a proof along the lines of chanwang2015. The resulting limiting distribution is non-standard because (1) the scaling matrix with convergence rates is non-diagonal and parameter dependent, and (2) second-order bias terms are present. Second, we propose a simulation-based approach to conduct inference. Monte Carlo simulations show clear benefits of this simulation-based approach in terms of size control compared to existing methods. Third, in the spirit of choisaikkonen2010, we report a multivariate KPSS-type test to verify the stationarity of the error process thus enabling researchers to avoid spurious results or misspecified cointegrating relations. Fourth and finally, in the empirical application, we investigate the EKC for Austria, Belgium, Finland, the Netherlands, Switzerland and the UK over the period 1870--2014. Nonparametric estimates and tests confirm that the global trend captures all nonlinearities in the data. Nonlinear effects in log per capita GDP (and thus also evidence for an EKC) are absent. Given such findings, we offer a clear recommendation to researchers to check whether their EKC conclusions are robust to the inclusion of power law trends. This paper is organized as follows. Section (ref) introduces the model and the estimation framework. Asymptotic properties of the estimators and parameter inference are discussed in Section (ref). The Monte Carlo simulations in Section (ref) compare asymptotic results to finite sample performance. An in depth discussion of the Environment Kuznets Curve can be found in Section (ref). Section (ref) concludes. The proofs of the main theorems are collected in the Appendix and further information is available in the Supplementary Material. Finally, some words on notation. The integer part of the number $a\in \SR^{+}$ is denoted by $[a]$. For a vector $\vx\in\SR^n$, its $p$-norm is denoted by $\|\vx\|_p=(\sum_{i=1}^{n}|x_i|^p)^{1/p}$. For a matrix $\mA$, say of dimension ($n\times m$), the induced $p$-norm is defined as $\|\mA\|_p=\sup_{\vx\neq \vzeros} \|\mA\vx\|_p/\|\vx\|_p$. We will omit the subscripts whenever $p=2$. The $(n\times n)$ identity matrix is denoted $\mI_n$ and $\vones_n$ signifies an $n$-dimensional column vector with all entries equal to 1. The block-diagonal matrix $\diag[\mA_1,\ldots,\mA_n]$ stacks the matrices $\mA_1,\ldots,\mA_n$ along its diagonal. We omit the integration bounds whenever the integration interval is $[0,1]$. The symbol “$\stackrel{d}{=}$” stands for equality in distribution, and “$\pto$” and “$\dto$” denote convergence in probability and in distribution. If convergence occurs conditionally on the sample, then we add a superscript “*” to the standard notation. The probabilistic Landau symbols are $O_p(\cdot)$ and $o_p(\cdot)$. Finally, the generic constant $C$ can change from line to line. \section{The Model and NLS Estimation} Our model specification enriches the Seemingly Unrelated Cointegrating Polynomial Regressions (SUCPRs) from wagnergrabarczykhong2019 with a flexible deterministic trend. That is, each individual series in the system is affected by specific deterministic variables and integrated regressors (and their integer powers) while a global flexible trend describes nonlinear behaviour that is prevalent across all series. The resulting \emph{Generalized Cointegrating Polynomial Regression (GCPR)} is given by: \begin{equation} y_{i,t} = \tau_g t^{\theta} + \tau_{1,i} + \tau_{2,i} t + \sum_{j=1}^{p_i} \phi_{j,i} x_{i,t}^j + u_{i,t}, \qquad i=1,\ldots,N,\qquad t=1,\ldots,T, \end{equation} where $\theta\in \Theta(\varepsilon)$ with $\Theta(\varepsilon)=\left\{\theta\in [\theta_L,\theta_U]:~ |\theta|>\epsilon,\, |\theta-1|>\epsilon \right\}$ and $-1<\theta_L\leq \theta_U<\infty$. Alternatively, we write $y_{i,t} = \tau_g t^\theta + \vz_{i,t}' \vbeta_i + u_{i,t}$, where $\vz_{i,t}=\big[1,t,x_{i,t},\ldots,x_{i,t}^{p_i}\big]'$ and $\vbeta_i=\big[\tau_{1,i},\tau_{2,i},\phi_{1,i},\ldots,\phi_{p_i,i}\big]'$. Finally, we stack all $N$ equations in (ref) in matrix form to retrieve \begin{equation} \vy_t= \tau_{g} t^{\theta} \vones_N+ \mZ_t'\vbeta + \vu_t, \qquad t=1,\ldots,T, \end{equation} with $\vy_t=\big[y_{1,t},\ldots,y_{N,t}\big]'$, $\mZ_t=\diag\big[\vz_{1,t},\ldots,\vz_{N,t}\big]$, and the vector $\vbeta=\big[\vbeta_1',\ldots,\vbeta_N'\big]'$ of length $p=2N+\sum_{i=1}^{N}p_i$ containing all local parameters. We consider nonlinear least squares (NLS) estimators of the unknown parameters in (ref). As such, we define the objective function $Q_T(\theta,\tau_g,\vbeta)= \frac{1}{2} \sum_{t=1}^{T} \big\| \vy_t - \tau_g t^\theta \vones_N - \mZ_t'\vbeta \big\|^2$ and compute \begin{equation} \left(\,\widehat{\theta}_T, \widehat\tau_{g,T},\widehat\vbeta_T \right) = \operatornamewithlimits{arg\;min}_{(\theta,\,\tau_g,\,\vbeta)\in \Theta(\varepsilon) \times\SR \times\SR^p} Q_T(\theta,\tau,\vbeta). \end{equation} The optimization problem in (ref) is easy to solve. For any given $\theta$, model (ref) is linear-in-parameters and the minimizers for $\tau_g$ and $\vbeta$ can be found from an OLS regression by constructing $\mZ_t'(\theta)= \big[ t^\theta\vones_N \;\mZ_t'\big]$ and computing $ \left[ \begin{smallmatrix} \tau_g(\theta) \\ \vbeta(\theta) \end{smallmatrix} \right] = \left( \sum_{t=1}^T \mZ_t^{}(\theta)\mZ_t'(\theta) \right)^{-1} \left( \sum_{t=1}^T \mZ_t(\theta) \vy_t \right) $. We subsequently minimize the concentrated criterion function $\widetilde{Q}_T(\theta)=Q_T\big(\theta,\tau_g(\theta),\vbeta(\theta)\big)$ to obtain $\widehat{\theta}_T$. At last, we plug in $\widehat{\theta}_T$ and recover $\widehat\tau_{g,T}$ and $\widehat\vbeta_T$ through a final OLS estimation. \begin{remark} Keeping the powers of $x_{i,t}$ fixed allows us to test for their significance and thereby distinguish between nonlinearities caused by deterministic and stochastic trends. This is important for our empirical application on the Environmental Kuznets Curve, see Section (ref). huphillipswang2019 study a model with a flexible power of the integrated regressor. That is, these authors derive the limiting distribution of the NLS estimators for $\beta$ and $\gamma$ when $y_t=\beta |x_t|^\gamma+u_t$ with $\beta\neq 0$. \end{remark} \begin{remark} The GCPR of (ref) can be extended in several directions. First, integer powers of deterministic trends can be added as long as $\Theta(\varepsilon)$ is adjusted accordingly (to avoid collinearity). Second, multiple explanatory variables can be included. Related to the EKC, the literature suggests examples such as: population density (seldensong1994), trade openness (jalilferdiun2011), energy prices (almulaliozturk2016), and educational level (maranzanobentomanera2021). For nonstationary variables, conditions similar to those on $\{x_{i,t}\}$ should be fulfilled (Assumption (ref) below). Stationary variables should be strictly exogenous. Avoiding the elaborate notation which would otherwise arise, we focus on the baseline specification in (ref). \end{remark} \section{Asymptotic Theory} We subsequently study the asymptotic properties of the NLS estimators. To this end we first collect all the unknown parameters in the vector $\vgamma= \big[\theta,\tau_g,\vbeta' \big]'$. This vector is assumed to be an element of the parameter space $\bm \Gamma = \Theta(\varepsilon)\times\SR^{1+p}$. The true parameter vector is $\vgamma_0= \big[\theta_0,\tau_{g,0},\vbeta_0' \big]'$. \begin{assumption} The global trend is relevant, i.e. $\tau_{g,0}\neq 0$. \end{assumption} \begin{assumption} Let $\vzeta_t=[\eta_t',\vepsi_t']'$ be a sequence of i.i.d. random vectors with $\E(\vzeta_t)=\vzeros$, $\mSigma=\E\big(\vzeta_t^{}\vzeta_t')\succ 0$, and $\E\left\|\vzeta_t\right\|^{2q}<\infty$ for some $q>2$. \begin{enumerate}[(a)] • $u_t = \sum_{k=0}^\infty \psi_k \eta_{t-k}$ with $\sum_{k=1}^\infty k | \psi_k | < \infty$. • $\vx_t = \sum_{s=1}^t \vv_s$, where $\vv_t = \sum_{k=0}^\infty \mPsi_k \vepsi_{t-k}$ with $\sum_{k=0}^\infty \|\mPsi_k\| < \infty$ and $\det\left( \sum_{k=0}^\infty \mPsi_k \right) \neq 0$. \end{enumerate} \end{assumption} The first assumption is needed to avoid identification issues. That is, if $\tau_{g,0}=0$, then $\theta$ is not identified and the Davies problem arrises when testing $H_0: \tau_i = 0$ (see davies1977,davies1987). Such complications are not investigated here and this is further reflected in our model specification (ref). That is, we consider \emph{flexible} powers of the deterministic trends but \emph{fixed} powers of the stochastic trends, hence allowing us to test zero restrictions on (elements of) $\vbeta$. This is of crucial importance in the EKC application while determining whether nonlinear effects in the economic growth variables ($x_{i,t}$) remain significant after nonlinear time trends have been added to the model. Assumption (ref) has been relaxed in the literature albeit for different models. baekchophillips2015 and chophillips2018 study the asymptotic behaviour of a quasi-likelihood ratio test when Assumption (ref) is violated and the conditional mean of the data contains strictly stationary regressors and a flexible time trend. Alternatively, one can use drifting parameter sequences with different identification strengths as in andrewscheng2012. Assumption (ref) excludes cointegration among elements of $\vx_t$ and defines this vector as the partial sum of a short memory process. The latter implies that $ \frac{1}{\sqrt{T}}\sum_{t=1}^{[rT]} \left[ \begin{smallmatrix} \vu_t\\ \vv_t \end{smallmatrix} \right] \longrightarrow_{d} \bm{B}(r)= \left[ \begin{smallmatrix} \bm{B}_u(r)\\ \bm{B}_v(r) \end{smallmatrix} \right] $ where $\bm{B}(r)$ denotes an $2N$-dimensional vector Brownian motion with covariance matrix $\mOmega = \left[\begin{smallmatrix} \mOmega_{uu} & \mOmega_{uv}\\ \mOmega_{vu} & \mOmega_{vv} \end{smallmatrix}\right]$. The one-sided long-run covariance matrix $ \mDelta= \sum_{h=0}^\infty \E\left( \left[\begin{smallmatrix} \vu_t \vu_{t+h} & \vu_t \vv_{t+h}' \\ \vv_t \vu_{t+h} & \vv_t^{}\vv_{t+h}' \end{smallmatrix} \right] \right) = \left[\begin{smallmatrix} \mDelta_{uu} & \mDelta_{uv}\\ \mDelta_{vu} & \mDelta_{vv} \end{smallmatrix}\right]$ is partitioned similarly. Subscripts refer to specific elements. For example, $\bm{B}_{v_i}$ and $\mDelta_{v_i u_j}$ denote the $i$\textsuperscript{th} and $(i,j)$\textsuperscript{th} elements of $\bm{B}_v$ and $\mDelta_{vu}$, respectively. A concise exposition of our results asks for additional notation. An enumeration of various definitions is presented below. \begin{enumerate}[(1)] • Introduce $\mD_{(i),T}=\diag\big[1,T,T^{1/2},T,\ldots,T^{p_i/2}\big]$ to scale the deterministic and stochastic trends within each equation. For the full system of equation, define $\mD_{Z,T}=\diag\left[\mD_{(1),T},\ldots,\mD_{(N),T}\right]$, $\mD_{\theta_{0},T}=\sqrt{T}\left[ \begin{smallmatrix} T^{\theta_{0}} & & \\ & T^{\theta_{0}} & \\ & & \mD_{Z,T} \end{smallmatrix}\right]$ and $\mL_{\tau_{g,0},T}=\left[ \begin{smallmatrix} 1 & -\tau_{g,0}\ln{T} & \\ 0 & 1 & \\ & & \mI_{p} \end{smallmatrix}\right]$. Finally, set $\mG_{\vgamma_0,T}=\mD_{\theta_0,T}^{}\mL_{\tau_{g,0},T}^{\prime-1}$. • Define $\vj_i(r)=\left[1,r,B_{v_i}(r),B_{v_i}^2(r),\ldots,B_{v_i}^{p_i}(r)\right]'$, $\mJ_Z(r)=\diag\left[\vj_1(r),\ldots,\vj_N(r)\right]$, and $\mJ(r;\vgamma_0)=\Big[\tau_{g,0} r^{\theta_0}\ln{r}\,\vones_N, r^{\theta_0}\vones_N, \mJ_Z'(r)\Big]'$. • For the second-order bias terms, we define $\vb_i= \Big[\vzeros_{1\times 2},1,2\int \bm{B}_{v_i}(r)dr,\ldots,p_i\int \bm{B}_{v_i}^{p_i-1}(r)dr \Big]'$ and $\bm{\calB}_{vu}=\big[\vzeros_{1\times 2},\vb_1' \mDelta_{v_1 u_1} ,\dots,\vb_N ' \mDelta_{v_N u_N}\big]'$. \end{enumerate} \begin{theorem} Under Assumptions (ref)-(ref), we have \begin{equation*} \mG_{\vgamma_0,T}\big(\widehat{\vgamma}_T-\vgamma_0\big) \longrightarrow_{d} \left(\int \mJ(r;\vgamma_0)\mJ'(r;\vgamma_0) dr\right)^{-1}\left(\int \mJ(r;\vgamma_0) d\bm{B}_{u}(r)+{\bm{\mathcal{B}}}_{vu}\right) =: {\bm \calJ}(\vgamma_0), \end{equation*} as $T\rightarrow\infty$ and $N$ fixed. \end{theorem} The proof of Theorem (ref) is closely related to the work by chanwang2015. These authors provide the asymptotic distribution of NLS estimators under a set of general conditions in univariate, nonstationary time series models (see their theorem 3.1). The results in chanwang2015 and wangwuzhu2018 suggest that Assumption (ref) can be replaced by a long memory specification for $\diff \vx_t$. However, long memory parameters will enter the limiting distribution and inference will be complicated further. We illustrate Theorem (ref) with two examples. These examples highlight the two mathematical features that complicate parameter inference. \begin{example} We consider $y_t= \tau t^\theta+ u_t$ with innovations satisfying Assumption (ref). The limiting distribution of the parameter estimators depends solely on the mean square Riemann-Stieltjes integrals $\int \tau_0 r^{\theta_0} \ln(r) dB_u$ and $\int r^{\theta_0} dB_u$, and is therefore normally distributed (e.g., section 2.3 in tanaka2017). We have \begin{equation} \begin{bmatrix} T^{\theta_0+\frac{1}{2}} & 0\\ T^{\theta_0+\frac{1}{2}} \tau_0 \ln(T) & T^{\theta_0+\frac{1}{2}} \end{bmatrix} \begin{bmatrix} \,\widehat{\theta}_T -\theta_0 \\ \,\widehat{\tau}_T - \tau_0 \end{bmatrix} \longrightarrow_{d} \rN \left( \vzeros, \Omega_{uu}(2\theta_0+1)^3 \begin{bmatrix} 2 \tau_0^2 & - \tau_0 (2\theta_0+1) \\ - \tau_0 (2\theta_0+1) & (2\theta_0+1)^2 \end{bmatrix}^{-1} \right). \end{equation} The scaling matrix in the LHS of (ref), $\left[ \begin{smallmatrix} T^{\theta_0+\frac{1}{2}} & 0\\ T^{\theta_0+\frac{1}{2}} \tau_0 \ln(T) & T^{\theta_0+\frac{1}{2}} \end{smallmatrix}\right]$, depends on $\theta_0$ and is non-diagonal. The dependence on $\theta_0$ is unavoidable but asymptotic results for the case of a diagonal scaling matrix are obtainable at the expense of a singular joint distribution. That is, noting that $ \left[ \begin{smallmatrix} T^{\theta_0+\frac{1}{2}} & 0\\ 0 & T^{\theta_0+\frac{1}{2}}/\ln(T) \end{smallmatrix} \right] = \left[ \begin{smallmatrix} 1 & 0 \\ -\tau_0 & 1/\ln(T) \end{smallmatrix} \right] \left[\begin{smallmatrix} T^{\theta_0+\frac{1}{2}} & 0\\ T^{\theta_0+\frac{1}{2}} \tau_0 \ln(T) & T^{\theta_0+\frac{1}{2}} \end{smallmatrix} \right] $ and since $\lim_{T\to \infty} \left[ \begin{smallmatrix} 1 & 0 \\ -\tau_0 & 1/\ln(T) \end{smallmatrix} \right] = \left[ \begin{smallmatrix} 1 & 0 \\ -\tau_0 & 0 \end{smallmatrix} \right] $, the continuous mapping theorem implies $ \left[\begin{smallmatrix} T^{\theta_0+\frac{1}{2}} & 0\\ 0 & T^{\theta_0+\frac{1}{2}}/\ln(T) \end{smallmatrix} \right] \left[ \begin{smallmatrix} \widehat{\theta}_T - \theta_0 \\ \widehat{\tau}_T - \tau_0 \end{smallmatrix} \right] \dto \left[\begin{smallmatrix} 1/\tau_0 \\ -1 \end{smallmatrix} \right] \times \rN \left( \vzeros, \Omega_{uu}(2\theta_0+1)^3 \right). $ This limiting distribution coincides with the result in theorem 6.3 of phillips2007. \end{example} \begin{example} If $y_t = \tau t^\theta+ \phi x_t + u_t$, then the limiting distribution of the NLS estimator is: \begin{equation*} \resizebox{1\hsize}{!}{$ \begin{bmatrix} T^{\theta_0+\frac{1}{2}} \\ T^{\theta_0+\frac{1}{2}} \tau_0 \ln(T) & T^{\theta_0+\frac{1}{2}} \\ & & T \end{bmatrix} \begin{bmatrix} \widehat{\theta}_T -\theta_0 \\ \widehat{\tau}_T - \tau_0 \\ \widehat{\phi}_T - \phi_0 \end{bmatrix} \longrightarrow_{d} \begin{bmatrix} \int \big(\tau_0 r^{\theta_0}\ln (r)\big)^2 dr & \int \tau_0 r^{2\theta_0}{\ln (r)} dr & \int \tau_0 r^{\theta_0} \ln (r) B_v dr \\ \int \tau_0 r^{2\theta_0}{\ln (r)} dr & \int r^{2\theta_0} dr & \int r^{\theta_0} B_v dr \\ \int \tau_0 r^{\theta_0} \ln (r) B_v dr & \int r^{\theta_0} B_v dr & \int B_v^2 dr \end{bmatrix}^{-1} \times \\ \left( \begin{bmatrix} \int \tau_0 r^{\theta_0} \ln(r) dB_u\\ \int r^{\theta_0} dB_u\\ \int B_v dB_u \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ \Delta_{vu} \end{bmatrix} \right). $} \end{equation*} This limiting distribution exhibits second order bias when $ \Delta_{vu}\neq 0$, or when $B_u$ and $B_v$ are correlated. \end{example} Two features of the limiting distribution of $\mG_{\vgamma_0,T}\big(\widehat{\vgamma}_T-\vgamma_0\big)$ deserve further comments. First, as emphasised in Examples (ref)--(ref), the scaling matrix $\mG_{\vgamma_0,T}$ features two less common properties: (1) this matrix depends on the true parameters $\vtau_{g,0}$ and $\theta_0$, and (2) $\mG_{\vgamma_0,T}$ is not diagonal. These peculiarities are caused by the nonlinearity and nonstationarity of the model. More specifically, these features can be traced back to the presence of the global trend. Limiting distributions with a similar mathematical structure can be found in the structural breaks literature, cf. model setting II.b of perronzhu2005 and its detailed analysis in beutnerlinsmeekes2019. Second, the nonstationary regressor $x_{i,t}$ enters the model (ref) through a polynomial transformation of the form $g(x_{i,t},\vphi_i)= \phi_{i,1} x_{i,t} + \ldots + \phi_{i ,p_i} x_{i,t}^{p_i}$ ($i=1,2,\ldots,N$). In the terminology of parkphillips2001 this part of the regression function is a linear combination of $H_0$-regular functions. It is well-documented in the literature, e.g. changparkphillips2001 and chanwang2015, that this leads to second-order bias terms and hence nonstandard inference (except for the special case of strictly exogenous nonstationary regressors). \subsection{Consistent Long-Run Covariance Matrix Estimation} Correcting for second-order bias terms typically involves estimating long-run variance (LRV) matrices. This subsection establishes that the NLS residuals can be used to construct consistent kernel estimators for the LRV matrices $\mDelta$ and $\mOmega$. Defining $\bm V_t(\vgamma)=\big[\vu_t(\vgamma)', \diff \vx_t']'$ with $\vu_t(\vgamma)=\vy_t - \tau_{g} t^{\theta} \viota_N - \mZ_t'\vbeta $, these LRV estimators are defined as \begin{equation} \widehat{\mDelta}_T = \frac{1}{T} \sum_{t=1}^T \sum_{s=1}^t k\left(\frac{|t-s|}{b_T} \right) \bm V_t(\,\widehat{\vgamma}_T) \bm V_t(\,\widehat{\vgamma}_T)', \qquad \widehat{\mOmega}_T = \frac{1}{T} \sum_{t=1}^T \sum_{s=1}^T k\left(\frac{|t-s|}{b_T} \right) \bm V_t(\,\widehat{\vgamma}_T) \bm V_t(\,\widehat{\vgamma}_T)', \end{equation} for some kernel function $k(\cdot)$ and bandwidth parameter $b_T$. The first $N$ elements of $\bm V_t(\,\widehat{\vgamma}_T)$ are the elements of the residual vector $\widehat{\vu}_t =\vy_t -\widehat{\tau}_{g,T}\,t^{\widehat{\theta}_T}\vones_N-\mZ_t'\widehat{\vbeta}_T^{}$. The remaining elements are $\diff \vx_t=\vv_t$. \begin{assumption} \begin{enumerate}[(a)] • • $k(0)=1$, $k(\cdot)$ is continuous at zero, and $\sup_{x\geq 0}\left| k(x) \right| < \infty$. • $\int_0^\infty \bar{k}(x)dx<\infty$, where $\bar{k}(x)= \sup_{y\geq x} \left| k(y) \right|$. • The bandwidth parameters $\{b_T : T\geq 1 \}$ satisfies $\{b_T\}\subseteq (0,\infty)$ and $\lim_{T\to\infty} \left( b_T^{-1} + T^{-1/2} b_T \ln T \right) = 0$. \end{enumerate} \end{assumption} The conditions on the kernel function $k(\cdot)$, Assumptions (ref)(a)--(b), are identical to those in jansson2002. jansson2002 remarks that these assumptions “\emph{would appear to be satisfied by any kernel in actual use}”. Commonly used kernel functions such as the Bartlett, Parzen, and Quadratic Spectral kernels indeed satisfy all these assumptions. Assumption (ref)(c) differs from the usual requirement, $\lim_{T\to\infty} \left( b_T^{-1} + T^{-1/2} b_T \right) = 0$, by a factor $\ln T$. The difference is caused by the estimation error in $\widehat{\theta}_T$. This error causes the residuals $\{ \widehat{\vu}_t \}$ to be less close to the innovations $\{ \vu_t\}$ and we balance this by including autocovariance matrices of higher lags at a slower pace. \begin{theorem} Under Assumptions (ref)-(ref), we have $\widehat{\mDelta}_T\pto \mDelta$ and $\widehat{\mOmega}_T \pto \mOmega$. \end{theorem} \subsection{Simulation-Based Inference} The limiting distribution in Theorem (ref) is nonpivotal and thus not directly suited for inference. Some popular solutions for linear-in-parameters cointegration models are: saikkonen1992's\ (saikkonen1992) dynamic least squares, the integrated modified OLS and fixed-$b$ approaches by vogelsangwagner2014, and the fully modified approach advocated in phillipshansen1990 and phillips1995. For a nonlinear-in-parameter model as in (ref), a preliminary Monte Carlo exercise\footnote{The details are available in Section (ref) of the Supplementary Material. The analytical results in that appendix also suggest that the convergence speed of $\widehat{\theta}_T$ to $\theta$ is too slow to recover the standard zero-mean Gaussian limiting distribution.} shows poor performance for fully modified inference but promising results for a simulation based approach. We pursue the latter method for the remainder of this paper. The main idea behind the simulation based approach is to replace nuisance parameters by consistent estimates and to rely on Monte Carlo (MC) simulations to approximate the limiting distribution. The empirical quantiles of these MC draws allow us to conduct inference. Clearly, this kind of approach will provide exact inference when the limiting distribution is invariant with respect to the nuisance parameters (e.g. dufourkhalaf2002 and dufour2006). In the absence of such invariance, wangwuzhu2018 and bergamellibianchikhalafurga2019 show that the simulation approach remains asymptotically justified in several model specification. We adapt the algorithm from wangwuzhu2018 to the current setting and prove its asymptotic validity. \begin{algorithm}[Simulation-Based Inference] \ \begin{enumerate}[\textsc{Step} 1:] • Estimate $\widehat{\vgamma}_T$ and use the residuals $\{ \widehat{\vu}_t \}$ to compute the estimators $\widehat{\mDelta}_T$ and $\widehat{\mOmega}_T$ from (ref). • Repeat for $j=1,\ldots,J$, \end{enumerate} \begin{enumerate}[(a)] • Draw random variables $\{ \ve_t \}_{t=1}^{M_T}$ i.i.d. from $\rN(\vzeros,\mI_{2N})$. • Compute $\left[ \begin{smallmatrix} \widehat{\vmu}_t \\ \widehat{\vv}_t \end{smallmatrix} \right]= \widehat{\mOmega}_T^{1/2} \ve_t $ and the partial sum $\widehat{\vchi}_t = \big[\widehat{\chi}_{1,t},\ldots,\widehat{\chi}_{N,t} \big]'= \sum_{s=1}^t \widehat{\vupsilon}_s$. • Let $\widehat{\mJ}\big(t;\widehat{\vgamma}_{T}\big)=\Big[\,\widehat{\tau}_{g,T}\,t^{\widehat{\theta}_T}\ln{t}\,\vones_N,t^{\widehat{\theta}_T}\vones_N,\widehat{\calZ}_t'\Big]'$, where $\widehat{\calZ}_t=\diag\big[\widehat{\vz}_{1,t},\ldots,\widehat{\vz}_{N,t}\big]$ with $\widehat{\vz}_{i,t}=\Big[1,t,\widehat{\chi}_{i,t},\ldots,\widehat{\chi}_{i,t}^{\,p_i}\Big]'$ (for $i=1,\ldots,N$). For a given $M_T$, construct the $j$\textsuperscript{th} simulated draw as \begin{equation*} \widehat{{\bm \calJ}}^{(j)} \left(\widehat{\vgamma}_T,\widehat{\mOmega}_T,\widehat{\mDelta}_{vu}^-\right) =\left\{\text{ $ \mG_{\widehat{\vgamma}_T,M_T}^{\prime-1}\left[\sum_{t=1}^{M_T}\widehat{\mJ}\big(t;\widehat{\vgamma}_{T}\big)\widehat{\mJ}\big(t;\widehat{\vgamma}_{T}\big)'\right]\mG_{\widehat{\vgamma}_T,M_T}^{-1} $} \right\}^{-1} \left\{\text{ $ \mG_{\widehat{\vgamma}_T,M_T}^{\prime-1}\left[\sum_{t=1}^{M_T}\widehat{\mJ}\big(t;\widehat{\vgamma}_{T}\big)\,\widehat{\vmu}_t\right]+\widehat{\bm{\calB}}_{vu}^{-} $} \right\}, \end{equation*} where $\widehat{\mDelta}_{vu}^-$ is a consistent estimator of the lower-left subblock of $\mDelta^{-}=\left[\begin{smallmatrix} \mDelta_{uu}^{-} & \mDelta_{uv}^{-}\\ \mDelta_{vu}^{-} & \mDelta_{vv}^{-} \end{smallmatrix}\right]=\mSigma-\mDelta'$, and $\widehat{\bm{\calB}}_{vu}^{-}=\left[\vzeros_{1\times 2},\widehat{\vb}_{1}' \widehat{\mDelta}_{v_1u_1}^{-},\dots,\widehat{\vb}_{N}'\widehat{\mDelta}_{v_Nu_N}^{-}\right]'$ with $\widehat{\vb}_{i}=\Bigg[\vzeros_{1\times 2},1,2\frac{1}{M_T}\sum_{t=1}^{M_T}\Big(\frac{\widehat{\chi}_{i,t}}{\sqrt{M_T}}\Big),\ldots,p_i\frac{1}{M_T}\sum_{t=1}^{M_T}\Big(\frac{\widehat{\chi}_{i,t}}{\sqrt{M_T}}\Big)^{p_i-1}\Bigg]'$. \end{enumerate} \begin{enumerate}[\textsc{Step} 1:] \setcounter{enumi}{2} • Use the empirical quantiles of elements of $\left\{\widehat{{\bm \calJ}}^{(1)},\ldots, \widehat{{\bm \calJ}}^{(J)} \right\}$ to conduct inference. \end{enumerate} \end{algorithm} Algorithm 1 uses a discretisation in $M_T$ steps to approximate the limiting distribution of the parameters. In practice, and in accordance with Theorem (ref), we can take $M_T=T$. Remark (ref) details how simulation-based inference can be used to test hypotheses concerning the model's parameters. Discussions on size and power are also presented there. \begin{theorem} Suppose Assumptions (ref)-(ref) hold, let $\{M_T\}\subseteq (0,\infty)$ with $\lim_{T\to\infty} \frac{M_T}{T} \leq \kappa$, for some $\kappa<\infty$, then we have \begin{equation} \begin{aligned} &\left\{\mG_{\widehat{\vgamma}_T,T}^{\prime-1}\left[\sum_{t=1}^{M_T}\widehat{\mJ}\big(t;\widehat{\vgamma}_{T}\big)\widehat{\mJ}\big(t;\widehat{\vgamma}_{T}\big)'\right]\mG_{\widehat{\vgamma}_T,T}^{-1} \right\}^{-1} \left\{ \mG_{\widehat{\vgamma}_T,T}^{\prime-1}\left[\sum_{t=1}^{M_T}\widehat{\mJ}\big(t;\widehat{\vgamma}_{T}\big)\,\widehat{\vmu}_t\right]+\widehat{\bm{\calB}}_{vu}^{-} \right\} \\ &\qquad\qquad\qquad\qquad\longrightarrow_{d^*} \left(\int \mJ(r;\vgamma_0)\mJ(r;\vgamma_0)' dr\right)^{-1}\left(\int \mJ(r;\vgamma_0) d\bm{B}_{u}(r)+{\bm{\mathcal{B}}}_{vu}\right), \end{aligned} \end{equation} in probability. \end{theorem} Theorem (ref) establishes the asymptotic validity of the simulation approach. That is, for a large enough $J$, the empirical quantiles of the simulated distribution will coincide with the asymptotic distribution. Two remarks are important. First, even though the simulation algorithm is adapted from wangwuzhu2018, the proof of Theorem (ref) is not. In particular, the method of proof is similar to Theorem (ref) and continues to allow for endogeneity of the regressors. Second, the simulation approach mimics the stochastic integrals in the limiting distribution directly. It therefore suffices to draw normally distributed random variables in Step 2(a) and use consistent long-run covariance estimates to replicate the covariance structure of the underlying Brownian motions. Compared to a bootstrap procedure, this simulation approach has the advantage of avoiding tedious NLS re-estimation on bootstrap samples but it comes at the cost of forsaking possible asymptotic refinements. \begin{remark} Step 3 in Algorithm 1 has been kept general for notational convenience. An illustrative example is as follows. Assume we are interested in $H_0:\phi_{2,1}=0$ (irrelevance of the regressor $x_{1,t}^2$) when $$ y_{i,t} = \tau_g t^{\theta} + \tau_{1,i} + \tau_{2t,i} t + \phi_{1,i} x_{i,t} + \phi_{2,i} x_{i,t}^2 + u_{i,t}, \qquad i=1,\ldots,N,\quad t=1,\ldots,T. $$ Under $H_0$, we have $T^{3/2}\widehat{\phi}_{2,1} = \ve_6' \mG_{\vgamma_0,T}\big(\widehat{\vgamma}_T-\vgamma_0\big) \longrightarrow_{d} \ve_6' {\bm \calJ}(\vgamma_0)$ with $\ve_k$ being the k\textsuperscript{th} basis vector in $\SR^{2+p}$. Denoting the empirical $\zeta$-quantiles of $\left\{\ve_6'\widehat{{\bm \calJ}}^{(1)},\ldots, \ve_6'\widehat{{\bm \calJ}}^{(J)} \right\}$ by $c_\zeta$, a test of size $\alpha$ will reject for $T^{3/2}\widehat{\phi}_{2,1}<c_{\alpha/2}$ or $T^{3/2}\widehat{\phi}_{2,1}>c_{1-\alpha/2}$. Under the alternative $\phi_{2,1} \neq 0$, we rewrite the test statistic as $T^{3/2}\widehat{\phi}_{2,1}=T^{3/2} \big(\, \widehat{\phi}_{2,1} - \phi_{2,1}\big) + T^{3/2} \phi_{2,1}$. Statistical power is guaranteed because the simulation approach mimics the asymptotic distribution and is thus bounded, whereas the second term diverges. \end{remark} \subsection{KPSS-Type Test for the Null of Cointegration} The correct specification of the nonlinear cointegrating relation will result in a stationary error process $\{u_t\}_{t\in\SZ}$. We consider a KPSS-type test statistic for the null of stationarity. The candidate statistic is $\widetilde K_T^+ = \frac{1}{T^2}\sum_{t=1}^T \left\|\widehat{\mOmega}_{u.v}^{-1/2}\sum_{i=\ell}^{t}\widehat{\vu}_i^+\right\|^2$, where $\widehat{\vu}_t^+=\vy_t - \widehat{\mOmega}_{uv}\widehat{\mOmega}_{vv}^{-1} \diff \vx_t -\widehat{\tau}_{g,T}\,t^{\widehat{\theta}_T}\vones_N-\mZ_t'\widehat{\vbeta}_T$ and $\widehat{\mOmega}_{u.v}$ is a consistent estimator of $\mOmega_{u.v}=\mOmega_{uu}-\mOmega_{uv}\mOmega_{vv}^{-1} \mOmega_{vu}$. This statistic is stochastically bounded under the null hypothesis but diverges under the alternative. Rejections of the null hypothesis are an indication of a spurious relationship and/or an incorrect functional form of the nonlinear cointegrating relationship. Several authors have reported model settings in which the asymptotic null distribution of $K_T^+$ is known, e.g. kwiatkowskiphillipsschmidtshin1992 and wagnerhong2016. The estimation of $\vtheta$ contaminates the limiting distribution of $\widetilde K_T^+$ with nuisance parameters.\footnote{Proposition 5 in wagnerhong2016 shows that the limiting distribution of $K_T^+$ is free of nuisance parameters if $\vtheta_0$ is known and only a single integrated regressor occurs with integer powers greater than one. This result does not carry over to the current setting because of the estimation error in $\widehat{\vtheta}_T$.} choisaikkonen2010, wagnerhong2016, jianglupark2019, and linreuvers2019, have shown that subsampling can resolve this issue. We will follow their approach and use subsamples of size $q_T$ to compute the test statistics. \begin{theorem} Under Assumptions (ref)-(ref) and if $\lim_{T\to\infty} \left( q_T^{-1} + (\ln T) \left(\frac{q_T}{T} \right)^{\theta_L+\frac{1}{2}} \right) = 0$, then for any $\ell\in\{1,\ldots,T-q_T+1 \}$, we have \begin{equation} K_{q_T,\ell}^{+}=\frac{1}{q_T}\sum_{t=\ell}^{\ell+q_T-1}\left\|\frac{1}{\sqrt{q_T}}\widehat{\mOmega}_{u.v}^{-1/2}\sum_{i=\ell}^{t}\widehat{\vu}_i^+\right\|^2\longrightarrow_{d} \int \left\| \mW(r) \right\|^2dr, \end{equation} where $\bm W(\cdot)$ denotes an $N$-dimensional standard Brownian motion. \end{theorem} Theorem (ref) does not provide any guidance on the choices for the starting value $\ell$ and the subsample size $q_T$. First, for a given $q_T$, choisaikkonen2010 argue that the use of a single subsample (instead of all $T$ observations) implies a significant loss of power. We follow their example and combine all $M=[T/q_T]$ subresidual series of length $q_T$ using a Bonferroni procedure. That is, we create subresiduals series by selecting adjacent blocks of $q_T$ residuals while alternating between the start and end of the sample. We calculate the KPSS-type test statistic for each subseries, say $K_1,\ldots, K_M$, and reject the null of stationarity at significance $\alpha$ whenever $\max\{K_1,\ldots,K_M \}$ exceeds $c_{\alpha/M}$ which is defined by $\mathbb{P}\left( \int \big\| \bm W(r) \big\|^2 dr \geq c_{\alpha/M} \right) = \alpha/M$ . Finally, we select the block size $q_T$ using romanowolf2001's\ (romanowolf2001) minimum volatility rule. The approach is now completely data-driven. \section{Simulations} This section lists various Monte Carlo simulations showing that the asymptotic approximations from Section (ref) provide useful guidance in finite samples. Further details on the implementation are as follows. The long-run covariance matrices in (ref) are computed using the Barlett kernel, $k(x)= 1 -|x|$ for $|x|\leq 1$ (and zero otherwise), and the bandwidth selection method described in andrews1991. Simulated limiting distributions are based on $J=299$ replicates and we set $M_T=T$. We test at 5% significance and report results based on $3,000$ Monte Carlo replications. Two data generating processes are studied: DGP1 and DGP2. \subsubsection*{DGP1: Empirical size and power of the coefficient tests} This DGP is inspired by the simulation study in wagnergrabarczykhong2019. It augments their quadratic seemingly unrelated cointegrating polynomial regression model with a global flexible trend. That is, we consider \begin{equation} y_{i,t} = \tau_g t^\theta + \tau_{1,i} + \tau_{2,i} t + \phi_{1,i} x_{i,t} + \phi_{2,i} x_{i,t}^2 + u_{i,t}, \qquad i=1,\ldots,N,\qquad t=1,\ldots,T, \end{equation} and compute $u_{i,t}$ and $\diff x_{i,t}=v_{i,t}$ recursively as \begin{equation*} u_{i,t} = \rho_1 u_{i,t-1} + \varepsilon_{i,t} + \rho_2 e_{i,t}, \qquad v_{i,t} = e_{i,t} + 0.5 e_{i,t-1}. \end{equation*} All recursions are initialized from zero, i.e. $x_{i,0}=u_{i,0}=e_{i,0}=0$ ($i=1,\ldots,N$). The innovations $\vepsi_t=[\varepsilon_{1,t},\ldots,\varepsilon_{N,t}]'$ and $\ve_t=[e_{1,t},\ldots,e_{N,t}]'$ are drawn independently as $\vepsi_t\stackrel{i.i.d.}{\sim}\rN(\vzeros,\mSigma_{\varepsilon\varepsilon})$ and $\ve_t\stackrel{i.i.d.}{\sim}\rN(\vzeros,\mSigma_{ee})$, where $$ \mSigma_{\varepsilon\varepsilon} = \begin{bmatrix} 1 & \rho_3 & \cdots & \rho_3\\ \rho_3 & 1 & \cdots & \rho_3\\ \vdots & \vdots & \ddots & \vdots\\ \rho_3 & \rho_3 & \cdots & 1\\ \end{bmatrix} ,\qquad\text{and}\qquad \mSigma_{ee} = \begin{bmatrix} 1 & \rho_4 & \cdots & \rho_4\\ \rho_4 & 1 & \cdots & \rho_4\\ \vdots & \vdots & \ddots & \vdots\\ \rho_4 & \rho_4 & \cdots & 1\\ \end{bmatrix}. $$ Regarding the global trend in (ref), we set $\tau_g=-0.2$ and consider $\theta\in\{0.8,1.3,1.8\}$. All other coefficient values are inspired by wagnergrabarczykhong2019. That is, $\tau_{1,i}=1$, $\tau_{2,i}=1$ and $\phi_{1,i}=5$ are identical across equations.\footnote{This homogenous parametrisation is particularly convenient to study the impact of the cross-sectional dimension. That is, we can vary $N$ without having to provide additional parameter values. DGP2 is directly inspired by the empirical application and thus more realistic.} Also, we let $\rho_1=\rho_2=\rho_3=\rho_4$ and redefine these four parameters as $\rho$. We vary $\rho\in\{0,0.3,0.6,0.8 \}$, $N\in\{3,5,10\}$, and $T\in\{150,300,600\}$. In line with the typical EKC application, we test for the significance of $x_{i,t}^2$. We set $\phi_{2,i}=0$ for $i=1,\ldots,N$ and report the empirical size of the single equation test for $H_0: \phi_{2,1}=0$ and the joint test for $H_0: \phi_{2,1}=\ldots=\phi_{2,N}=0$. \begin{table}[t] \caption{The empirical size (in %) of the single-equation tests $H_0:\phi_{2,1}=0$ and the joint test for $H_0: \phi_{2,1}=\ldots=\phi_{2,N}=0$ with $\phi_{2,i}$ denoting the coefficient in front of $x_{i,t}^2$. The Monte Carlo results are based on: simulated inference with $\theta$ estimated by NLS (SimNLS), simulated inference with known $\theta=1.3$ (SimNLS($\theta_0$)), and two Fully Modified estimators for systems developed by wagnergrabarczykhong2019 with known $\theta=1.3$ (FM-SOLS($\theta_0$) and FM-SUR($\theta_0$)).} \resizebox{\textwidth}{!}{ \begin{tabular}{crrrrlrrrrlrrrr} \toprule $\theta_0=1.3$ & \multicolumn{4}{c}{$N=3$} & \multicolumn{1}{c} & \multicolumn{4}{c}{$N=5$} & \multicolumn{1}{c} & \multicolumn{4}{c}{$N=10$} \\ \midrule $\rho$ & \multicolumn{1}{c}{SimNLS} & \multicolumn{1}{c}{SimNLS($\theta_0$)} & \multicolumn{1}{c}{FM-SOLS($\theta_0$)} & \multicolumn{1}{c}{FM-SUR($\theta_0$)} & & \multicolumn{1}{c}{SimNLS} & \multicolumn{1}{c}{SimNLS($\theta_0$)} & \multicolumn{1}{c}{FM-SOLS($\theta_0$)} & \multicolumn{1}{c}{FM-SUR($\theta_0$)} & & \multicolumn{1}{c}{SimNLS} & \multicolumn{1}{c}{SimNLS($\theta_0$)} & \multicolumn{1}{c}{FM-SOLS($\theta_0$)} & \multicolumn{1}{c}{FM-SUR($\theta_0$)} \\ \midrule \multicolumn{15}{l}{\textbf{Panel A: Single-equation test}} \\ \midrule \multicolumn{15}{l}{$T=150$} \\ \midrule 0 & 4.77 & 4.63 & 9.53 & 10.77 & & 4.43 & 4.90 & 10.03 & 12.80 & & 4.37 & 4.47 & 9.70 & 16.63 \\ 0.3 & 4.80 & 4.77 & 10.47 & 11.97 & & 4.53 & 4.63 & 10.00 & 13.00 & & 4.23 & 4.27 & 11.60 & 19.00 \\ 0.6 & 4.83 & 4.53 & 11.90 & 12.93 & & 3.93 & 4.17 & 11.87 & 16.50 & & 5.13 & 4.90 & 14.23 & 31.70 \\ 0.8 & 4.50 & 4.67 & 14.00 & 19.10 & & 4.90 & 4.60 & 15.23 & 26.53 & & 5.27 & 4.43 & 16.83 & 56.47 \\ \midrule \multicolumn{15}{l}{$T=300$} \\ \midrule 0 & 4.43 & 4.10 & 8.03 & 8.63 & & 4.43 & 4.20 & 7.33 & 8.33 & & 4.50 & 4.67 & 8.60 & 11.73 \\ 0.3 & 4.20 & 4.60 & 8.13 & 9.20 & & 4.37 & 4.77 & 8.97 & 9.80 & & 4.43 & 4.23 & 8.67 & 12.83 \\ 0.6 & 5.80 & 5.70 & 10.23 & 11.80 & & 4.97 & 4.97 & 10.17 & 12.73 & & 4.40 & 4.30 & 10.83 & 18.77 \\ 0.8 & 4.97 & 4.43 & 11.17 & 13.40 & & 4.60 & 4.53 & 12.07 & 19.37 & & 4.03 & 3.53 & 13.67 & 36.53 \\ \midrule \multicolumn{15}{l}{$T=600$} \\ \midrule 0 & 4.67 & 4.67 & 7.00 & 7.20 & & 4.77 & 4.73 & 6.57 & 7.53 & & 4.27 & 4.23 & 7.37 & 9.20 \\ 0.3 & 4.77 & 4.37 & 7.67 & 7.77 & & 4.43 & 4.73 & 7.37 & 7.40 & & 4.80 & 4.67 & 7.90 & 9.57 \\ 0.6 & 5.60 & 5.07 & 9.43 & 9.50 & & 5.43 & 5.43 & 8.77 & 10.30 & & 5.30 & 4.83 & 9.57 & 14.10 \\ 0.8 & 4.70 & 4.50 & 8.73 & 9.20 & & 5.10 & 5.23 & 9.27 & 13.87 & & 5.80 & 5.20 & 11.47 & 24.63 \\ \midrule \multicolumn{15}{l}{\textbf{Panel B: Joint test}} \\ \midrule \multicolumn{15}{l}{$T=150$} \\ \midrule 0 & 4.23 & 4.03 & 12.70 & 15.10 & & 4.10 & 4.07 & 14.50 & 21.57 & & 3.63 & 3.47 & 25.67 & 51.23 \\ 0.3 & 4.80 & 4.57 & 14.37 & 17.33 & & 3.80 & 3.70 & 19.50 & 26.90 & & 3.70 & 3.77 & 31.63 & 59.73 \\ 0.6 & 4.27 & 4.03 & 18.03 & 22.07 & & 3.80 & 3.57 & 23.67 & 37.77 & & 3.13 & 3.03 & 40.27 & 82.17 \\ 0.8 & 3.13 & 2.93 & 23.37 & 30.47 & & 3.33 & 2.73 & 31.60 & 57.20 & & 2.00 & 1.50 & 49.73 & 82.60 \\ \midrule \multicolumn{15}{l}{$T=300$} \\ \midrule 0 & 4.80 & 4.73 & 10.07 & 10.87 & & 4.30 & 4.37 & 10.63 & 14.83 & & 3.77 & 3.80 & 18.17 & 31.80 \\ 0.3 & 4.77 & 4.97 & 11.60 & 12.93 & & 4.67 & 4.67 & 14.03 & 17.70 & & 3.40 & 3.47 & 19.13 & 36.53 \\ 0.6 & 4.93 & 4.40 & 14.33 & 14.87 & & 3.77 & 3.90 & 17.50 & 25.63 & & 3.10 & 3.10 & 29.20 & 59.97 \\ 0.8 & 3.87 & 3.40 & 17.40 & 20.97 & & 3.33 & 2.77 & 22.97 & 38.80 & & 2.60 & 2.13 & 37.87 & 86.33 \\ \midrule \multicolumn{15}{l}{$T=600$} \\ \midrule 0 & 4.27 & 4.53 & 7.37 & 8.03 & & 4.33 & 4.13 & 8.50 & 10.83 & & 3.97 & 4.07 & 12.90 & 19.23 \\ 0.3 & 4.93 & 5.17 & 9.23 & 10.00 & & 4.80 & 4.53 & 10.57 & 12.30 & & 4.60 & 4.50 & 14.87 & 24.03 \\ 0.6 & 4.07 & 3.80 & 11.63 & 12.43 & & 4.57 & 4.73 & 12.30 & 16.97 & & 4.30 & 4.23 & 21.73 & 38.80 \\ 0.8 & 5.00 & 4.53 & 12.77 & 14.37 & & 3.57 & 3.80 & 15.67 & 25.33 & & 3.60 & 3.57 & 26.43 & 66.63\\ \bottomrule \end{tabular} } \end{table} For $\theta_0=1.3$, the empirical size of various tests are displayed in Table (ref).\footnote{The results for $\theta=0.8$ and $\theta=1.8$ are qualitatively the same. For brevity, we do not include these results in the main paper. The interested reader can find such simulation results in Section (ref) of the Supplementary Material.} These tests are based on four estimators: (1) the NLS estimator with simulated critical values as described in Section (ref) (SimNLS); (2) the NLS estimator with simulated critical values and the true value for $\theta_0=1.3$ being provided (SimNLS($\theta_0$)); (3) the FM-SOLS estimator based on $\theta_0=1.3$ (FM-SOLS($\theta_0$)); and (4) the FM-SUR estimator based on $\theta_0=1.3$ (FM-SUR($\theta_0$)). The main findings are as follows: \begin{enumerate}[(a)] • The simulation-based approaches SimNLS and SimNLS($\theta_0$) offer better size control. The size improvements are particularly pronounced when $T=150$ and $\rho=0.8$. The differences in the empirical size of SimNLS and SimNLS($\theta_0$) are small. • Size distortions are more severe when $N$ increases and/or a joint test is performed. The same observation was made in wagnergrabarczykhong2019. The behaviour of the simulation-based and fully modified tests is opposite in these cases. SimNLS and SimNLS($\theta_0$) tend to become conservative whereas FM-SOLS and FM-SUR are oversized. \end{enumerate} We subsequently simulate power curves.\footnote{Power curves are computationally more intensive. We economize computational time by (1) reducing the number of Monte Carlo replicates to 1,000 and (2) investigating a subset of all possible parameter configurations.} The specification of the single equation test and joint test are as before but we now vary $\phi_{2,1}=\ldots=\phi_{2,N}$ over the set $[-0.008,-0.007,\ldots, 0]$. We take $\rho=0.3$, $\theta=1.3$ and $N=3$ as the baseline scenario and subsequently vary these quantities one-by-one. Figures (ref)--(ref) show the results. As expected, power increases with increasing sample size, and as $\phi_{2,i}$ moves away from zero. \begin{figure}[h] \begin{subfigure}{.9\textwidth} \caption \end{subfigure} \\ \begin{subfigure}{.6\textwidth} \caption \end{subfigure} \\ \begin{subfigure}{.6\textwidth} \caption \end{subfigure} \caption{The power curves for the \emph{single equation test} $H_0:\phi_{2,1}=0$. The reference model is DGP1 with $\rho=0.3$, $\theta=1.3$ and $N=3$. We vary the parameters of this reference specification one-by-one while keeping the remaining two parameters fixed at their baseline values. Specifically, we study changes in: \textbf{(a)} the serial correlation and endogeneity parameter $\rho$, \textbf{(b)} the nonlinear deterministic time trend power, and \textbf{(c)} the cross-sectional dimension.} \end{figure} \begin{figure}[h] \begin{subfigure}{.9\textwidth} \caption \end{subfigure} \\ \begin{subfigure}{.6\textwidth} \caption \end{subfigure} \\ \begin{subfigure}{.6\textwidth} \caption \end{subfigure} \caption{The power curves for the \emph{joint test} for $H_0: \phi_{2,1}=\ldots=\phi_{2,N}=0$. The reference model is DGP1 with $\rho=0.3$, $\theta=1.3$ and $N=3$. We vary the parameters of this reference specification one-by-one while keeping the remaining two parameters fixed at their baseline values. Specifically, we study changes in: \textbf{(a)} the serial correlation and endogeneity parameter $\rho$, \textbf{(b)} the nonlinear deterministic time trend power, and \textbf{(c)} the cross-sectional dimension.} \end{figure} \subsubsection*{DGP2: Illustrative simulations in line with the empirical application} Our second set of simulations is tailored towards the empirical application. That is, we employ parametrizations that mimic the distributional properties of data. Generally speaking, we first estimate the baseline model specification on the data and subsequently fit a VAR(1) specification on the stacked vector of residuals and first-differenced explanatory variables.\footnote{All details on the simulation designs for DGP2(a)--(c) are available in Section (ref) of the Supplementary Material.} In line with the empirical application, these simulations use $N=6$ and $T=145$. All results are displayed in Figure (ref). Below, we motivate the simulation settings in view of the EKC application and draw conclusions. \begin{enumerate}[(a)] • \textbf{Correctly specified model}: The specification $ y_{i,t} = \tau_g t^{\theta} + \tau_{1,i} + \tau_{2,i} t + \phi_{1,i} x_{i,t} + \phi_{2,i} x_{i,t}^2 + u_{i,t}$ with $\phi_{2,i}=0$ is estimated on the data. We subsequently move $\phi_{2,1}=\ldots=\phi_{2,6}$ away from zero in the DGP and check whether we can detect the resulting curvature caused by the integrated variable. Power curves for the individual and joint test for the coefficients in front of $x_{i,t}^2$ are found in Figures (ref)(a) and (ref)(b), respectively. Clearly, nonlinear effects due to $x_{i,t}^2$ are detectable. The statistical power varies across units because (contrary to DGP1) time series properties are now heterogenous across equations. • \textbf{Redundant global trend}: Assumption (ref) requires the global trend to be relevant. This simulation DGP investigates how violations of this assumption affect the typical EKC coefficient test. We obtain parameter values by fitting the model $y_{i,t}= \tau_{1,i} + \tau_{2,i} t + \phi_{1,i} x_{i,t} + \phi_{2,i} x_{i,t}^2 + u_{i,t}$ with $\phi_{2,i}=0$. As in (a), we vary $\phi_{2,1}=\ldots=\phi_{2,6}$ and test for the significance of these parameters. The solid lines in Figures (ref)(c) and (ref)(d) are power curves obtained using the correctly specified DGP whereas markers indicate the power when a redundant global trend is estimated as well. The redundant trend has virtually no influence on the statistical power of the coefficient tests of the first five series. There is a power loss for $i=6$. An inspection of the coefficients offers an explanation. The estimated coefficients in front of the global trend are mostly small ($10^{-10}$ to $10^{-9}$) and thus irrelevant. However, in a fraction of cases the flexible trend mimics the curvature in the $6$\textsuperscript{th} series causing the quadratic stochastic trend to become insignificant. As reported in the introduction, the power of the joint test does not suffer from the inclusion of a redundant trend. • \textbf{KPSS test}: Nonstationary residuals are an indication of model misspecification. That is, either the regression is spurious or the functional form of the cointegrating relation is misspecified. We look at the latter situation. The simulation DGP is the quadratic GCPR as in DGP2(a) but the quadratic component is missing in the fitted model. The empirical rejection frequency of the KPSS test (Figure (ref)) is signalling that there are specification issues. However, a comparison with Figures (ref)(a) and (ref)(b) also reveals that if the source of misspecification is known, then a dedicated coefficient test leads to higher power. \end{enumerate} \begin{figure}[h] \begin{subfigure}{.5\textwidth} \caption \end{subfigure} \begin{subfigure}{.5\textwidth} \caption \end{subfigure} \begin{subfigure}{.5\textwidth} \caption \end{subfigure} \begin{subfigure}{.5\textwidth} \caption \end{subfigure} \center \begin{subfigure}{.5\textwidth} \caption \end{subfigure} \caption{An overview of various power curves. \textbf{(a)} Unit-specific power curves when testing $H_0 : \phi_{2,i}=0$ ($i=1,\ldots,6$) for a correctly specified model. \textbf{(b)} The power curve when testing $H_0: \phi_{2,1}=\ldots=\phi_{2,6}=0$ for a correctly specified model. \textbf{(c)} The empirical rejection frequencies for a correctly specified model (lines) and an estimation with a redundant global time trend (dots). Individual coefficients are tested. \textbf{(d)} As in (c), but now for the joint test $H_0: \phi_{2,1}=\ldots=\phi_{2,6}=0$. \textbf{(e)} The empirical power of the KPSS for a misspecified linear cointegrating relation.} \end{figure} \section{Empirical Application} We examine the evidence for an EKC for a collection of 18 countries over the period 1870--2014 ($T=145$). Economic growth is measured by GDP and we use carbon dioxide (CO\textsubscript{2}) emissions as a proxy for air pollution. The origin of these data is as follows. We use population and GDP data from the Maddison Project (see https://www.rug.nl/ggdc/historicaldevelopment/maddison/). Our carbon dioxide observations are fossil-fuel CO\textsubscript{2} emissions as made available by the Carbon Dioxide Information Analysis Center (CDIAC, see https://cdiac.ess-dive.lbl.gov). The CDIAC database ceased operation in 2017 causing these time series to be available until 2014. Both GDP and CO\textsubscript{2} emissions are expressed per capita and subsequently log-transformed. In accordance with the notation of this paper, we will denote them by $x_{i,t}$ and $y_{i,t}$, respectively. The same data (or subsets thereof) have also been studied by wagner2015, chanwang2015, wangwuzhu2018, wagnergrabarczykhong2019, and linreuvers2019.\footnote{The stationarity properties of the series have been extensively studied and discussed in these papers. We will not repeat this analysis but refer the interested reader to Section (ref) of the Supplement. The exact numbers may show (minor) differences from previously reported results due to differences in: (1) the time span of the data, (2) the implemented long-run covariance estimator, and (3) the scaling of the data. Related to scaling, we follow the official guidelines and multiply by 3.667 and $10^3$ to convert thousand of metric tons of carbon into units of carbon dioxide. Since the data will be expressed in logarithms, this rescaling effectively amounts to a change of intercept.} This conveniently allows us to compare results. All user choices (kernel specification, bandwidth selection, etc.) are kept the same as during the simulation study (see page (ref)). \subsection{An Illustration using Belgian Data} Prior to the analysis of a multivariate specification, we will first discuss several features of the individual time series (hence omitting subscripts “$i$”). The example throughout this narrative is Belgium (Figure (ref)).\footnote{The data for Austria, Belgium, and Finland are mentioned in both wagner2015 and wagnergrabarczykhong2019 to behave in line with the EKC. We discuss Belgium in the main text but the interested reader can find the same figures for Austria and Finland in Section (ref). Qualitatively, the findings for these other two countries are the same.} An inverted U-shaped relationship between GDP and CO\textsubscript{2} (both in log per capita) is clearly visible in Figure (ref)(a) and behavior like this has triggered research on the Environmental Kuznets Curve. However, the time heat map also shows that time is almost monotonically increasing along the curve. Time effects -- e.g. increasing global environmental awareness, worldwide advances in sustainable technologies -- can be valid alternative explanations for these nonlinearities and their omission can (falsely) exaggerate the influence of GDP. It is for this reason that we develop and analyse the Generalized Cointegrating Polynomial Regression (GCPR). \begin{figure}[h] \begin{subfigure}{.5\textwidth} \caption \end{subfigure} \begin{subfigure}{.5\textwidth} \caption \end{subfigure} \begin{subfigure}{.5\textwidth} \caption \end{subfigure} \begin{subfigure}{.5\textwidth} \caption \end{subfigure} \begin{subfigure}{.5\textwidth} \caption \end{subfigure} \begin{subfigure}{.5\textwidth} \caption \end{subfigure} \caption{Overview graphs for Belgium over 1870-2014. \textbf{(a)} $\log(\text{GDP})$ versus $\log(\text{CO}_2)$ (both per capita). \textbf{(b)} The same series as in subfigure (a), but now using detrended variables. \textbf{(c)} The log per capita CO\textsubscript{2} emissions time series for Belgium over time. \textbf{(d)} The residual sum of squares (RSS) for the nonlinear model specification $y_t=\tau_1+\tau_2 t + \phi_1 x_t+ \phi_2 x_t^\theta+u_t$ for various values of $\theta$. \textbf{(e)} The RSS as a function of $\theta$ for the flexible nonlinear trend specification $y_t=\tau_1 + \tau_2 t + \tau_3 t^\theta + \phi x_t+u_t$. \textbf{(f)} The relation between $x_t$ and $y_t$ after partialling out the constant, linear trend, and flexible deterministic trend.} \end{figure} More evidence for the importance of time effects is available in Figure (ref)(b). This figure depicts the same per capita series after detrending.\footnote{The perronyabu2009 test allows us to test for the presence of a deterministic trend irrespectively of the series being trend-stationary or having an unit root. The results of this test (see supplement) indicate that log per capita GDP is likely to have a deterministic trend component. It is thus recommended to have a deterministic trend in the model for log per capita $\mathrm{CO}_2$ emissions and the visual inspection of the relationship between GDP and $\mathrm{CO}_2$ emissions (in log per capita) should take place after partialling out this deterministic trend.} The inverted U-shape is now (visually) less pronounced or even absent. Finally, let us depart from a traditional linear cointegration specification: $y_{t} = \tau_1 + \tau_2 t + \phi_1 x_{t}+u_{t}$. This model cannot incorporate any nonlinear behaviour over time and is therefore ill-suited to fit the data displayed in Figure (ref)(c). Cointegrating polynomial regressions use integer powers of $x_{t}$ to describe the curvature over time. More general, as in huphillipswang2019, we can allow for an integrated regressor with a flexible power and estimate $y_{t}=\tau_1+\tau_2 t + \phi_1 x_{t}+ \phi_2 x_{t}^\theta+u_{t}$. The residual sum of squares (RSS) of the NLS estimator for this specification is shown in Figure (ref)(d). The absence of a minimum at $\theta=2$ casts doubt on the commonly used quadratic specification in $x_{t}$. Additionally, the lack of any minimum might be interpreted as a sign that log per capita GDP is not the source of nonlinearity. This finding is not specific for Belgium. There are no minima in the RSS for 15 out of 18 countries (see Section (ref)). For the remaining three countries -- Denmark, France and the Netherlands -- minima are found at $\widehat\theta_{DK}=1.46$, $\widehat\theta_{FR}=3.61$ and $\widehat\theta_{NL}=1.28$, respectively. Alternatively, we can describe the nonlinearity in the data using a flexible deterministic trend as in $y_{t}=\tau_1 + \tau_2 t + \tau_3 t^\theta + \phi_1 x_{t}+u_{t}$. The RSS in Figure (ref)(e) now exhibits a clear minimum. Further empirical analysis on individual countries (see Section (ref) of the Supplement) suggests that: (1) the inclusion of a flexible time trend renders all quadratic effects in squared log per capita GDP insignificant, and (2) models remain well-specified after removing quadratic income effects from the model. These results suggest -- albeit in a univariate setting -- that flexible time trends gives a more satisfactory (or at least competing) description of the nonlinearities in the data. \subsection{Seemingly Unrelated Regression} The interpretation of a country-specific flexible deterministic trend is complicated because of its high collinearity with GDP per capita. The multivariate analysis of this section allows us to separate country-specific environmental improvements caused by national income growth from global environmental improvements. We study the following six countries ($N=6$): Austria, Belgium, Finland, the Netherlands, Switzerland, and the UK. The motivation behind this choice is as follows. First, based on data series to ours, piaggiopadilla2012, mazzantimusolesi2013, and wagnergrabarczykhong2019 report considerable evidence of parameter heterogeneity across countries.\footnote{Parameter heterogeneity is also reported for other data sets. Examples are listgallet1999, cole2005, and dijkgraafvollebergh2005.} The evidence in mazzantimusolesi2013 is anecdotal in the sense that these authors consider groups of similar countries and find different results for different groups. The lack of overlap among confidence intervals of country-specific parameters has also been interpreted as a sign of heterogeneity (section 4.2 in piaggiopadilla2012). wagnergrabarczykhong2019 explicitly test for various forms of poolability and conclude that pooling is (at most) appropriate for small subgroups of countries. This lack of parameter homogeneity justifies a multivariate approach with small $N$ rather than a panel setting. Admittedly, in the current time-series setting, studying “large $N$” is also infeasible since consistent estimators for $(2N\times 2N)$ long-run covariance matrices are required. Second, prior studies already refute the existence of a carbon dioxide EKC for several countries and little seems lost by excluding these countries from the outset.\footnote{Most of the parameters in the Generalized Cointegrating Polynomial Regression are country-specific. The estimation accuracy of these parameters should deteriorate little when focussing attention on a subset of countries. Losses will occur in the precision of the estimators for $\tau_g$ and $\theta$. There is thus a trade-off between accurate global trend estimation (improving with large $N$) and accurate LRV estimation (deteriorating with large $N$). To strike a balance and to connect to the recent literature, we continue the analysis of wagnergrabarczykhong2019 and take $N=6$.} That is, we consider the same countries as in wagnergrabarczykhong2019, who decide on these countries because their prior cointegration analysis “\emph{leads to evidence for a quadratic cointegrating EKC including a constant and linear trend}”. \begin{table}[t] \caption{Parameter estimates and test results for Models (ref)--(ref). The joint $p$-value refers to the test with null hypothesis $H_0: \phi_{2,1}=\ldots=\phi_{2,6}=0$ and is thus inapplicable for Model (ref).} \resizebox{\textwidth}{!}{ \begin{threeparttable} \begin{tabular}{c d{2.5} d{3.5} d{2.5} d{2.5} d{2.5} d{2.4} c d{2.5} d{2.2} c d{1.3} c} \toprule & \multicolumn{6}{c}{Omitted Global Trend} & & \multicolumn{5}{c}{Global Trend} \\ \cmidrule{2-13} Model & \multicolumn{6}{c}{(M1)} & & \multicolumn{2}{c}{(M2)} & & \multicolumn{2}{c}{(M3)} \\ \cmidrule{2-7}\cmidrule{9-10}\cmidrule{12-13} & \multicolumn{2}{c}{FM-SOLS} & \multicolumn{2}{c}{FM-SUR} & \multicolumn{2}{c}{SimNLS} & & \multicolumn{2}{c}{SimNLS} & & \multicolumn{2}{c}{SimNLS} \\ \cmidrule{1-13} & \multicolumn{1}{c}{$\phi_{1,i}$} & \multicolumn{1}{c}{$\phi_{2,i}$} & \multicolumn{1}{c}{$\phi_{1,i}$} & \multicolumn{1}{c}{$\phi_{2,i}$} & \multicolumn{1}{c}{$\phi_{1,i}$} & \multicolumn{1}{c}{$\phi_{2,i}$} & & \multicolumn{1}{c}{$\phi_{1,i}$} & \multicolumn{1}{c}{$\phi_{2,i}$} & & \multicolumn{2}{c}{$\phi_{1,i}$} \\ \midrule Austria & 9.37^{***} & -0.43^{***} & 3.96^{*} & -0.16 & 6.42^{***} & -0.28 & & 3.08^{***} & -0.09 & & \multicolumn{2}{c}{$1.73^{***}$} \\ Belgium & 11.78^{***} & -0.59^{***} & 9.92^{***} & -0.50^{***} & 12.36^{***} & -0.62^{**} & & 7.68^{***} & -0.36 & & \multicolumn{2}{c}{$1.01^{***} $}\\ Finland & 16.00^{***} & -0.72^{***} & 15.07^{***} & -0.68^{***} & 17.18^{***} & -0.78^{*} & & 15.19^{***} & -0.65 & & \multicolumn{2}{c}{$2.22^{***}$} \\ Netherlands & 10.68^{***} & -0.51^{***} & 9.58^{***} & -0.46^{***} & 9.27^{***} & -0.44^{*} & & 4.97^{***} & -0.20 & &\multicolumn{2}{c}{$ 1.33^{***}$} \\ Switzerland & 8.17^{***} & -0.27^{***} & 7.29^{***} & -0.23^{***} & 8.11^{***} & -0.28 & & 0.58^{*} & 0.10 & & \multicolumn{2}{c}{$2.55^{***}$} \\ UK & 9.28^{***} & -0.47^{***} & 7.93^{***} & -0.40^{***} & 9.16^{***} & -0.46^{*} & & 4.93^{***} & -0.21 & & \multicolumn{2}{c}{$1.33^{***}$} \\ \midrule Joint $p$-value & \multicolumn{2}{c}{0.00} & \multicolumn{2}{c}{0.00} & \multicolumn{2}{c}{0.16} & & \multicolumn{2}{c}{0.39} & & \multicolumn{2}{c}{---} \\ \addlinespace[0.1cm] KPSS-statistic & \multicolumn{2}{c}{3.45} & \multicolumn{2}{c}{5.10} & \multicolumn{2}{c}{3.46} & & \multicolumn{2}{c}{3.48} & & \multicolumn{2}{c}{3.78}\\ \midrule $\widehat{\tau}\, t^{\widehat{\theta}}$ & & & & & & & & \multicolumn{2}{c}{$-0.012\, t^{1.263}$} & & \multicolumn{2}{c}{$-1.374\cdot 10^{-5} t^{2.450}$} \\ \bottomrule \end{tabular} \begin{tablenotes} • Note: Asterisks denote rejection of the null hypothesis at the $^{***}1\%$, $^{**}5\%$, and $^{*}10\%$ significance level. Depending on the specific table entry, the null hypothesis refers to coefficient(s) being zero or a well-specified cointegrating relation. \end{tablenotes} \end{threeparttable} } \end{table} Having decided on the set of countries, we subsequently study the effect of the global flexible trend on EKC evidence. Table (ref) shows the estimation results of the quadratic EKC specification \begin{equation} y_{i,t} = \tau_{1,i} + \tau_{2,i} t + \phi_{1,i} x_{i,t} + \phi_{2,i} x_{i,t}^2 + u_{i,t}. \tag{M1} \end{equation} This setting (possibly with the additional constraint $\tau_{2,i}=0$) has been explored in numerous papers, for example: seldensong1994, piaggiopadilla2012, chanwang2015, wagner2015, wangwuzhu2018, and wagnergrabarczykhong2019. For Model (ref), an inverted-U relationship results when $\phi_{1,i}>0$ and $\phi_{2,i}<0$ and empirical evidence hereof is traditionally interpreted as the existence of an EKC. If these coefficients have the correct signs, then the country's turning point -- the level of economic growth at which environmental improvement starts -- can be computed as $\exp\left( -\phi_{1,i}/2\phi_{2,i} \right)$. We assess the parameter values and their significance using FM-SOLS and FM-SUR (repeating the analysis of wagnergrabarczykhong2019 for ease of comparison) and the simulated approach of Section (ref). Regardless of estimation method and country, all coefficient signs are in agreement with the EKC hypothesis. The parameters $\phi_{1,i}$ are generally significantly different from zero but the significance of $\phi_{2,i}$ does vary across estimation methods. FM-SOLS and FM-SUR typically (strongly) reject $H_0: \phi_{2,i}=0$ ($i=1,\ldots,6$) whereas evidence against these null-hypotheses is less pronounced for the simulation-based approach. The same behaviour emerges when testing $\phi_{2,1}=\ldots=\phi_{2,6}=0$ jointly. This pattern reminds of the simulation results in Table (ref) where the cross-sectional dimensions $N=5$ and $N=10$ cause over-sized tests for FM-SOLS and FM-SUR and conservative tests for simulation-based inference. The KPSS test does not indicate any signs of misspecification. Overall, Model (M1) leads to considerable evidence in favour of a quadratic cointegrating EKC. The reported evidence in favour of the EKC should not come as surprise. First, the set of countries was selected based on this criteria. Second, the visualisations of the data clearly suggest nonlinear effects (recall Figures (ref)(a) and (ref)(c) for the case of Belgium). With Model (M1) being restrictive in the sense that nonlinearities over time are solely incorporable through $x_{i,t}^2$, we expect this variable to be important. In line with our proposed Generalized Cointegrating Polynomial Regression (GCPR) framework, we subsequently add a global flexible trend and estimate \begin{equation} y_{i,t} = \tau_g t^{\theta} + \tau_{1,i} + \tau_{2,i} t + \phi_{1,i} x_{i,t} + \phi_{2,i} x_{i,t}^2 + u_{i,t}. \tag{M2} \end{equation} From a statistical perspective, the term $\tau_g t^{\theta}$ opens a different channel through which nonlinearities can be described. We refer back to the introduction for a further elaboration on this point. From an economic perspective, $\tau_g t^{\theta}$ captures changes in CO\textsubscript{2} emissions that are common across series and thus unrelated to changes in national GDPs. Parameter inference for Model (ref) is also reported in Table (ref). The contributions of $x_{i,t}^2$ are insignificant for both individual countries and all countries jointly. How about the significance of the global trend? The standard Wald test for $\tau_g=0$ is invalid because $\theta$ is unidentified under the null hypothesis (see Assumption (ref) and the related discussion). As a heuristic alternative, we vary $\theta$ over the interval $[0,2.5]$ and compute Wald statistics while assuming $\theta$ to be fixed. Comparing these Wald statistics to the 95% quantile of a $\chi^2(1)$-distributed random variable (critical value: 3.842), the range of $\theta$-values from about 0.5 to 1.75 implies a significant global trend (Figure (ref)). Having estimated $\widehat\theta = 1.263$, our analysis suggests that the global trend and not GDP per capita is the source of nonlinearity. Before interpreting this result, we first verify whether the model with $\phi_{2,1}=\ldots=\phi_{2,6}=0$ shows signs of misspecification. \begin{figure}[t] \center \caption{The magnitude of the Wald test for fixed values of $\theta$ when testing $H_0:\tau_g = 0$ under Model (ref). Dash lines display the 95% quantile of a chi-squared distributed random variable with 1 degree of freedom (red) and the NLS estimate $\widehat\theta=1.263$ for specification (ref).} \end{figure} Omitting insignificant parameters from the previous model specification, we arrive at \begin{equation} y_{i,t} = \tau_g t^{\theta} + \tau_{1,i} + \tau_{2,i} t + \phi_{1,i} x_{i,t} + u_{i,t}. \tag{M3} \end{equation} Model (ref) is linear in log per capita GDP. The positive parameter estimates for $\phi_{1,i}$ imply that \emph{at a given point in time} increases in economic growth imply increases in CO\textsubscript{2} emissions. However, as $\widehat\tau_g= -1.374\times 10^{-5}$ and $\widehat\theta= 2.45$, there will be common emission reductions over time. Also, the omission of the quadratic terms in log per capita GDP do not seem to result in a misspecified model. First, the KPSS test does not reject the null of cointegration. Second, there is no (visual) evidence that the linear functional form of (ref) is inappropriate. To arrive at this last conclusion, we compute $\widetilde y_{i,t} = y_{i,t}-\widehat\tau_g t^{\widehat\theta} - \widehat{\tau}_{1,i}-\widehat{\tau}_{2,i} t$ and employ the nonparametric kernel estimator from wangphillips2009 to estimate $\widetilde y_{i,t}= f(x_{i,t})+\widetilde u_{i,t}$ for each individual country.\footnote{The properties of nonparametric kernel estimators in nonlinear cointegration models have been studied by wangphillips2009, gaokanayalitjostheim2015 and wangphillips2016, among others. The latter reference is particularly relevant because it establishes that kernel estimators remain consistent and asymptotically (mixed) normal under serially correlated errors and endogeneity. None of these papers includes deterministic trends in the DGP. However, we conjecture that detrending does not affect the asymptotic properties of the kernel estimator due to the high convergence rates of the trend parameters in comparison to the slow convergence rates of the nonparametric estimator. Our bandwidth choice is $h=T^{-1/3}$.} Figure (ref) shows the nonparametric estimate in blue and the fit of Model (ref) in red. After removal of the global trend, there are some temporary departures from linearity but there is little curvature overall and certainly no visual turning point. In Table (ref), we formally test the null of linearity using the model specification test as outlined in section 3 of wangphillips2016. Based on the full sample, linearity is rejected for Austria only. A comparison with the 95% confidence intervals of the kernel estimate (Figure (ref)) suggests that this rejection is caused by the sharp decline in CO\textsubscript{2} emissions during World War II. We subsequently repeat the analysis using the $T=69$ observations after 1945. Linearity is never rejected.\footnote{The high $p$-values in Table (ref) are caused by visually small deviations from the linear trend (see Section (ref) in the Supplement for detailed graphs). Also, the relatively small sample size (for nonparametric settings) might adversely affect power. Matlab functions for nonparametric kernel regression and specification test are available at \texttt{https://github.com/HannoReuvers}. For remarks on bandwidth choice and detrending, we refer to footnote (ref).} All this align well with our earlier findings of a relevant global trend and irrelevant quadratic effects in log GDP per capita. \begin{figure}[h] \center \caption{The 95% (point-wise) confidence intervals of the non-parametric kernel estimate for the relationship between GDP and CO\textsubscript{2} emissions (blue) after removal of the country-specific and joint flexible deterministic trends. The red line indicates the linear fit implied by the estimation results of Model (ref). As the sample covers the years 1870--2014 there are several observations during World War I and World War II. The affected ranges of GDP are indicated in grey.} \end{figure} \begin{table}[t] \caption{Linearity test results. The linearity test is based on the model specification test documented in section 3 of wangphillips2016. The test is based on the integrated weighted squared deviations between the data and the linear model fit. We report the integration range, the (standardized) test statistic, and the $p$-value. Under $H_0$, the relationship between (detrended) log CO\textsubscript{2} emissions and log GDP per capita is linear.} \begin{threeparttable} \begin{tabular}{c c c c c c c c c c} \toprule & \multicolumn{3}{c}{Full Sample} & & \multicolumn{3}{c}{After World War II} \\ \cmidrule{2-4} \cmidrule{6-8} & range & $\frac{\phi}{\tau_0 \sqrt{n} h} T_n$ & $p$-value && range & $\frac{\phi}{\tau_0 \sqrt{n} h} T_n$ & $p$-value \\ \midrule Austria & [8.003,10.635] & 8.987 & 0.000 && [8.129,10.635] & 0.162& 0.871 \\ Belgium & [8.389,10.553] & 0.299 & 0.765 && [8.923,10.553] & 0.076& 0.939 \\ Finland & [7.494,10.602] & 0.105 & 0.916 && [8.693,10.602] & 0.020& 0.984 \\ Netherlands & [8.469,10.728] & 0.046 & 0.964 && [8.990,10.728] & 0.022& 0.982 \\ Switzerland & [8.708,10.993] & 0.028 & 0.977 && [9.925,10.993] & 0.009& 0.993 \\ UK & [8.641,10.510] & 0.022 & 0.982 && [9.242,10.510] & 0.021& 0.984 \\ \bottomrule \end{tabular} \begin{tablenotes} • Note: The asymptotic properties of $\frac{\phi}{\tau_0 \sqrt{n} h} T_n$ are established in wangphillips2016. That is, under suitable conditions, $\frac{\phi}{\tau_0 \sqrt{n} h} T_n\to L_W(1,0)$ as $n\to\infty$ with $L_W(1,0)$ denoting the sojourning time of a standard Brownian motion around zero during the time interval $[0,1]$. The $p$-values are computed using the cumulative distribution function of $L_W(1,0)$, see (2.11) in donggaotjostheimyin2017. \end{tablenotes} \end{threeparttable} \end{table} The preceding analysis suggests that the global flexible trend captures omitted determinants of CO\textsubscript{2} emission levels that have been decreasing over time. In their analysis, grossmankrueger1995 already included a global deterministic trend in their model because they “\emph{did not want to attribute to national income growth any improvements in local environmental quality that might actually be due to global advances in the technology for environmental preservation or to an increased global awareness of the severity of environmental problems}”. Indeed, since reliable data on green technology adaptation\footnote{nordhaus2014 discusses the link between climate change and technological changes. As another example, Figure 2 in gillinghamstock2018 reports a steady decline in the price of solar panels and a steady growth in solar panel sales. Cheaper solar energy can substitute fossil energy thereby reducing pollution.} and global awareness is scarcely available (certainly for time horizons allowing for a cointegration analysis), these variables are likely missing and thus requiring a proxy. Similar remarks are applicable to variables such as pollution control policies\footnote{A policy variable, `Repudiation of Contracts by Government', was included by panayotou1997 to proxy the quality of environmental policies and institutions.}. In reduced-form models, an EKC finding is typically explained by national income being the proxy for these omitted variables. That is, at higher levels of national income, countries have access to cleaner technologies and its citizens show greater appreciation for the environment and pollution legislation. The current analysis contradicts these income effects and points towards improvements being captured by a global trend. Our final model specification, Model (ref), is linear in log GDP per capita. Moreover, for a given year, the coefficient estimates suggest that increasing national income by 1% implies an \emph{increase} in carbon-dioxide emissions of about 1%--2.5% (depending on the country). This result seems plausible for non-carbon-neutral economies. However, CO\textsubscript{2} emissions in Austria, Belgium, Finland, the Netherlands, Switzerland, and the UK are jointly reducing at the end of the sample. What cause these global emission reductions? mazzantimusolesi2013 suggest that conglomerates of countries anticipate and respond to international climate agreements such as the Rio convention (1992) and Kyoto protocol (adopted in 1997; operational since 2005). Interestingly, the latter agreement contains emission reduction targets to be reached in 2020 and such “working-towards-a-common-reduction-deadline” does point towards a time effect.\footnote{According to the Doha amendment of the Kyoto protocol, the reduction commitments were 92% (over the period 2008--2012) and 80% (over the period 2013--2020) of 1990 emission levels for Austria, Belgium, Finland, the Netherlands and the UK. For Switzerland, the reduction target was also 92% (over the period 2008--2012) but 84.2% (over the period 2013--2020). (source: https://unfccc.int/files/kyoto_protocol/application/pdf/kp_doha_amendment_english.pdf).} Alternatively, given our sample of European countries, EU coordinated emission reduction efforts like the EU Emissions Trading System (ETS) can be a driving force behind these common emission decreases. \section{Summary and Conclusion} In this paper we have extended the Seemingly Unrelated Cointegrating Polynomial Regression (SUCPR) model of wagnergrabarczykhong2019 with a global power law deterministic trend. This multivariate specification allows us to disentangle national income and unobserved time effects. The importance of this separation is well-documented (see, e.g. volleberghmelenbergdijkgraaf2009 and mazzantimusolesi2013) but a methodological approach accounting for nonstationary regressors is currently unavailable. We fill this gap. The unknown powers of the global trend are estimated jointly with the parameters in the cointegrating relation. The limiting distribution is nonstandard due to a non-diagonal scaling matrix and second order bias terms. We therefore suggest a simulation-based approach to conduct inference. The usual subsampling KPSS-type for stationarity of the innovations of the nonlinear cointegrating relation remains valid. Our results are supported by Monte Carlo simulation. The empirical application on the Environmental Kuznets Curve shows that a flexible trend can fully capture the nonlinearity in the data thereby making higher order powers of log per capita GDP redundant. Our resulting model is linear in log per capita GDP and suggests an alternative explanation in which time effects -- e.g. technological progress, increasing environmental awareness, tightening pollution policy -- rather than economic growth cause the recent slowdown in CO\textsubscript{2} emissions. Contrary to the opening quote in the introduction, our analysis suggests that CO\textsubscript{2} emissions increase with economic growth. Carbon dioxide emissions do decrease due to time effects. \section{Acknowledgements} Earlier versions of this paper have been presented at the 2019 CFE meeting in London, the Econometrics Internal Seminar (EIS) at Erasmus University Rotterdam, and the Brownbag Seminar at Vrije Universiteit Amsterdam. We gratefully acknowledge the comments by the participants. We extend our thanks to Eric Beutner, Dick van Dijk, Stephan Smeekes, and Xiaohu Wang for their valuable feedback. All remaining errors are our own.

spacing{1.52}