Monitoring Cointegrating Polynomial Regressions: Theory and Application to the Environmental Kuznets Curves for Carbon and Sulfur Dioxide Emissions

Abstract

This paper develops residual-based monitoring procedures for cointegrating polynomial regressions (CPRs), i.e., regression models including deterministic variables and integrated processes, as well as integer powers, of integrated processes as regressors. The regressors are allowed to be endogenous, and the stationary errors are allowed to be serially correlated. We consider five variants of monitoring statistics and develop the results for three modified least squares estimators for the parameters of the CPRs. The simulations show that using the combination of self-normalization and a moving window leads to the best performance. We use the developed monitoring statistics to assess the structural stability of environmental Kuznets curves (EKCs) for both CO${}_2$ and SO${}_2$ emissions for twelve industrialized countries since the first oil price shock.

1. Introduction

This paper develops residual-based monitoring procedures for structural change in cointegrating polynomial regressions (CPRs), using the terminology of Wagner and Hong (2016). CPRs are regression models that include as explanatory variables deterministic terms, integrated processes, and integer powers of integrated processes. The regressors are allowed to be endogenous, and the stationary errors are allowed to be serially correlated. Structural change—at an unknown point in time—can occur in two facets: First, the relationship may turn into a spurious relationship.1 Second, the parameters of the relationship may change. The developed monitoring statistics extend those of Wagner and Wied (2017) in two dimensions. First, a variety of monitoring statistics is considered, including self-normalized versions and moving window detectors.2 Second, the approach is extended from cointegrating linear to cointegrating polynomial regressions.

All considered monitoring statistics are based on parameter estimation for the CPR relationship over a calibration period known to be—or at least assumed to be—free of structural change, an approach to monitoring inspired by Chu et al. (1996). With regressors that are potentially endogenous and errors that are potentially serially correlated, appropriately modified least squares estimators have to be employed to allow for the construction of nuisance parameter-free limiting distributions of the detectors obtained by scaling out a scalar long-run variance parameter. We consider the CPR-versions of three well-known estimators: Fully Modified Ordinary Least Squares (FM-OLS) considered in the CPR context in, e.g., Wagner and Hong (2016), Dynamic Ordinary Least Squares (D-OLS) considered for more general functions in Choi and Saikkonen (2010), and Integrated Modified Ordinary Least Squares (IM-OLS) considered in Vogelsang and Wagner (2014b).3 In the general CPR case, however, even the usage of the mentioned modified least squares estimators is not sufficient for nuisance parameter-free limiting distributions, and the assumption of full design, using the terminology of Vogelsang and Wagner (2014b), is required. Full design means that the limiting distributions of the modified estimators can be written as the product of a regular non-random matrix and a functional of standard Brownian motions and deterministic components. This allows to construct monitoring statistics, based on the residual limit processes, which are proportional up to a scalar long-run variance to functionals of standard Brownian motions and deterministic components. Scaling out the long-run variance, either by self-normalization or by standardization, then leads to nuisance parameter-free limiting distributions of the monitoring statistics. Note that full design, albeit not generally prevalent in the CPR case, holds in a variety of empirically relevant settings, including cointegrating linear regressions, cointegrating polynomial regressions where only one of the integrated regressors occurs as regressor also with higher order powers, and Translog-type relationships (see, e.g., Christensen et al. 1971). Environmental Kuznets curves (EKCs), which are the focus of our paper, typically include only one integrated regressor and its powers and are, therefore, of full design.

We perform a detailed simulation study of the five considered variants of the monitoring procedure in combination with the three mentioned parameter estimation methods. The performance dimensions considered include null rejection probabilities, (empirical) size-corrected power, as well as detection delays. It turns out that the combination of self-normalization and a moving window leads almost throughout to the best performance in terms of lowest over-rejections, whilst exhibiting very favorable size-corrected power properties and short delays. The choice of the estimator, FM-, D-, or IM-OLS, affects the results in particular for small samples. IM-OLS mostly leads to smaller over-rejections under the null hypothesis than FM-OLS, which in turn outperforms D-OLS. With respect to size-corrected power IM-OLS is outperformed, as expected, by both FM-OLS and D-OLS. In addition, the delays are often a bit smaller for FM-OLS than for IM-OLS. Since the differences are often relatively small, there is no clear choice between IM-OLS and FM-OLS.

We use the developed monitoring tools to assess the stability of environmental Kuznets curves (EKCs) for CO${}_2$ and SO${}_2$ emissions for twelve countries using a calibration period 1946–1973 and a monitoring period 1974–2016. The EKC hypothesis postulates an inverted U-shaped relationship between the level of economic development and pollution or emissions.4 Brock and Taylor (2005) or Kijima et al. (2010) provide survey discussions of the links between economic activity or growth and the environment.5 We, thus, use our monitoring tools to assess whether and if so when the relationship between emissions and economic activity has changed after the first oil price shock, which has led to fundamental changes in economic activity triggered not least by changing energy prices, but also by changes in environmental legislation that has been put in place in the 1970s.6 Taking into account that the polynomial functional form should probably be interpreted more as an approximation to an underlying relationship of unknown form rather than as a true relationship, of course, implies a trade-off (at least in-sample) between finding structural breaks and approximation with a higher polynomial degree.7 Against this background, it may be interesting to use monitoring tools to study whether, and at which point in time, the EKC relationship needs to be modeled with more curvature: Over the calibration period 1946–1973, for most countries (as minimum polynomial degree), a cointegrating linear relationship prevails.8 This is not too big of a surprise since, about until the mid 1970s, economic activity expanded roughly in line with emissions in many countries. Only thereafter, and triggered—as mentioned—by price, technical, and legislative changes the economic activity and pollution start to be decoupled to a certain extent. From a CPR perspective, this could mean either a structural change in the parameters of a relationship of given degree or a change to a relationship with a higher polynomial degree, e.g., from a linear to a quadratic relationship. For both CO${}_2$ and SO${}_2$ emissions, for nine of the twelve countries structural breaks are detected. The detected break points are, in some cases, quite late, which most likely reflects the delays inherent in monitoring procedures. The evidence, when considering also the full sample, is mixed with respect to structural change in the parameters but unchanged polynomial degree or structural change also with respect to the polynomial degree. The monitoring decisions lead to, as a simple empirical cross-check, good results in the following sense: For those country-pollutant combinations where no structural break is detected, using the specification and parameter values from the calibration period leads to good fit also for the full period until 2016, with obviously even better fit when re-estimating the relationship over the full sample.

The paper is organized as follows: Section 2 contains the setting, assumptions, monitoring statistics, and asymptotic results. Section 3 discusses finite sample simulation results. Section 4 presents the monitoring application to CO₂ and SO₂ emissions data. Section 5 briefly summarizes and concludes. Appendix A contains all proofs. In addition, four supplementary appendices are available: Supplementary Appendix B discusses local asymptotic power properties. Supplementary Appendix C contains additional finite sample simulation results, and Supplementary Appendix D presents additional empirical results. Supplementary Appendix E includes tables with critical values for the detectors for a broad variety of specifications relevant for EKC-type analysis. MATLAB code—including the critical values for the mentioned variety of specifications—for the monitoring statistics developed in this paper is available upon request.9

We use the following notation: Definitional equality is signified by := and ⇒ denotes weak convergence. $\lfloor x\rfloor$ denotes the integer part of $x\in \mathbb{R}$ and diag$(·)$ denotes a diagonal matrix. For a vector $x\in {\mathbb{R}}^n$ we denote by ${\parallel x\parallel }^2:={\sum }_{i=1}^n{x}_i^2$ and for a matrix $M\in {\mathbb{R}}^{m\times n}$ we denote by $\parallel M\parallel :={sup}_x\frac{\parallel Mx\parallel }{\parallel x\parallel }$. $\mathbb{E}(·)$ denotes the expected value and L is the backward-shift operator, i.e., $L{\left\{x_t\right\}}_{t\in \mathbb{Z}}:={\left\{x_{t-1}\right\}}_{t\in \mathbb{Z}}$. The first-difference operator is denoted with $\Delta :=1-L$. We denote the k- dimensional identity matrix with $I_k$. A Brownian motion with covariance matrix specified in the context is denoted by $B\left(r\right)$, and $W\left(r\right)$ denotes a standard Brownian motion.

2. Theory

2.1. Model, Assumptions, and Parameter Estimation

We consider monitoring—using the terminology of Wagner and Hong (2016)—a cointegrating polynomial regression (CPR), i.e., a regression of the form:

$$\begin{array}{cc}\hfill y_t& =\left\{\begin{array}{cc}{D}_t^{\prime }{\theta }_D+X_t^{\prime }{\theta }_X+u_t,\hfill & t=1,\dots ,\lfloor rT\rfloor ,\hfill \\ D_t^{\prime }{\theta }_{D,1}+X_t^{\prime }{\theta }_{X,1}+u_t,\hfill & t=\lfloor rT\rfloor +1,\dots ,T,\hfill \end{array}\right.\hfill \end{array}$$

(1)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill x_t& =x_{t-1}+v_t,t=1,\dots ,T,\hfill \end{array}\end{array}$$

(2)

with $x_t:={[x_{1t},\dots ,x_{kt}]}^{\prime }\in {\mathbb{R}}^k$ and $X_t:={[x_t^{\prime },x_{kt}^2,\dots ,x_{kt}^{p_k}]}^{\prime }\in {\mathbb{R}}^p$ with $p=k-1+p_k$, the deterministic trend function $D_t\in {\mathbb{R}}^q$, the parameter vectors ${\theta }_D,{\theta }_{D,1}\in {\mathbb{R}}^q$ and ${\theta }_X,{\theta }_{X,1}\in {\mathbb{R}}^p$. Furthermore, we define the combined parameter vectors $\theta :={[{\theta }_D^{\prime },{\theta }_X^{\prime }]}^{\prime },{\theta }_1:={[{\theta }_{D,1}^{\prime },{\theta }_{X,1}^{\prime }]}^{\prime }\in {\mathbb{R}}^{q+p}$. Note that, as, e.g., discussed in Wagner (2012) and Stypka and Wagner (2019), polynomial transformations of integrated processes—in our case, $\left\{y_t\right\}$—are not integrated processes (of any order). This is due to the basic fact that summation (integration) and polynomial transformation do not commute.

Under the null hypothesis, no structural change occurs, that is ${\theta }_1=\theta$ and ${\left\{u_t\right\}}_{t\in \mathbb{Z}}$ is an I(0) process, with detailed assumptions specified below, throughout. Under the alternative hypothesis, either some parameters change or the relation turns spurious, i.e., ${\left\{u_t\right\}}_{t\in \mathbb{Z}}$ turns into an I(1) process for every $\theta \in {\mathbb{R}}^{q+p}$, or both at a sample fraction $\lfloor rT\rfloor$ that has to be—as discussed in the introduction—larger than $\lfloor mT\rfloor$ for some $0<m<1$.10 In formal terms:

$$\begin{array}{cc}\hfill & H_0:\left\{\begin{array}{c}{\theta }_1=\theta \forall r\in [m,1)\mathrm{and}\hfill \\ u_t\mathrm{is}I\left(0\right)\mathrm{for}t=1,\dots ,T\hfill \end{array}\right.\hfill \\ \\ \hfill & H_1:\left\{\begin{array}{c}\exists r\in [m,1):{\theta }_1\ne \theta \mathrm{or}\hfill \\ \exists r\in [m,1):u_t,t=1,\dots ,\lfloor rT\rfloor \mathrm{is}I\left(0\right)\mathrm{and}\hfill \\ u_t,t=\lfloor rT\rfloor +1,\dots ,T,\mathrm{is}I\left(1\right).\hfill \end{array}\right.\hfill \end{array}$$

Remark 1.

Note that the regression model given in (1) is a special case of the CPR model considered in Wagner and Hong (2016) since only one of the integrated regressors, w.l.o.g. $x_{kt}$, is allowed to enter the regression model with powers larger than one. Whilst this is, obviously, restrictive compared to the case where higher order powers of all elements of $x_t$ can be present as regressors, this special case covers EKCs and MKCs and similar applications. The mathematical reason for considering this special case is that it allows for—potentially up to a scalar nuisance parameter that can be consistently estimated and scaled out—nuisance parameter-free limiting distributions of the considered detectors that can be simulated. In the terminology of Vogelsang and Wagner (2014b), a cointegrating polynomial regression needs to exhibit full design to allow for asymptotic standard inference by scaling out a scalar long-run variance. This is the case for EKC-type relationships with only one integrated regressor present with powers larger than one but also for some other economically relevant more general cases of CPRs, e.g., Translog functions (see, e.g., Christensen et al. 1971). Given our focus on EKCs, we abstain from formulating the results here in the most general form; the required extensions are straightforward. For even more general specifications that do not exhibit full design, a sub-sampling approach may be considered relying upon similar arguments as discussed in Wagner and Hong (2016, Proposition 6). The performance of sub-sampling-based procedures, however, suffers particularly from short sample periods, as also illustrated by the simulations reported in Wagner and Hong (2016). Consequently, a sub-sampling-based approach cannot be expected to perform well in a monitoring context.

The results developed below rest upon the following assumptions:

Assumption 1.

There exists a sequence of $q\times q$ scaling matrices $G_D\left(T\right)$ and a q-dimensional vector of functions D(z), with ${\int }_0^{s}D\left(z\right)D{\left(z\right)}^{\prime }dz<\infty$ for $0\le s\le 1$, such that, for $0\le s\le 1$, it holds that:

$$\underset{T\to \infty }{lim}\sqrt{T}{G}_D^{-1}\left(T\right)D_{\lfloor sT\rfloor }=D\left(s\right).$$

(3)

If, e.g., $D_t:={[1,t,t^2,\dots ,t^{q-1}]}^{\prime }$, then $G_D\left(T\right):=\mathrm{diag}(T^{1/2},T^{3/2},T^{5/2},\dots ,T^{q-1/2})$ and $D\left(z\right)={[1,z,z^2,\dots ,z^{q-1}]}^{\prime }$. In relation to the integrated regressors and the powers, we need a scaling matrix $G_X\left(T\right):=\mathrm{diag}(T\times I_k,T^{3/2},\dots ,T^{\frac{p_k+1}{2}})$ later.

Assumption 2.

The process ${\left\{{\eta }_t\right\}}_{t\in \mathbb{Z}}:={\left\{{[u_t,v_t^{\prime }]}^{\prime }\right\}}_{t\in \mathbb{Z}}$ is generated under the null hypothesis as:

$${\eta }_t=C\left(L\right){\xi }_t=\sum _{j=0}^{\infty }{C}_j{\xi }_{t-j},t\in \mathbb{Z},$$

(4)

with $C_j\in {\mathbb{R}}^{(k+1)\times (k+1)},j=0,1,\dots ,$ and the conditions:

$$det\left(C\left(1\right)\right)\ne 0,\sum _{j=0}^{\infty }j\parallel C_j\parallel <\infty and\parallel C_0\parallel <\infty .$$

(5)

Furthermore, we assume that the process ${\left\{{\xi }_t\right\}}_{t\in \mathbb{Z}}$ is a strictly stationary and ergodic martingale difference sequence with natural filtration ${\mathcal{F}}_t:=\sigma \left({\left\{{\xi }_s\right\}}_{-\infty }^t\right)$, $\mathbb{E}\left({\xi }_t{\xi }_t^{\prime }\right|{\mathcal{F}}_{t-1})=:{\mathsf{\Sigma }}_{\xi \xi }>0$ with, in addition ${sup}_{t\ge 1}\mathbb{E}(\parallel {\xi }_t{\parallel }^a\left|{\mathcal{F}}_{t-1}\right)<\infty$ a.s., for some $a>4$.

The assumptions on the deterministic components, the regressors, and error terms are similar to the assumptions used in Wagner and Hong (2016) and, more implicitly, in Wagner and Wied (2017). In particular, Assumption 2 is sufficient for a functional central limit theorem to hold for ${\left\{{\eta }_t\right\}}_{t\in \mathbb{Z}}$:

$$\begin{array}{c}\hfill \frac{1}{\sqrt{T}}\sum _{t=1}^{\lfloor sT\rfloor }{\eta }_t\Rightarrow B\left(s\right)=\left[\begin{array}{c}{B}_u\left(s\right)\\ B_v\left(s\right)\end{array}\right]={\Omega }^{1/2}W\left(s\right),0\le s\le 1,\end{array}$$

(6)

with the positive definite long-run covariance matrix $\Omega :={\sum }_{j=-\infty }^{\infty }\mathbb{E}\left({\eta }_{t-j}{\eta }_t^{\prime }\right)$ and $W\left(s\right):={[W_{u·v}\left(s\right),W_v{\left(s\right)}^{\prime }]}^{\prime }$ a $(k+1)$-dimensional vector of standard Brownian motions. We also define the half long-run covariance matrix $\Delta :={\sum }_{j=0}^{\infty }\mathbb{E}\left({\eta }_{t-j}{\eta }_t^{\prime }\right)$. The matrices $\Omega$ and $\Delta$ are partitioned according to the partitioning of $B\left(s\right)$, i.e.,

$$\begin{array}{c}\hfill \Omega =\left[\begin{array}{cc}{\Omega }_{uu}& {\Omega }_{uv}\\ {\Omega }_{vu}& {\Omega }_{vv}\end{array}\right],\Delta =\left[\begin{array}{cc}{\Delta }_{uu}& {\Delta }_{uv}\\ {\Delta }_{vu}& {\Delta }_{vv}\end{array}\right].\end{array}$$

(7)

Using, e.g., the Cholesky decomposition of $\Omega$ yields:

$$\begin{array}{c}\hfill {\Omega }^{1/2}=\left[\begin{array}{cc}{\omega }_{u·v}& {\lambda }_{uv}\\ 0& {\Omega }_{vv}^{1/2}\end{array}\right],\end{array}$$

(8)

where ${\omega }_{u·v}^2:={\Omega }_{uu}-{\Omega }_{uv}{\Omega }_{vv}^{-1}{\Omega }_{vu}$ and ${\lambda }_{uv}:={\Omega }_{uv}{\left({\Omega }_{vv}^{1/2}\right)}^{-1}$. The conditional long-run variance ${\omega }_{u·v}^2$ is a key quantity that needs to be estimated for all but the two self-normalized detectors.

Where needed, consistent long-run covariance estimation is performed non-parametrically, requiring the choice of both a kernel function and a bandwidth parameter. The inputs in the non-parametric estimation are the OLS residuals from estimating (1) over the calibration period and the first difference of $x_t$ over the same period. For consistent long-run covariance estimation, it suffices to assume (following, e.g., Jansson 2002):

Assumption 3.

The kernel function $k(·)$ satisfies:

(i) 

$k\left(0\right)=1$, $k(·)$ is continuous at 0 and $\overline{k}\left(0\right):={sup}_{x\ge 0}\left|k\left(x\right)\right|<\infty$,

(ii) 

${\int }_0^{\infty }\overline{k}\left(x\right)dx<\infty$, where $\overline{k}\left(x\right)={sup}_{y\ge x}\left|k\left(y\right)\right|$.

Assumption 4.

The bandwidth parameter $M_T\subseteq (0,\infty )$ satisfies ${lim}_{T\to \infty }(M_T^{-1}+T^{-1/2}{M}_T)=0$.

All our monitoring statistics discussed in the following subsection are based on consistent parameter estimators that are required to lead to limiting distributions that are nuisance parameter free up to a scalar parameter, the conditional long-run variance ${\omega }_{u·v}^2$, that can (asymptotically) be scaled out, either by scaling by a consistent estimator, which we refer to later as standardized, or by self-normalization.

As indicated, estimation takes place on the calibration sample $t=1,\dots ,\lfloor mT\rfloor$ for some $0<m<1$ that is known to be generated under the null hypothesis. This approach to monitoring, based on parameter estimation on a calibration sample known to be—or at least assumed to be—free of structural change, has been popularized in the econometrics community by the seminal work of Chu et al. (1996).

The cointegration literature provides a variety of modified ordinary least squares estimators of $\theta$ with the required asymptotic properties; see, e.g., Wagner (2018), for a survey. All these estimators commence from the fact that the OLS estimator of $\theta$ is consistent with—in case of regressor endogeneity and error serial correlation—a limiting distribution that is contaminated by second order bias terms. These second order bias terms are removed, one way or another, by the various modifications of OLS. In this paper, we consider three modified OLS estimators: Fully Modified OLS (FM-OLS), Dynamic OLS (D-OLS), and Integrated Modified OLS (IM-OLS). These three estimators were originally developed for cointegrating linear regressions: FM-OLS in Phillips and Hansen (1990); D-OLS in Saikkonen (1991), Phillips and Loretan (1991), or Stock and Watson (1993); and IM-OLS in Vogelsang and Wagner (2014a). The extensions to the CPR setting are discussed for FM-OLS in Wagner and Hong (2016), for D-OLS in Choi and Saikkonen (2010), and for IM-OLS in Vogelsang and Wagner (2014b).11

Our brief discussion of the three estimators first necessitates the definition of a few more quantities, i.e., $Z_t:={[D_t^{\prime },X_t^{\prime }]}^{\prime }$ and $y_{t,m}^{+}:=y_t-\Delta x_t^{\prime }{\widehat{\Omega }}_{vv,m}^{-1}{\widehat{\Omega }}_{vu,m}$, with the second subscript m indicating that estimation of the long-run covariances is—as mentioned—also based on the calibration sample $t=1,\dots ,\lfloor mT\rfloor$. Furthermore, define:

$$\begin{array}{cc}\hfill A_m^{*}:=& \left[\begin{array}{c}{0}_{q\times 1}\\ \lfloor mT\rfloor {\widehat{\Delta }}_{vu,m}^{+}\\ M_m^{*}\end{array}\right],M_m^{*}:={\widehat{\Delta }}_{v_{k}u,m}^{+}\left[\begin{array}{c}2\sum _{t=1}^{\lfloor mT\rfloor }{x}_{kt}\\ ⋮\\ p_k\sum _{t=1}^{\lfloor mT\rfloor }{x}_{kt}^{p_k-1}\end{array}\right],\hfill \end{array}$$

(9)

with ${\widehat{\Delta }}_{vu,m}^{+}:={\widehat{\Delta }}_{vu,m}{\widehat{\Omega }}_{vv,m}^{-1}{\widehat{\Delta }}_{vv,m}$ and ${\widehat{\Delta }}_{v_{k}u,m}^{+}{\widehat{\Delta }}_{v_{k}u,m}{\widehat{\Omega }}_{vv,m}^{-1}{\widehat{\Delta }}_{v{v}_k,m}$.

Long-run covariance estimation uses the OLS residuals of (1) from estimation over the calibration period $t=1,\dots ,\lfloor mT\rfloor$ in conjunction with the first differences of the integrated regressors, i.e., ${\widehat{\eta }}_{t,m}:={[{\widehat{u}}_{t,m},v_t^{\prime }]}^{\prime }$, with ${\widehat{u}}_{t,m}$ denoting the OLS residuals here:

$$\begin{array}{cc}\hfill {\widehat{\Delta }}_m:=& \frac{1}{\lfloor mT\rfloor }\sum _{h=0}^{\lfloor mT\rfloor -1}k\left(\frac{h}{M_T}\right)\sum _{t=1}^{\lfloor mT\rfloor -h}{\widehat{\eta }}_{t,m}{\widehat{\eta }}_{t+h,m}^{\prime },\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\widehat{\mathsf{\Sigma }}}_m:=& \frac{1}{\lfloor mT\rfloor }\sum _{t=1}^{\lfloor mT\rfloor }{\widehat{\eta }}_{t,m}{\widehat{\eta }}_{t,m}^{\prime },\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\widehat{\Omega }}_m:=& {\widehat{\Delta }}_m+{\widehat{\Delta }}_m^{\prime }-{\widehat{\mathsf{\Sigma }}}_m.\hfill \end{array}$$

(12)

In both the finite sample simulations, as well as the application, we use for long-run covariance estimation, in line with Assumptions 3 and 4, the Bartlett kernel with bandwidth chosen according to Newey and West (1994).

With all required quantities defined, the FM-OLS estimator computed over the calibration sample is given by:

$${\widehat{\theta }}_m^F:={\left(\sum _{t=1}^{\lfloor mT\rfloor }{Z}_t{Z}_t^{\prime }\right)}^{-1}\left(\sum _{t=1}^{\lfloor mT\rfloor }{Z}_t{y}_{t,m}^{+}-A_m^{*}\right).$$

(13)

While FM-OLS is based on a two-part nonparametric transformation to remove endogeneity and serial correlation related bias terms from the limiting distribution of the OLS estimator, D-OLS is based on a more “classical projection and orthogonalization” argument by performing OLS estimation in an augmented version of the CPR regression (1), with leads and lags of the first differences of $x_t$ added as regressors to “clean the limiting distribution”. The D-OLS regression—estimated by OLS over the calibration sample—is given by:

$$y_t=Z_t^{\prime }\theta +\sum _{j=-d_1}^{d_2}\Delta x_{t-j}^{\prime }{\mathsf{\Theta }}_j+u_t,$$

(14)

with the number of leads ($d_1$) and lags ($d_2$) chosen to ensure consistent parameter estimation of $\theta$ with a limiting distribution that is—up to a scalar—nuisance parameter-free. In general, this requires that $d_1,d_2\to \infty$ at suitable rates. More specifically, in our finite sample simulations and the application, we choose leads and lags using the AIC-type criterion of Choi and Kurozumi (2012). The resultant OLS estimators of $\theta$ and ${\widehat{\mathsf{\Theta }}}_j$ from (14), estimated over the calibration sample, are referred to as ${\widehat{\theta }}_m^{\mathrm{D}}$ and ${\widehat{\mathsf{\Theta }}}_{j,m}^{\mathrm{D}}$, respectively.

The third estimation principle addresses endogeneity correction by partial summation. Define for a sequence $z_t$, $t=1,\dots ,T$, the partial summed variable by $S_t^z:={\sum }_{j=1}^t{z}_j,t=1,\dots ,T$. Then, the IM-OLS regression—estimated by OLS over the calibration sample—is given by:

$$S_t^y=S_t^{Z\prime }\theta +x_t^{\prime }\phi +S_t^u.$$

(15)

The OLS estimators of $\theta$ and $\phi$ from (15) estimated over the calibration sample are referred to as ${\widehat{\theta }}_m^{\mathrm{I}}$ and ${\widehat{\phi }}_m^{\mathrm{I}}$, respectively. Note that endogeneity correction in the IM-OLS estimator does not require any lead-lag or kernel-bandwidth choices, as it suffices to simply add the original integrated regressor vector $x_t$ to the partial summed regression.

The key input for the monitoring statistics discussed in the following subsection are the residual (processes) obtained with these three estimators. In particular, the asymptotic null behavior of the residual (partial sum) processes is the key ingredient to derive asymptotic properties of any of our variance-ratio type monitoring statistics. This result is formalized in the following lemma, for which formulation requires to define some additional (asymptotic) quantities first, i.e., ${\mathbf{W}}_v\left(s\right):={[W_{v_1}\left(s\right),W_{v_2}\left(s\right),\dots ,W_{v_k}\left(s\right),W_{v_k}^2\left(s\right),\dots ,W_{v_k}^{p_k}\left(s\right)]}^{\prime }$, $J\left(s\right):={[D{\left(s\right)}^{\prime },{\mathbf{W}}_v{\left(s\right)}^{\prime }]}^{\prime }$, $f\left(s\right):={[{\int }_0^{s}D{\left(z\right)}^{\prime }dz,{\int }_0^s{\mathbf{W}}_v{\left(z\right)}^{\prime }dz,W_v{\left(s\right)}^{\prime }]}^{\prime }$ and $F\left(s\right):={\int }_0^{s}f\left(z\right)dz$.

Lemma 1.

Let the data be generated according to (1) and (2) with Assumptions 1 and 2 in place. Furthermore, let long-run covariance estimation be performed under Assumptions 3 and 4 and let lead-lag choices be made as discussed in Choi and Kurozumi (2012).

Based on estimation over the calibration period $t=1,\dots ,\lfloor mT\rfloor$—with the estimators ${\widehat{\theta }}_m^{\mathrm{F}}$, ${\widehat{\theta }}_m^{\mathrm{D}}$ and ${\widehat{\theta }}_m^I$ as discussed before—denote the FM-OLS residuals by:

$${\widehat{u}}_{t,m}^{\mathrm{F}}:=y_{t,m}^{+}-Z_t^{\prime }{\widehat{\theta }}_m^{\mathrm{F}},$$

(16)

the D-OLS residuals by:

$${\widehat{u}}_{t,m}^{\mathrm{D}}:=y_t-Z_t^{\prime }{\widehat{\theta }}_m^{\mathrm{D}}-\sum _{j=-d_1}^{d_2}\Delta x_{t-j}^{\prime }{\widehat{\mathsf{\Theta }}}_{j,m}^{\mathrm{D}},$$

(17)

and the IM-OLS residuals by:

$${\widehat{S}}_{t,m}^{u,\mathrm{I}}:=S_t^y-S_t^{Z\prime }{\widehat{\theta }}_m^{\mathrm{I}}-x_t^{\prime }{\widehat{\phi }}_m^{\mathrm{I}}.$$

(18)

For $T\to \infty$, it holds under the null hypothesis for $m\le s\le 1$ that:12

$$\begin{array}{cc}\hfill \frac{1}{\sqrt{T}}{\widehat{S}}_{\lfloor sT\rfloor ,m}^{u,\mathrm{F}}:=\frac{1}{\sqrt{T}}\sum _{t=2}^{\lfloor sT\rfloor }{\widehat{u}}_{t,m}^{\mathrm{F}}& \Rightarrow {\omega }_{u·v}\left(W_{u·v}\left(s\right)-{\int }_0^{s}J{\left(z\right)}^{\prime }dz{\left({\int }_0^{m}J\left(z\right)J{\left(z\right)}^{\prime }dz\right)}^{-1}{\int }_0^{m}J\left(z\right)d{W}_{u·v}\left(z\right)\right)\hfill \\ \hfill & =:{\omega }_{u·v}{\tilde{W}}_{u·v,m}\left(s\right),\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \frac{1}{\sqrt{T}}{\widehat{S}}_{\lfloor sT\rfloor ,m}^{u,\mathrm{D}}:=& \frac{1}{\sqrt{T}}\sum _{t=d_2+2}^{min\{\lfloor sT\rfloor ,T-d_1\}}{\widehat{u}}_{t,m}^{\mathrm{D}}\Rightarrow {\omega }_{u·v}{\tilde{W}}_{u·v,m}\left(s\right),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \frac{1}{\sqrt{T}}{\widehat{S}}_{\lfloor sT\rfloor ,m}^{u,\mathrm{I}}:=\frac{1}{\sqrt{T}}\sum _{t=2}^{\lfloor sT\rfloor }\Delta {\widehat{S}}_{t,m}^{u,\mathrm{I}}& \Rightarrow {\omega }_{u·v}\left(W_{u·v}\left(s\right)-f{\left(s\right)}^{\prime }{\left({\int }_0^{m}f\left(z\right)f{\left(z\right)}^{\prime }dz\right)}^{-1}{\int }_0^m[F\left(m\right)-F\left(z\right)]d{W}_{u·v}\left(z\right)\right)\hfill \\ \hfill & =:{\omega }_{u·v}{\tilde{P}}_{u·v,m}\left(s\right).\hfill \end{array}$$

(21)

The lemma shows that indeed all three partial sum processes of the residuals converge to processes that are (i) functionals of standard Brownian motions, $W_v\left(r\right)$ and $W_{u·v}\left(r\right)$, and (ii) proportional to ${\omega }_{u·v}$, a scalar nuisance parameter that can be consistently estimated and hence scaled out from the limit processes or that can be eliminated by self-normalization. The limiting null distributions of test statistics based on the (normalized) limit processes consequently can be obtained by simulating the corresponding functionals of standard Brownian motions. Note that the FM-OLS and D-OLS residual partial sum processes converge to the same limiting process, which is a consequence of these two estimators having identical limiting distributions.

2.2. The Monitoring Statistics

Similarly to Wagner and Wied (2017), the starting point of our monitoring statistics is to combine the approach of Chu et al. (1996) with variance-ratio statistics that diverge under the alternative. More specifically, the underlying variance-ratio (full sample) statistic motivating the construction of our monitoring statistics is the Kwiatkowski et al. (1992) stationarity test, respectively, the related Shin (1994) cointegration test.13 Using our notation, the (full sample) Shin-statistic is given by:

$$\begin{array}{cc}\hfill T_{Shin}:=& \frac{1}{{\widehat{\omega }}_{u·v}^2}\left(\frac{1}{T}\sum _{t=1}^T{\left(\frac{1}{\sqrt{T}}\sum _{i=1}^t{\widehat{u}}_i\right)}^2\right)\hfill \\ \hfill =& \frac{1}{{\widehat{\omega }}_{u·v}^2}\left(\frac{1}{T}\sum _{t=1}^T{\left(\frac{1}{\sqrt{T}}{\widehat{S}}_t^u\right)}^2\right),\hfill \end{array}$$

(22)

with ${\widehat{u}}_t$ denoting the residuals from (full sample) estimation with, e.g., FM-OLS or D-OLS. In case that IM-OLS is used for estimation, the resultant residuals are already partial summed quantities, i.e., one immediately obtains (by construction) quantities to insert into the expression in the second line of (22). The test statistic given above converges under the null hypothesis to a functional of standard Brownian motions, which is as expected when considering (22), where convergence to the squared integral of a standard Brownian motion follows immediately from our assumptions if instead of ${\widehat{u}}_t$ the errors $u_t$ were used (and scaling would take place by a consistent estimator of ${\omega }_u^2$). Using ${\widehat{u}}_t$ instead of $u_t$ leads to a similar result, but with a different (specification dependent) functional of standard Brownian motions after scaling out ${\omega }_{u·v}^2$ rather than ${\omega }_u^2$. To be precise, when using FM-OLS or D-OLS for parameter estimation, the limiting null distribution will be a function of ${\tilde{W}}_{u·v,m}\left(r\right)$. When using IM-OLS for parameter estimation, the limiting null distribution will be a function of ${\tilde{P}}_{u·v,m}\left(r\right)$; see Proposition 1 below.

The above test statistic (22) can be easily seen to diverge under the alternative hypothesis of a structural break occurring after the calibration period. Consider, e.g., the FM-OLS residuals (with the argument entirely analogous for all three considered estimators) using our already established notation:

$$\begin{array}{cc}\hfill {\widehat{u}}_{t,m}^{\mathrm{F}}:=& y_{t,m}^{+}-Z_t^{\prime }{\widehat{\theta }}_m^{\mathrm{F}}\hfill \\ \hfill =& u_t-v_t^{\prime }{\widehat{\Omega }}_{vv,m}^{-1}{\widehat{\Omega }}_{vu,m}-D_t^{\prime }({\widehat{\theta }}_{D,m}^{\mathrm{F}}-{\theta }_D)-X_t^{\prime }({\widehat{\theta }}_{X,m}^{\mathrm{F}}-{\theta }_X),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \frac{1}{\sqrt{T}}\sum _{t=2}^{\lfloor sT\rfloor }{\widehat{u}}_{t,m}^{\mathrm{F}}=& \frac{1}{\sqrt{T}}\sum _{t=2}^{\lfloor sT\rfloor }{u}_t-\frac{1}{\sqrt{T}}\sum _{t=2}^{\lfloor sT\rfloor }{v}_t^{\prime }{\widehat{\Omega }}_{vv,m}^{-1}{\widehat{\Omega }}_{vu,m}-\frac{1}{\sqrt{T}}\sum _{t=2}^{\lfloor sT\rfloor }{D}_t^{\prime }({\widehat{\theta }}_{D,m}^{\mathrm{F}}-{\theta }_D)\hfill \\ \hfill & -\frac{1}{\sqrt{T}}\sum _{t=2}^{\lfloor sT\rfloor }{X}_t^{\prime }({\widehat{\theta }}_{X,m}^{\mathrm{F}}-{\theta }_X).\hfill \end{array}$$

(24)

Now, suppose that at some time point $\lfloor rT\rfloor >\lfloor mT\rfloor$ a structural change occurs. If, e.g., ${\left\{u_t\right\}}_{t\in \mathbb{Z}}$ turns from being I(0) to I(1), then the first term in (24) diverges for $s>r$. Similarly, the third or fourth term (or both) diverge in case of change in the parameter vector, i.e., when ${\theta }_1\ne \theta$, as, of course, ${\widehat{\theta }}_m^{\mathrm{F}}\to \theta$, because of parameter estimation on the calibration sample.

We consider, for each of the three considered estimators—neglecting for notational brevity the dependence of the residuals and, thus, the test statistics on the estimation method—, five monitoring statistics:

$$\begin{array}{cc}\hfill {\widehat{H}}^m\left(s\right):=& \frac{1}{{\widehat{\omega }}_{u·v,m}^2}\left(\frac{1}{T}\sum _{i=\lfloor mT\rfloor +1}^{\lfloor sT\rfloor }{\left(\frac{1}{\sqrt{T}}{\widehat{S}}_{i,m}^u\right)}^2\right),\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill {\widehat{H}}_d^m\left(s\right):=& \frac{1}{{\widehat{\omega }}_{u·v,m}^2}\left(\frac{1}{T}\sum _{i=\lfloor mT\rfloor +1}^{\lfloor sT\rfloor }{\left(\frac{1}{\sqrt{T}}{\widehat{S}}_{i,m}^u\right)}^2-\frac{1}{T}\sum _{i=1}^{\lfloor mT\rfloor }{\left(\frac{1}{\sqrt{T}}{\widehat{S}}_{i,m}^u\right)}^2\right),\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\widehat{H}}_{sn}^m\left(s\right):=& \frac{{\sum }_{i=\lfloor mT\rfloor +1}^{\lfloor sT\rfloor }{\left({\widehat{S}}_{i,m}^u\right)}^2}{{\sum }_{i=1}^{\lfloor mT\rfloor }{\left({\widehat{S}}_{i,m}^u\right)}^2},\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\widehat{H}}_{mov}^{m,n}\left(s\right):=& \frac{1}{{\widehat{\omega }}_{u·v,m}^2}\left(\frac{1}{T}\sum _{i=max\{1,\lfloor sT\rfloor -\lfloor nT\rfloor +1\}}^{\lfloor sT\rfloor }{\left(\frac{1}{\sqrt{T}}{\widehat{S}}_{i,m}^u\right)}^2\right),\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\widehat{H}}_{mov,sn}^{m,n}\left(s\right):=& \frac{{\sum }_{i=max\{1,\lfloor sT\rfloor -\lfloor nT\rfloor +1\}}^{\lfloor sT\rfloor }{\left({\widehat{S}}_{i,m}^u\right)}^2}{{\sum }_{i=1}^{\lfloor mT\rfloor }{\left({\widehat{S}}_{i,m}^u\right)}^2}.\hfill \end{array}$$

(29)

The monitoring statistic ${\widehat{H}}^m\left(s\right)$ given in (25) is of the same form as the monitoring statistic used in Wagner and Wied (2017) considered here in the CPR context. The monitoring statistic ${\widehat{H}}_d^m\left(s\right)$ given in (26)—with a term calculated only over the calibration sample subtracted—is of a similar form as used in Chu et al. (1996). The third variant ${\widehat{H}}_{sn}^m\left(s\right)$ given in (27) is a self-normalized statistic, for which, under the null hypothesis, both the numerator and denominator converge (appropriately scaled) to functionals of standard Brownian motions proportional to ${\omega }_{u·v}$, which is hence scaled out in the ratio. Long-run covariance estimation is known to be a notoriously problematic aspect in unit root and cointegration analysis; therefore, test statistics that do not require this step may exhibit better performance. The fourth considered variant is a moving window statistic ${\widehat{H}}_{mov}^{m,n}\left(s\right)$ given in (28) with n denoting the moving window (sample fraction or) length. The key difference between the moving window detector and the expanding window detectors is that ${\widehat{H}}_{mov}^{m,n}\left(s\right)$ is based on a constant number of residual partial sums for all values of s. This construction increases, under the alternative hypothesis, the impact of post-break residuals on the test statistic, which is ex ante expected to lead to faster detection of structural breaks. The performance of the fourth variant will depend on the length of the moving window, to be chosen in applications. Finally, the fifth monitoring statistic $H_{mov,sn}^{m,n}\left(s\right)$ given in (29) combines self-normalization and moving window estimation, with the performance as for the fourth variant expected to depend upon the moving window length.14 The following proposition summarizes the asymptotic behavior of the monitoring statistics under the null hypothesis.

Proposition 1.

Let the data be generated according to (1) and (2) with Assumptions 1 and 2 in place. Furthermore, let long-run covariance estimation be performed under Assumptions 3 and 4 and let lead-lag choices be made as discussed in Choi and Kurozumi (2012). In case parameter estimation is performed with FM-OLS or D-OLS, the limiting process ${\tilde{Q}}_{u·v,m}\left(s\right)$ below equals ${\tilde{W}}_{u·v,m}\left(s\right)$, and, in case parameter estimation is performed by IM-OLS, it equals ${\tilde{P}}_{u·v,m}\left(s\right)$. The defined monitoring statistics converge under the null hypothesis for $T\to \infty$; in particular, it holds that:

$$\begin{array}{cc}\hfill {\widehat{H}}^m\left(s\right)& \Rightarrow {\int }_m^s{\tilde{Q}}_{u·v,m}^2\left(z\right)dz=:{\mathcal{H}}^m({\tilde{Q}}_{u·v,m},s),\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\widehat{H}}_d^m\left(s\right)& \Rightarrow {\int }_m^s{\tilde{Q}}_{u·v,m}^2\left(z\right)dz-{\int }_0^m{\tilde{Q}}_{u·v,m}^2\left(z\right)dz=:{\mathcal{H}}_d^m({\tilde{Q}}_{u·v,m},s),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\widehat{H}}_{sn}^m\left(s\right)& \Rightarrow \frac{{\int }_m^s{\tilde{Q}}_{u·v,m}^2\left(z\right)dz}{{\int }_0^m{\tilde{Q}}_{u·v,m}^2\left(z\right)dz}=:{\mathcal{H}}_{sn}^m({\tilde{Q}}_{u·v,m},s),\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\widehat{H}}_{mov}^{m,n}\left(s\right)& \Rightarrow {\int }_{max\{0,s-n\}}^s{\tilde{Q}}_{u·v,m}^2\left(z\right)dz=:{\mathcal{H}}_{mov}^{m,n}({\tilde{Q}}_{u·v,m},s),\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\widehat{H}}_{mov,sn}^{m,n}\left(s\right)& \Rightarrow \frac{{\int }_{max\{0,s-n\}}^s{\tilde{Q}}_{u·v,m}^2\left(z\right)dz}{{\int }_0^m{\tilde{Q}}_{u·v,m}^2\left(z\right)dz}=:{\mathcal{H}}_{mov,sn}^{m,n}({\tilde{Q}}_{u·v,m},s).\hfill \end{array}$$

(34)

It is widely-used practice in monitoring to base the decision not on monitoring statistics as just defined, but on monitoring statistics divided by a weighting function, say $g\left(s\right)$. For chosen weighting function $g\left(s\right)$—with $0<g\left(s\right)<\infty$—, the null hypothesis is rejected, if the weighted monitoring statistic $\left|\frac{\widehat{H}\left(s\right)}{g\left(s\right)}\right|$ is larger than a critical value c for the first time. We denote this point in time as detection time ${\tau }_m$, i.e.,

$${\tau }_m:=\underset{s:\lfloor mT\rfloor +1\le \lfloor sT\rfloor \le T}{min}\left\{\left|\frac{\widehat{H}\left(s\right)}{g\left(s\right)}\right|>c\right\},$$

(35)

with $\widehat{H}\left(s\right)$ short-hand notation for any of the considered detectors.15 In case no structural change is detected, i.e., $\left|\frac{\widehat{H}\left(s\right)}{g\left(s\right)}\right|\le c$ for all $m\le s\le 1$, we set ${\tau }_m=\infty$. A finite value of ${\tau }_m$ not only indicates a structural break but also contains information about the location of the potential break point.

Weighting function and critical value have to be chosen so that, under the null hypothesis, it holds that:

$$\begin{array}{ccc}\hfill \underset{T\to \infty }{lim}\mathbb{P}({\tau }_m<\infty )& =& \underset{T\to \infty }{lim}\mathbb{P}\left(\underset{s:\left[mT\right]+1\le \left[sT\right]\le T}{min}\left\{\left|\frac{\widehat{H}\left(s\right)}{g\left(s\right)}\right|>c\right\}<\infty \right)\hfill \\ & =& \underset{T\to \infty }{lim}\mathbb{P}\left(\underset{s:\left[mT\right]+1\le \left[sT\right]\le T}{sup}\left|\frac{\widehat{H}\left(s\right)}{g\left(s\right)}\right|>c\right)\hfill \\ & =& \mathbb{P}\left(\underset{m\le s\le 1}{sup}\left|\frac{\mathcal{H}\left(s\right)}{g\left(s\right)}\right|>c\right)=\alpha ,\hfill \end{array}$$

(36)

with $\alpha$ denoting the chosen significance level, and $\mathcal{H}\left(s\right)$ short-hand notation for the limit corresponding to the considered monitoring statistic. Considering only positive and bounded weighting functions, also seen in Aue et al. (2012, Assumption 3.6), allows us to derive the required result given above based on the developed asymptotic null behavior of the monitoring statistics and the continuous mapping theorem.

Proposition 2.

Let the data be generated according to (1) and (2) with Assumptions 1 and 2 in place. Furthermore, let long-run covariance estimation be performed under Assumptions 3 and 4 and let lead-lag choices be made as discussed in Choi and Kurozumi (2012). In addition, assume that the weighting function $g\left(s\right)$ is continuous and bounded. Then, there exist critical values $c=c(\alpha ,\widehat{H},g,m,n)$ such that, for any $0<\alpha <1$, it holds that:16

$$\underset{T\to \infty }{lim}\mathbb{P}\left({\tau }_m(\widehat{H},g,c)<\infty \right)=\alpha .$$

(37)

Clearly, the choice of a weighting function $g\left(s\right)$ impacts the performance of monitoring procedures and has to combine two opposing goals: (i) small size distortions under the null hypothesis and (ii) small delays under the alternative hypothesis, that is, detection of a break as soon as possible after the break. The discussion in Chu et al. (1996, Section 3) makes clear that it is, in general, even in more standard regression models, impossible to derive analytically tractable optimal weighting functions (from a certain class of functions), e.g., with respect to minimal expected delay.17

Given the lack of analytical results concerning optimal choices of weighting functions, we have performed a large number of preliminary simulations using a range of candidate weighting functions.18 The starting point of these considerations is Wagner and Wied (2017), who choose the weighting function in relation to the expected value of the monitoring statistics, resulting in $g\left(s\right)=s^3$ in case $D_t=1$ (intercept only) and $g\left(s\right)=s^5$ in case $D_t={[1,t]}^{\prime }$ (intercept and linear trend). In case of a linear trend, we have, in addition, experimented with $g\left(s\right)\in \{1,s^{10},s^{5(0.5+m)},s^{5{(0.85+m)}^2},\sqrt{m}(1+\frac{s-m}{m}),\frac{s}{\sqrt{m}}{\left(\frac{s-m}{s}\right)}^{1/2}\}$, with the last two functions inspired by Horváth et al. (2004).19 It turns out that no weighting, i.e., $g\left(s\right)=1$ does not lead to favorable performance compared to $g\left(s\right)=s^3$ or $s^5$. The function $s^{10}$ is chosen by “extrapolation” of the fact that $s^5$ works better than $s^0=1$. The idea of the third and fourth functions is—merely the result of some experimentation and heuristics—to make the detector more sensitive by increasing the value of the statistic, whilst at the same leaving the critical values effectively unchanged. The effects are to a certain extent as expected, without, however, leading to overall better performance. Taking $s^6={\left(s^3\right)}^2$ or $s^{10}={\left(s^5\right)}^2$ as weighting functions does indeed lead, e.g., to earlier detection times, however, often also to detections in cases when there is no structural change, i.e., these functions lead to larger over-rejections under the null hypothesis. The two functions $s^{5(0.5+m)},s^{5{(0.85+m)}^2}$, where we have also experimented with other powers and values, try to strike a balance between earlier rejections and size distortions. Altogether, however, the simple functions $s^3$ and $s^5$ perform most stably over a variety of configurations, in terms of comparably low over-rejections under the null hypothesis as first priority and short delays in the detection times.20 Finally, the two functions inspired by Horváth et al. (2004), where we have also experimented with different powers, lead to essentially the same null rejection probabilities and size corrected power as, e.g., $s^3$ or $s^5$, but lead to partly substantially bigger delays than the other weighting functions. Therefore, we stick to the weighting functions already used also in Wagner and Wied (2017), i.e., the end point of the considerations is the starting point. As mentioned, it remains an open challenge to make progress on finding optimal weighting functions for the monitoring problem and detectors considered in this paper, or more generally when monitoring cointegrating relationships.

It remains to characterize the asymptotic behavior of the proposed monitoring procedures under the relevant alternatives in our setting: First, the error process ${\left\{u_t\right\}}_{t\in \mathbb{Z}}$ changes its behavior from I(0) to I(1), i.e., it changes to being an integrated process, and, second, there are breaks in (some of) the parameter values. For both cases, we consider the asymptotic behavior against fixed and local alternatives, with the local alternatives having to be specified, as always, in line with the convergence rates of parameter estimation.

Proposition 3.

Let the data be generated for $t=1,\dots ,\lfloor rT\rfloor$ according to (1) and (2) with Assumptions 1 and 2 in place. Furthermore, let long-run covariance estimation be performed under Assumptions 3 and 4 and let lead-lag choices be made as discussed in Choi and Kurozumi (2012). In addition, assume that the weighting function $g\left(s\right)$ is continuous, positive, and bounded. Furthermore, $\widehat{H}\left(s\right)$ again denotes any of the considered monitoring statistics.

(a) 

Let

(i) 

${\left\{u_t\right\}}_{t\in \mathbb{Z}}$ be an $I\left(1\right)$ process from $\lfloor rT\rfloor +1$ onwards, or

(ii) 

${\theta }_1\ne \theta$, with the condition

$$\underset{T\to \infty }{lim}\frac{1}{\sqrt{T}}\sum _{t=\lfloor rT\rfloor +1}^T{D}_t^{\prime }({\theta }_D-{\theta }_{D,1})=\pm \infty$$

(38)

fulfilled.

Then, the monitoring procedures are consistent, i.e., for any $0<c<\infty$, it holds that

$$\underset{T\to \infty }{lim}\mathbb{P}({\tau }_m(\widehat{H},g,c)<\infty )=1.$$

(39)

(b) 

Let

(i) 

${\left\{u_t\right\}}_{t\in \mathbb{Z}}={\left\{u_t^0\right\}}_{t\in \mathbb{Z}}$ for all $t\le \lfloor rT\rfloor$, with ${\left\{u_t^0\right\}}_{t\in \mathbb{Z}}$ satisfying Assumption 2 and

$$u_t=u_t^0+\frac{\delta }{T}\sum _{i=\lfloor rT\rfloor +1}^t{\gamma }_i,$$

(40)

for all $t>\lfloor rT\rfloor$, where ${\left\{u_t^0\right\}}_{t\in \mathbb{Z}}$ and ${\left\{{\gamma }_t\right\}}_{t\in \mathbb{Z}}$ are independent processes and where ${\left\{{\gamma }_t\right\}}_{t\in \mathbb{Z}}$ fulfills an invariance principle with long-run variance ${\omega }_{\gamma }^2>0$ and $\delta >0$;

(ii) 

${\theta }_{D,1}={\theta }_D+G_D^{-1}\left(T\right){\Delta }_{{\theta }_D}$ from $\lfloor rT\rfloor +1$ onwards with $G_D\left(T\right)$ as in Assumption 1 and ${\Delta }_{{\theta }_D}$ fulfilling

$${\int }_r^{1}D{\left(z\right)}^{\prime }dz{\Delta }_{{\theta }_D}\ne 0;$$

(41)

or

(iii) 

${\theta }_{X,1}={\theta }_X+G_X^{-1}\left(T\right){\Delta }_{{\theta }_X}$ from $\lfloor rT\rfloor +1$ onwards with ${\Delta }_{{\theta }_X}\ne 0$.

Then, for any $0<\epsilon \le 1-\alpha$ and the critical value $0<c<\infty$ from Proposition 2, there exists a $\delta =\delta (c,g)$, ${\Delta }_{{\theta }_X}={\Delta }_{{\theta }_X}(c,g)$ or ${\Delta }_{{\theta }_D}={\Delta }_{{\theta }_D}(c,g)$ such that

$$\underset{T\to \infty }{lim}\mathbb{P}({\tau }_m(\widehat{H},g,c)<\infty )\ge 1-\epsilon .$$

(42)

The local asymptotic power (LAP) properties of the procedures, with respect to break type, estimation method, self-normalization, and window size of the moving window detectors, are discussed in some detail in Supplementary Appendix B.21 The LAP results carry over by and large to similar relative performance findings in the finite sample simulations presented in Section 3, e.g., in case of I(1) breaks, a small moving window leads to both highest LAP, as well as highest size-corrected finite sample power. Another example is power differences across estimation methods, with IM-OLS dominated by FM-/D-OLS in terms of LAP and also in finite sample size-corrected power for large sample sizes.

The critical values depend on the specification in several ways: on the specification of the deterministic component, the number of I(1) regressors, the highest power of the single integrated regressor that enters the CPR with higher order powers, the estimation method, and—as discussed in Proposition 2—the detector, the weighting function, the calibration fraction, and for the moving window detectors, in addition the window size. Supplementary Appendix E provides tables containing critical values for the usual deterministic components, i.e., an intercept only and intercept and linear trend, and up to four integrated regressors, of which one regressor enters the model with up to power three. Furthermore, the critical values are available for all considered detectors, FM-/D-OLS and IM-OLS estimation and significance levels of $0.01,0.025,0.05,0.1$. The weighting function $g\left(s\right)$ is set to $s^3$ and $s^5$, depending on the deterministic specification. Critical values are available for a fine grid of the calibration fraction m ranging from $0.1,0.11,\dots ,0.9$. For the moving window detectors, critical values are available for the window sizes equal to $10\%,20\%$, and $30\%$ of the sample size.22

It remains to clarify the “meaning” of m and T for the monitoring procedures (also see p. 967 in Wagner and Wied 2017). It is convenient to interpret T as the sample size including the out-of-sample monitoring period. Let $T_0$ denote the length of the actually available sample—in our application, in Section 4, $T_0=71$ with annual data from 1946–2016. Then, denote $T=T_0+H$, with $H>0$, indicating that out of sample monitoring is intended, and $H=0$, indicating that monitoring takes place on a historic data set, as in Section 4. The fact that the critical values depend on m, and the moving detectors on n, means that a decision has to be made about the length of the calibration period, potentially the length of the moving window and about the out-of-sample monitoring period H prior to the analysis. The latter necessity renders our procedure a closed-end monitoring procedure. The calibration period will be chosen as large as possible (as a sub-sample $1,\dots ,T_C$ of $1,\dots ,T_0$) to increase the precision of the parameter estimates while avoiding the risk of having a structural break in the calibration period. Now, m is given by $m=\frac{T_C}{T_0+H}$. Thus, choosing H larger implies that m is smaller, which in turn implies that the critical values are larger (since they are decreasing in m). This decreases ceteris paribus, despite asymptotic size control, the empirical rejection probabilities under both the null and the alternative. This is the reason why one should choose the monitoring period as short as possible, a calibration period as large as possible, and an out-of-sample monitoring period as short as possible.

Finally, similar to Wagner and Wied (2017), the procedures are consistent also against a variety of other forms of structural changes, with all results following more or less straightforwardly from the construction principle. Of course, the finite sample performance may be relatively poor in some cases.

Remark 2.

The developed monitoring procedures are, in addition to the results provided in Proposition 3, also consistent against the following types of structural change:

(i) 

The process ${\left\{u_t\right\}}_{t\in \mathbb{Z}}$ changes its behavior from I(0) to being a near-integrated process, compare Phillips (1987), from $\lfloor rT\rfloor +1$ onwards. In this case, effectively, under the alternative functionals of Wiener processes will be replaced by functionals of Ornstein-Uhlenbeck processes. The rates of divergence are the same as for the I(1) alternative.

(ii) 

Similarly, consistency also prevails in case ${\left\{u_t\right\}}_{t\in \mathbb{Z}}$ changes its behavior from $\lfloor rT\rfloor +1$ onwards to being fractionally integrated, compared to Davidson and de Jong (2000), with fractional integration parameter $0<f<1/2$. In this case, contrary to item (i), the divergence rate under the alternative changes and depends upon f since, under this alternative, $\frac{1}{T^{1/2+f}}{\sum }_{t=\lfloor rT\rfloor +1}^{\lfloor sT\rfloor }{u}_t$ converges to a fractional Brownian motion. Thus, the smaller f, the more difficult it will be be to detect this form of structural change.

(iii) 

The approach can also be employed for detecting bubbles. In the recent literature, a bubble is often characterized as a period where the behavior of a time series has switched to explosive behavior, compare, e.g., Phillips et al. (2011), Phillips et al. (2015a, 2015b), and Phillips and Shi (2018). Thus, our procedure allows us to detect (the beginning of) a bubble by considering the first difference of the series since, in the absence of a bubble, the first differences are stationary, whereas, in case of explosive behavior, the first differences also exhibit explosive behavior.

(iv) 

In relation to the previous item, with bubbles typically considered to be temporary rather than permanent phenomena, it has to be noted that our procedures will be consistent in detecting episodes of I(1) or explosive behavior, as long as these episodes have asymptotically positive length. In, e.g., the case of only one period under the alternative, it has to hold that this period occurs over a sub-sample of the form $\lfloor r_{1}T\rfloor ,\dots ,\lfloor r_{2}T\rfloor$ with $r_1<r_2$. It is immediate that consistency generalizes to multiple periods of this form.

3. Finite Sample Performance

Under the null hypothesis of no structural change, we consider the same data generating process as Wagner and Hong (2016, Section 3), i.e., we consider a quadratic cointegrating polynomial regression model:

$$y_t={\theta }_{D_0}+{\theta }_{D_1}t+{\theta }_{X_{1,1}}{x}_t+{\theta }_{X_{1,2}}{x}_t^2+u_t,$$

(43)

with the errors $u_t$ and $\Delta x_t=v_t$ generated as:

$$\begin{array}{cc}\hfill u_t& ={\rho }_1{u}_{t-1}+e_{1,t}+{\rho }_2{e}_{2,t},u_0=0,\hfill \\ \hfill v_t& =e_{2,t}+0.5e_{2,t-1},\hfill \end{array}$$

(44)

where ${(e_{1,t},e_{2,t})}^{\prime }\sim \mathcal{N}(0,I_2)$. The parameter ${\rho }_1$ controls the extent of serial correlation in ${\left\{u_t\right\}}_{t\in \mathbb{Z}}$ and is set to ${\rho }_1=1$ after $\lfloor rT\rfloor$ under the alternative of I(1) errors, whereas ${\rho }_2$ controls the extent of regressor endogeneity. The parameter values are ${\theta }_{D_0}={\theta }_{D_1}=1$, ${\theta }_{X_{1,1}}=5$ and ${\theta }_{X_{1,2}}=-0.3$, with the values for ${\theta }_{X_{1,1}}$ and ${\theta }_{X_{1,2}}$ inspired by the FM-OLS EKC coefficient estimates for Austria (see Wagner 2015).23

We provide simulation results for $T\in \{200,500\}$ and ${\rho }_1={\rho }_2\in \{0,0.3,0.6,0.9\}$.24 For the moving window and self-normalized moving window detectors, we use window sizes $n\in \{0.1,0.2,0.3\}$. As indicated in the previous section, long-run covariance estimation is performed with the Bartlett kernel with bandwidth chosen according to Newey and West (1994). For D-OLS, estimation leads and lags choices are performed using the AIC-type criterion of Choi and Kurozumi (2012). The number of replications is $\mathrm{10,000}$ throughout. All monitoring decisions are performed at the nominal $5\%$ significance level.

We start the analysis by considering empirical null rejection probabilities. In doing so, we vary the calibration fraction over a grid of 81 values in the range $m=0.1,0.11,\dots ,0.9$. Two main observations that allow us to zoom in subsequently, mostly on the self-normalized detectors, emerge. Figure 1 clearly shows that with respect to null rejection probabilities, the detectors are separated in two groups, with the better performance offered by self-normalization. Figure 2 shows that, when using moving window detectors, the choice of the window size has no visible impact on null rejection probabilities.25

To assess the performance improvement that can be realized by self-normalization in some detail, Figure 3 compares ${\widehat{H}}_{mov}^{m,0.1}$ and ${\widehat{H}}_{mov,sn}^{m,0.1}$ for $T=200$.26 This figure illustrates also other more generally observed patterns: First, using IM-OLS for parameter estimation leads to the smallest over-rejections under the null hypothesis. Often, and particularly for small values of m, D-OLS estimation leads to the largest over-rejections. The differences across estimation methods widen for increasing ${\rho }_1,{\rho }_2$, with D-OLS most strongly negatively affected. In addition, with the exception of D-OLS, self-normalization attenuates the detrimental impact of increasing ${\rho }_1,{\rho }_2$. An increasing sample size, of course, reduces over-rejections by and large.

We turn to size-corrected power.27 For brevity, the main text focuses on size-corrected power against I(1) breaks, i.e., the situation where ${\left\{u_t\right\}}_{t\in \mathbb{Z}}$ changes its behavior from I(0) to I(1) after $\lfloor rT\rfloor$. The other two breaks dealt with in Proposition 3, trend and slope breaks, with changes in ${\theta }_D$ or ${\theta }_X$, respectively, are discussed in Supplementary Appendix C.28

Size-corrected power simulations are performed for $m,r\in \{0.25,0.5,0.75\}$, which includes, therefore, cases where $r<m$, i.e., where a break occurs in the calibration period. We report results for all nine detectors considered, for all three estimation methods, ${\rho }_1={\rho }_2=0.3$ until $\lfloor rT\rfloor$ and $T=200$ in

Table 1. Size-corrected power against I(1) breaks for $T=200$ and ${\rho }_1={\rho }_2=0.3$.

${\mathit{\rho }}_1={\mathit{\rho }}_2=0.3$			$\mathit{m}=0.25$				$\mathit{m}=0.5$				$\mathit{m}=0.75$
		r	$0.25$	$0.5$	$0.75$		$0.25$	$0.5$	$0.75$		$0.25$	$0.5$	$0.75$
	FM		0.19	0.07	0.05		0.32	0.64	0.14		0.34	0.46	0.80
${\widehat{H}}^m$	D		0.13	0.06	0.05		0.30	0.59	0.11		0.34	0.46	0.76
	IM		0.09	0.06	0.05		0.26	0.47	0.07		0.32	0.41	0.66
	FM		0.19	0.07	0.05		0.32	0.64	0.14		0.33	0.45	0.79
${\widehat{H}}_d^m$	D		0.13	0.06	0.05		0.30	0.59	0.11		0.33	0.45	0.76
	IM		0.09	0.06	0.05		0.26	0.47	0.07		0.32	0.41	0.66
	FM		0.20	0.07	0.05		0.17	0.61	0.14		0.16	0.26	0.77
${\widehat{H}}_{sn}^m$	D		0.14	0.06	0.05		0.20	0.56	0.11		0.19	0.30	0.73
	IM		0.11	0.05	0.05		0.19	0.49	0.08		0.21	0.31	0.66
	FM		0.18	0.07	0.05		0.32	0.66	0.17		0.34	0.46	0.80
${\widehat{H}}_{mov}^{m,0.1}$	D		0.12	0.06	0.05		0.29	0.61	0.13		0.34	0.46	0.77
	IM		0.09	0.05	0.05		0.25	0.48	0.09		0.31	0.40	0.67
	FM		0.18	0.07	0.05		0.32	0.65	0.15		0.34	0.46	0.80
${\widehat{H}}_{mov}^{m,0.2}$	D		0.13	0.06	0.05		0.30	0.60	0.12		0.34	0.46	0.76
	IM		0.09	0.06	0.05		0.26	0.48	0.08		0.32	0.41	0.66
	FM		0.18	0.07	0.05		0.32	0.64	0.14		0.34	0.46	0.80
${\widehat{H}}_{mov}^{m,0.3}$	D		0.13	0.06	0.05		0.30	0.59	0.11		0.34	0.46	0.76
	IM		0.09	0.06	0.05		0.26	0.47	0.08		0.32	0.41	0.66
	FM		0.20	0.07	0.05		0.17	0.63	0.18		0.16	0.26	0.77
${\widehat{H}}_{mov,sn}^{m,0.1}$	D		0.14	0.06	0.05		0.20	0.58	0.13		0.20	0.30	0.74
	IM		0.10	0.05	0.05		0.18	0.49	0.10		0.20	0.30	0.66
	FM		0.20	0.07	0.05		0.17	0.62	0.15		0.16	0.26	0.77
${\widehat{H}}_{mov,sn}^{m,0.2}$	D		0.14	0.06	0.05		0.20	0.57	0.11		0.19	0.30	0.73
	IM		0.10	0.05	0.05		0.19	0.49	0.09		0.21	0.31	0.66
	FM		0.20	0.07	0.05		0.17	0.62	0.14		0.16	0.26	0.77
${\widehat{H}}_{mov,sn}^{m,0.3}$	D		0.14	0.06	0.05		0.20	0.56	0.11		0.19	0.30	0.73
	IM		0.10	0.05	0.05		0.19	0.49	0.08		0.21	0.31	0.66

.29 The following observations emerge: Grosso modo, highest size-corrected power is achieved by the moving window detector with $n=0.1$, either with or without self-normalization; with as reported above self-normalization leading to smaller over-rejections under the null hypothesis. The differences in size-corrected power are typically relatively or even very small. This, in conjunction with the performance improvement of self-normalization under the null hypothesis, makes the case in favor of self-normalization clear. In line with standard asymptotic theory concerning estimator efficiency, size-corrected power is typically lower for IM-OLS than for FM-/D-OLS, with the difference between the latter two often rather small, also in line with asymptotic theory.30 These results, of course, have to be seen in conjunction with the smaller over-rejections of IM-OLS under the null hypothesis.

Next, consider the impact of m and r on size-corrected power. In case $r<m$, i.e., when a break occurs in the calibration period, size-corrected power is often low, which is related to inconsistency of parameter estimation due to the structural break in the calibration period. It turns out that, in this case, the self-normalized detectors lead to lower size-corrected power than the standardized detectors; see Table 1. For $m\le r$, size-corrected power behavior is as expected, i.e., for fixed m, size-corrected power decreases with increasing r, and, for fixed r, size-corrected power increases with increasing calibration period m. The case $m=r$ leads to the highest size-corrected power. In addition, as expected, increasing T leads to higher size-corrected power, whereas increasing ${\rho }_1,{\rho }_2$ leads to lower size-corrected power. For trend breaks—investigated in some more detail in Supplementary Appendix C—, many observations are qualitatively similar. One difference is that, in case of trend breaks, the detectors ${\widehat{H}}^m$ or ${\widehat{H}}_d^m$ lead to highest size-corrected power in some configurations. This finding, however, has to be considered in light of the larger over-rejections exhibited by these two detectors under the null hypothesis. Slope breaks lead to very similar results as I(1) breaks.

It remains to investigate the detection times and delays. We again consider the case of I(1) breaks in the main text, for $T=200$ and for ${\rho }_1={\rho }_2=0.3$ in Figure 4.31 With respect to m and r we only display results for $m=0.5$ and $r=0.5$; additional configurations leading to qualitatively very similar results are displayed in Figures S21–S26 in Supplementary Appendix C. The figures display the results in the form of box-whiskers plots—for completeness—for all nine considered detectors. The numbers below the detector labels indicate the corresponding null rejection probabilities, i.e., size-corrected power, given in Table 1. Thus, the different box-whiskers plots are based on different numbers of replications because of different numbers of rejections across different detectors, sample sizes, and $\rho$ values.

One relatively clear observation that emerges from the figures, both in the main text and in Supplementary Appendix C, is that the choice of the estimation method does not exhibit major impacts on detection times and delays. Given the standard asymptotic properties of the estimators, as expected, IM-OLS leads to slightly larger delays than FM-OLS and D-OLS in many cases. The moving window detectors lead to the shortest delays, both standardized and self-normalized, with the best performance achieved with $n=0.1$. One interesting observation is that, for the delay, the choice of the window size does matter and, in fact, exerts bigger influence on the results than the choice of standardizing or self-normalizing. Furthermore, an increasing sample size leads to a—ceteris paribus—more concentrated distribution of the estimated detection times (based on a larger number of observations) but does not, throughout, lead to smaller average delays.

The results of the finite sample simulations can essentially be summarized as follows: First, when comparing self-normalized monitoring statistics with their standardized versions, self-normalization leads to smaller over-rejections, without having any systematic or sizeable negative impacts on size-corrected power (for $m\le r$), and without leading to larger delays when compared to standardization. Second, within the group of self-normalized detectors, the moving window detector with $n=0.1$ leads almost throughout to highest size-corrected power without detrimental effects on null rejection probabilities. For the moving window detectors, size-corrected power decreases with increasing window size. For all moving window detectors, size-corrected power is at least as good or higher than for the expanding window detector. Some exceptions occur in case of trend breaks, in which case the expanding window detector sometimes outperforms the moving window detectors in terms of size-corrected power. In addition, the window size exhibits some impact on the delays, with the moving window detector with $n=0.1$ leading to the shortest delays. Third, with respect to estimator choice, the usual trade-off between IM-OLS, on the one hand, and FM-/D-OLS, on the other, occurs. IM-OLS leads to lower over-rejections at the cost of lower size-corrected power and slightly larger delays than FM-OLS, with D-OLS outperformed by FM-OLS in terms of larger over-rejections by D-OLS under the null hypothesis. Size-corrected power and detection times based on D-OLS are by and large very similar compared to those based on FM-OLS. Self-normalization in conjunction with IM-OLS is particularly beneficial for small samples, as in this case no long-run covariance estimates are required. The poor performance of FM-OLS compared to IM-OLS for small samples can be traced back to long-run covariance estimation required for FM-OLS parameter estimation, even when using self-normalized detectors. For larger sample sizes, with better properties of long-run covariance estimation, the asymptotic efficiency advantage of FM-OLS over IM-OLS becomes visible, particularly with respect to higher size-corrected power.

Altogether, we recommend to use the self-normalized moving window detector with $n=0.1$, i.e., ${\widehat{H}}_{mov,sn}^{m,0.1}$. The choice of the estimator is dictated by the trade-off between null rejection probabilities on target and high size-corrected power.

4. The Environmental Kuznets Curves (EKCs) for Carbon and Sulfur Dioxide Emissions

We now apply the developed monitoring procedures to EKCs for carbon and sulfur dioxide (CO${}_2$ and SO${}_2$) emissions. We commence from a sample of 18 industrialized countries, listed in

Table 2. List of countries included in the empirical analysis. The sample range is 1946–2016. Italic country names indicate that the augmented Dickey-Fuller test rejects the unit root null hypothesis for log GDP per capita on the calibration period (1946–1973) at the 10% level, and bold country names indicate rejections at the 5% level. Intercept and linear trend are included.

Australia	Austria	Belgium	Canada	Denmark	Finland
France	Germany	Italy	Japan	New Zealand	Norway
Portugal	Spain	Sweden	Switzerland	United Kingdom	United States

, over the period 1946–2016.32 The highest polynomial degree we consider is the cubic cointegrating polynomial relationship, i.e.,

$$y_t={\theta }_{D_0}+{\theta }_{D_1}t+{\theta }_{X_{1,1}}{x}_t+{\theta }_{X_{1,2}}{x}_t^2+{\theta }_{X_{1,3}}{x}_t^3+u_t,$$

(45)

with $x_t$ the logarithm of real per capita GDP and $y_t$ the logarithm of per capita CO₂ or SO₂ emissions. As will be seen below, over the calibration period 1946–1973, for many countries, in fact, a cointegrating quadratic or even cointegrating linear relationship prevails, i.e., is not rejected.

The GDP and CO${}_2$ emissions data are similar to those used in Wagner et al. (2020). The main difference is the extension of the sample period from 2013 to 2016; additionally, there are some marginal changes in the GDP data in the newer vintage. More precisely, the GDP (in 2011 prices) and population data stem from the Maddison project database (in the 2018 version of Bolt et al. 2018). The CO${}_2$ emissions data—which cover CO${}_2$ emissions from fossil fuel usage—are taken from the Carbon Dioxide Information and Analysis Center (CDIAC) of the U.S. Department of Energy; see Boden et al. (2018). The SO${}_2$ emissions have been combined from two sources. The data for 1946–2005 are from the NASA Socioeconomic Data and Applications Center (SEDAC); see Smith et al. (2011). The data for the period 2006–2016 are from the OECD (2020).33

With the developed methods resting upon a calibration period, a corresponding choice has to be made. We choose the calibration period 1946–1973, reflecting that the first oil price shock of 1974 is considered a major event for changes in energy consumption patterns.34 In our notation, this amounts to $m=28/71\approx 0.4$. Given this choice, the first step is to perform the CPR modeling cycle for the calibration period. The augmented Dickey and Fuller (1981) test results, obtained on the calibration period for the null hypothesis of a unit root in log real per capita GDP against the alternative of (linear-)trend stationarity, are contained in the country list Table 2.35 These results lead to the following twelve countries being included in the subsequent analysis (based on a 5% significance level): Australia, Belgium, Canada, Denmark, Finland, Italy, Japan, Portugal, Spain, Sweden, the United Kingdom (UK), and the United States (U.S.). The next step is to test for the prevalence of a CPR relationship over the monitoring period—for both CO${}_2$ and SO${}_2$ emissions for the countries with I(1) log real per capita GDP. Given that the polynomial degree is ex ante unclear, we perform a testing sequence with polynomial degrees ranging from three (the cubic specification) to one (linear cointegration).36 The deterministic specification consists of intercept and linear trend throughout. The lowest degree polynomial for which a CPR relationship is not rejected is considered as starting point for monitoring in the following subsections; see

Table 3. Minimal polynomial degrees for cointegrating environmental Kuznets curves (EKCs) over the calibration period 1946–1973.

	Polynomial Degree
Country	CO2	SO2
Australia	1	1
Belgium	1	1
Canada	1	1
Denmark	1	2
Finland	2	2
Italy	1	1
Japan	1	2
Portugal	1	1
Spain	1	1
Sweden	1	2
United Kingdom	1	1
United States	2	3

.

The most striking feature of Table 3 is that, in most countries—for CO${}_2$ even more than for SO${}_2$—, a cointegrating linear relationship appears to be present over the period 1946–1973. This seems to be at odds, at first sight, with the EKC hypothesis of an inverted U-shaped relationship. However, these results reflect the fact that, until the early 1970s, per capita GDP and per capita emissions developed very similarly on a “log-linear extension path”; see the scatter plots in Figures S27 and S28 in Supplementary Appendix D. Only starting in the mid 1970s, the oil price shock, as well as environmental legislation—particularly with respect to SO${}_2$ and “acid rain”—, lead to a decoupling of the two quantities, to a certain extent, in many countries.37

Given these preliminary results, we turn in the following two subsections to monitoring, using the moving window detector ${\widehat{H}}_{mov,sn}^{m,0.1}$ with window size $n=0.1$, i.e., seven years.38

4.1. Results for Carbon Dioxide Emissions

The detection times for the twelve considered countries using, with the exceptions of Finland and the United States, where the quadratic specification is estimated, the linear specification are displayed in

Table 4. Detection times when monitoring CO${}_2$ emissions using the detector ${\widehat{H}}_{mov,sn}^{m,0.1}$ and both Fully Modified Ordinary Least Squares (FM-OLS) and Integrated Modified Ordinary Least Squares (IM-OLS). The column p indicates the polynomial degree, the calibration period is 1946–1973, the monitoring period is 1974–2016. The nominal significance level is 5%.

Country	p	${\widehat{\mathit{H}}}_{\mathbf{mov},\mathbf{sn}}^{\mathit{m},0.1}$
		FM-OLS	IM-OLS
Australia	1	1993	2001
Belgium	1	1988	1992
Canada	1	∞	∞
Denmark	1	1991	2011
Finland	2	1989	1990
Italy	1	1981	1982
Japan	1	1982	1980
Portugal	1	1998	∞
Spain	1	∞	∞
Sweden	1	1982	1983
United Kingdom	1	1984	1987
United States	2	1988	1992

.

With the exception of Canada and Spain for all countries the linear relationship found over the calibration period appears to break down. Furthermore, Portugal is a “mixed” case, with unsurprisingly similar behavior to Spain in many ways, as discussed below. In line with the simulation results, the detection times are—with the exception of Japan—earlier with FM-OLS than with IM-OLS, with the differences often just a few years. Mostly, breaks are dated in the 1980s or early 1990s, with a few exceptions: Australia (IM-OLS: 2001), Denmark (IM-OLS: 2011), Portugal (FM-OLS: 1998; IM-OLS: no break, which is the only difference in finding a break point or not between the two methods). This can be interpreted, admittedly unsurprisingly, as strong evidence against a continued log-linear co-movement, i.e., aligned expansion, of log per capita GDP and CO${}_2$ emissions from the—given the delays seen in the simulation section at latest—early 1990s onwards. Or, the other way around, these findings are indicative of curvature picking up in CO${}_2$ EKC-type relationships in the 1980s, in most of the considered countries, with the exception of Canada, and potentially Portugal and Spain, which is discussed further below.39

As an illustration, Figure 5 displays the mechanics of the monitoring procedure for Finland.40 The upper graph displays the FM- and IM-OLS residuals over both the calibration and the monitoring period. By definition, since an intercept is included, the residuals have zero mean over the calibration period and then turn systematically (more and more) negative thereafter on the monitoring period. This, by definition, means that the estimated—in case of Finland quadratic—relationship systematically over-predicts actual log per capita CO${}_2$ emissions, i.e., the slope of an actual relationship between output and emissions is estimated as too high; also see Figure 6. The statistical monitoring procedures need to collect enough signals until a break is detected; in this example, this takes until 1989 for FM-OLS and until 1990 for IM-OLS, which is clearly too late. Nevertheless, it may be interesting to note that Finland was severely adversely affected by the collapse of the Soviet Union in the early 1990s, which, in fact, for Finland, led to a deeper recession than the Great Depression in the 1930s.

Figure 6 shows the fitted values obtained when estimating the quadratic EKC with FM-OLS over the calibration sample, the sample ending prior to the detected break-point 1988 and when estimation takes place over the full sample.41 The dashed red line shows that, in fact, more or less immediately after the oil price shock, the relationship between GDP and CO${}_2$ emissions appears to have changed, with the detector taking about 15 years to collect enough signals to declare the null hypothesis rejected at the 5% significance level.42 Parameter estimation until 1988 “catches the turn” of the mid 1970s and leads to good fit. However, using the full sample for parameter estimation of the quadratic EKC leads to the best fit, particularly for the period after 1989. This leads to some support for the interpretation that maybe the relationship is after all quadratic and one needs to have a sample that is (i) larger and (ii) covers the inverted U behavior to be able to estimate the relationship with sufficient precision; compare Footnote 37.43

To complete the analysis for CO${}_2$ emissions, we now turn to the countries for which no structural change is indicated by the monitoring procedure, Canada, the mixed case Portugal, and Spain.

Table 5. FM-OLS and IM-OLS estimation results for a cointegrating linear relationship between log per capita GDP and log per capita CO${}_2$ emissions for Canada, Portugal and Spain for the calibration period 1946–1973 (left panel) and the full sample period 1946–2016 (right panel). Italic entries indicate significance of coefficients at the 10% level, and bold entries significance of coefficients at the 5% level.

	1946–1973		1946–2016
	${\mathit{\theta }}_{{\mathit{D}}_{\mathbf{1}}}$	${\mathit{\theta }}_{{\mathit{X}}_{\mathbf{1},\mathbf{1}}}$	${\mathit{\theta }}_{{\mathit{D}}_{\mathbf{1}}}$	${\mathit{\theta }}_{{\mathit{X}}_{\mathbf{1},\mathbf{1}}}$
Canada
FM-OLS	−0.056	2.841	−0.027	1.661
IM-OLS	−0.059	2.990	−0.028	1.714
Portugal
FM-OLS	0.000	1.003	−0.008	1.341
IM-OLS	0.008	0.869	0.001	1.156
Spain
FM-OLS	−0.023	1.519	−0.033	1.818
IM-OLS	−0.039	1.789	−0.036	1.911

displays the estimated coefficients for these three countries when the cointegrating (linear) relationships between log per capita GDP and log per capita CO${}_2$ emissions are estimated over the calibration sample and over the whole sample period.

The results reported in Table 5 are a bit mixed, and no uniform pattern across all three countries emerges, with some of the smoke clearing when considering Figure 7 below. For Canada, the slope coefficient corresponding to log per capita GDP becomes much smaller when estimated over the full sample, but the trend coefficient also becomes less negative. For Portugal, the trend coefficient is not significant in both periods, but turns from non-negative point estimates in the calibration period to smaller, and for FM-OLS negative, values in the full period. The slope coefficient, however, is bigger, when estimated over the full sample. For Spain, the trend coefficient does become more negative and significant when estimated over the full sample period and the slope coefficient becomes bigger. Thus, qualitatively, the results are similar for the two countries on the Iberian peninsula, which not only share many economic similarities but also political similarities, with the end of the Franco and Salazar regimes in these countries in the mid 1970s.

Probably more informative, Figure 7 displays similar results for these three countries as displayed for Finland in Figure 6.44 The message is clear, particularly in comparison with Figure 6. In case no structural break is detected, estimation over the larger sample does lead—unsurprisingly—to better fit, but the differences in fit are, compared to the differences observed in Figure 6, minor. From a standard inverted U quadratic EKC perspective, the excellent fit obtained for these three countries with the linear specification is to a certain extent surprising—and, which might be bad news, of course, there is no turning point.

Note that, when considering the full sample period 1946–2016, we find a cointegrating EKC relationship for all nine countries, for which monitoring detects a structural break.45 For about a third of the countries, the polynomial degree of the EKC appears to be larger for the full sample than for the calibration period, supporting the notion of more curvature. For most countries, the evidence, with an unchanged minimal polynomial degree, points, however, more towards structural change of the parameters for a given form of the relationship. For Finland and the U.S., it appears that the polynomial degree is lower over the full sample than for the calibration period, albeit the coefficients to squared log per capita GDP are significant when estimated over the full sample, which leads to an unclear picture. Altogether, the results in this CO${}_2$ subsection indicate that the procedures work even on such small samples, albeit, of course, for many applications the sizeable delay is most certainly problematic.

4.2. Results for Sulfur Dioxide Emissions

We now turn to monitoring the EKC for SO${}_2$ emissions, with the results displayed in

Table 6. Detection times when monitoring SO${}_2$ emissions using the detector ${\widehat{H}}_{mov,sn}^{m,0.1}$ and both FM-OLS and IM-OLS. The column p indicates the polynomial degree, the calibration period is 1946–1973, the monitoring period is 1974–2016. The nominal significance level is 5%.

Country	p	${\widehat{\mathit{H}}}_{\mathbf{mov},\mathbf{sn}}^{\mathit{m},0.1}$
		FM-OLS	IM-OLS
Australia	1	∞	∞
Belgium	1	1985	1985
Canada	1	1980	1981
Denmark	2	1997	2006
Finland	2	1988	1991
Italy	1	1982	1982
Japan	2	1991	2004
Portugal	1	2015	∞
Spain	1	2005	2012
Sweden	2	1985	1988
United Kingdom	1	1983	1983
United States	3	∞	∞

for the detector ${\widehat{H}}_{mov,sn}^{m,0.1}$ for both FM-OLS and IM-OLS.46

Qualitatively, the results for SO${}_2$ are similar to the results for CO${}_2$ discussed in the previous subsection. With the exception of Australia, the United States and the mixed case (as for CO${}_2$ emissions) Portugal monitoring indicates structural breaks. Again, in line with the simulation findings, the FM-OLS break points are in no case later than the IM-OLS break points, with equal break dates from both estimators for three countries. This similarity of detected break points can be tentatively interpreted as evidence for underlying structural changes in the relation between economic activity and emissions in—or due to the delays prior to—the 1980s, potentially as a consequence of tighter environmental legislation.47

Table 7. FM-OLS and IM-OLS estimation results for a cointegrating linear relationship between log per capita GDP and log per capita SO${}_2$ emissions for Australia and Portugal, and a cointegrating cubic relationship for the United States, for the calibration period 1946–1973 (left panel) and the full sample period 1946–2016 (right panel). Italic entries indicate significance of coefficients at the 10% level, and bold entries significance of coefficients at the 5% level.

	1946–1973				1946–2016
	${\mathit{\theta }}_{{\mathit{D}}_{\mathbf{1}}}$	${\mathit{\theta }}_{{\mathit{X}}_{\mathbf{1},\mathbf{1}}}$	${\mathit{\theta }}_{{\mathit{X}}_{\mathbf{1},\mathbf{2}}}$	${\mathit{\theta }}_{{\mathit{X}}_{\mathbf{1},\mathbf{3}}}$	${\mathit{\theta }}_{{\mathit{D}}_{\mathbf{1}}}$	${\mathit{\theta }}_{{\mathit{X}}_{\mathbf{1},\mathbf{1}}}$	${\mathit{\theta }}_{{\mathit{X}}_{\mathbf{1},\mathbf{2}}}$	${\mathit{\theta }}_{{\mathit{X}}_{\mathbf{1},\mathbf{3}}}$
Australia
FM-OLS	0.043	−0.819			−0.042	2.541
IM-OLS	0.042	−0.770			−0.043	2.550
Portugal
FM-OLS	−0.005	1.047			−0.104	3.491
IM-OLS	0.002	0.873			−0.092	3.215
United States
FM-OLS	−0.013	−1,831.120	183.982	−6.158	−0.133	−305.981	31.361	−1.052
IM-OLS	−0.005	−2,190.198	220.246	−7.380	−0.133	−61.963	7.512	−0.275

displays the estimation results for Australia, Portugal and the U.S., i.e., the three countries where no structural break was detected, or—to be precise—where, for Portugal, only FM-OLS estimation leads to a break being detected. The results lead to some interesting observations. For Australia, estimation on the calibration period leads to (with all coefficients not significantly different from zero) positive trend coefficient and negative slope coefficient in a cointegrating linear relationship.48 This is surprising to a certain extent, as a negative trend slope is typically considered to capture autonomous energy efficiency increases, admittedly more important for CO${}_2$ emissions rather than SO${}_2$ emissions. The expected signs, with the trend coefficient negative and the slope coefficient positive (and significant), emerge only over the full sample. For Portugal, the coefficient signs are as expected, with the trend coefficient only significant for the full sample. The slope coefficient to log per capita GDP becomes much larger when estimated over the full sample. The results for the cubic specification estimated for the U.S. are hard to interpret with the enormous fluctuations in the coefficient values. However, the implied turning points do change in expected ways, meaning that the first one becomes smaller and the second one larger, e.g., for FM-OLS and estimation over the short sample the two turning points are at about US$ 16,000 and 27,000, whereas full sample estimation leads to turning points at about US$ 5400 and 79,000.

The results in Figure 8, however, indicate that indeed the relationship specified over the calibration period, when estimated over the full sample, leads to very good fit, either with a cointegrating linear relationship for Australia and Portugal or with a cubic relationship for the U.S.49

Considering the full sample period leads—similar to CO${}_2$ emissions—to a CPR relationship for all nine countries where a structural break has been detected. In addition, similar to CO${}_2$, the evidence is mixed with respect to either change in the parameters or a change in the minimal polynomial degree of a cointegrating EKC for SO${}_2$ emissions.50

5. Summary and Conclusions

The paper extends the residual-based monitoring procedure for cointegrating relationship developed in Wagner and Wied (2017) in two dimensions. First, in addition to the detector studied in that paper, we consider a number of detectors that consider two aspects not considered in detail in the earlier paper: self-normalization and moving window detectors. Second, the approach is extended from cointegrating linear to cointegrating polynomial regressions (CPRs).

The starting point is the extension of the Shin (1994) test statistic for the null hypothesis of linear cointegration to the CPR setting (see, e.g., Wagner 2020). The full sample test statistic is then used as the basis for developing monitoring procedures, similar to Wagner and Wied (2017) for the case of cointegrating linear regressions. As also usual in the cointegrating regression literature, the regressors are allowed to be endogenous and the stationary errors are allowed to be serially correlated. Consequently, parameter estimation is based on modified OLS estimators, to be precise, FM-OLS, D-OLS, and IM-OLS, in their versions adapted to the CPR setting to allow for the construction of nuisance parameter-free monitoring statistics by scaling out a scalar long-run variance. However, even the usage of suitably modified least squares regressors is not, in the general case, sufficient to arrive at nuisance parameter-free monitoring statistics. In the CPR case, as mentioned in Remark 1, additionally, full design is required to arrive at nuisance parameter-free test statistics, either by standardization or by self-normalization. By construction, and as in the cointegrating linear case, the limiting distributions of the monitoring statistics coincide for FM-OLS and D-OLS, whilst IM-OLS leads to different limiting distributions for the monitoring statistics.

Both self-normalization, as well as the consideration of moving window detectors, turn out to be beneficial. The combination of these two new elements leads to the, grosso modo, best performance. More precisely, this means that self-normalization leads to smaller over-rejections under the null hypothesis, without leading to strong disadvantages in terms of either size-corrected power or detection delays. The performance differences between a standardized and a self-normalized detector reflect the (well-known small-sample) problems associated with long-run covariance estimation. Using a moving window rather than an expanding window also contributes positively to performance. The idea behind moving window detectors is to reduce the impact of pre-break observations on the monitoring statistics in the post-break period, and the simulation evidence indicates that this is indeed what happens. The finite sample performance differs to a certain extent across the three estimators. IM-OLS leads to, with the differences more pronounced for small samples, smaller over-rejections under the null hypothesis than FM-OLS, which in turn leads to lower over-rejections than D-OLS. This, however, comes in conjunction with slightly smaller size-corrected power and slightly bigger delays of IM-OLS compared to FM-OLS. The evidence is a bit mixed, and the differences are often a bit too small to give an unequivocal recommendation concerning the usage of either FM-OLS or IM-OLS.

The monitoring statistics, precisely (in the main text only) the best performing variant using self-normalization in conjunction with a short moving window, are then used to investigate the stability of cointegrating (polynomial) EKCs for both carbon and sulfur dioxide emissions for a sample of twelve countries over the period 1946–2016, with a calibration period ranging from 1946–1973, i.e., until prior to the first oil price shock. One interesting observation is that, over the calibration period, especially for CO${}_2$ for the majority of countries, in fact, a cointegrating linear relationship is present, in line with effectively “balanced” growth of economic activity and emissions (from burning fossil fuel), until the mid 1970s. The monitoring procedures indicate, for both CO${}_2$ and SO${}_2$, structural changes in nine out of twelve countries. Of course, especially for the cases where a linear cointegrating relationship prevails over the calibration period, this may not be too big a surprise. No structural break is detected for Canada, Portugal, and Spain for CO${}_2$ emissions, and Australia, Portugal, and the U.S. for SO${}_2$ emissions. For the country-pollutant combinations where no structural break is detected, the type of cointegrating relationship specified over the calibration period leads to very good fit also over the full sample period, and, of course, even better fit when the parameters are estimated over the full sample. This is not the case for countries where breaks are detected, indicating that structural breaks are indeed detected when present. Altogether, with the caveat of, in some cases, late detection times, remember the delays illustrated in the simulation study, the monitoring procedures lead to “plausible” findings in our application.

Future work needs to address inter alia the following issues left open in this paper: First, of course, large delays are problematic, even in an ex-post exercise one would want to date the breaks as precisely as possible, e.g., to gauge the time it takes until behavioral or legislative changes translate into changes in the economic activity-emissions nexus. For certain applications (when, e.g., legislation comes into place in a group of countries at the same time), it may be possible to rely also upon the cross-sectional dimension to date break points. The seemingly unrelated cointegrating polynomial regression analysis put forward in Wagner et al. (2020) or the classical panel EKC analysis of de Jong and Wagner (2018) might serve as starting points for developing monitoring statistics for small and large cross-sectional dimension, respectively. Second, if one wants to use the monitoring procedures in a real time manner—rather than for a historical analysis as in this paper—, it may be important to consider extending end-of-sample cointegrating break point tests to the setting considered here (one such test for the cointegrating linear case is discussed in Andrews and Kim 2006). Extensions to more general settings are, however, not obvious. Third, it may be interesting to monitor and detect structural breaks in the other direction, i.e., from a spurious regression to a cointegrating polynomial regression. This could be achieved, e.g., by suitably extending the results of Sakarya et al. (2015) from cointegrating linear to cointegrating polynomial regressions.