Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
100,597 characters · 13 sections · 66 citation commands
Autoregressive Wild Bootstrap Inference for Nonparametric Trends
JEL classifications: C14, C22.\\ Keywords: autoregressive wild bootstrap, nonparametric estimation, time series, simultaneous confidence bands, trend estimation.
The analysis of smoothly evolving trends is of interest in many fields such as economics and climatology. For instance, trend analysis in environmental variables is of major importance, as it can often directly be linked to climate change. Given that trends typically do not evolve in a linear way, fitting linear trends to the data does not uncover actual change accurately. Instead, one would like to use more flexible trend models that avoid making parametric assumptions on the form of the trend. A large body of statistical and econometric research therefore focuses on nonparametric trend modeling and estimation.
In addition, when modeling trends in temperature and emission data, researchers also need to take into account that serial dependence and heteroskedasticity may be present in the data, see for example FV and MV who study parametric trend modeling in temperature series in the presence of serial dependence. Bootstrap methods provide an easy and powerful way to account for heteroskedasticity and autocorrelation. Buhlmann shows the validity of the autoregressive sieve bootstrap for nonparametric trend modeling under general forms of dependence. Neumann uses a wild bootstrap method to achieve robustness to heteroskedasticity for a similar model.
The wild bootstrap approach is also suitable for dealing with missing data, as advocated for example by Shao. It does not require any resampling and therefore the missing data points can keep their original date in a bootstrap sample. In particular in climatology, this feature of the wild bootstrap offers an important benefit over other methods, since there is no need of imputing missing data points. Missing data are a prominent feature in many climatological datasets due to instrument failure or adverse weather and measurement conditions; for instance, in our application, data are missing when cloud cover prevents measurements from being taken. The wild bootstrap, however, relies on independence of the error terms, which is a situation rarely encountered in practice. To relax this strong assumption, dependent versions of wild bootstrap methods have been proposed - see Shao, LN and SU - but not in the context of nonparametric trend estimation. Moreover, so far no theory exists on the validity of such wild bootstrap methods in the presence of serial dependence, heteroskedasticty and missing data. In this paper we address this issue and propose an autoregressive wild bootstrap method that provides valid inference under general conditions for nonparametric trend modeling.
Next to the basic pointwise confidence intervals, we also study simultaneous confidence bands, which are often more informative about trend shapes than pointwise confidence intervals. Research questions, like whether upward trends are present over a certain period of time, should be addressed with simultaneous confidence bands as they involve multiple points in time at once. WuZhao derive such bands for the nonparametric trend model that have asymptotically correct coverage probabilities, but do not consider bootstrap methods. Buhlmann proposes sieve bootstrap-based simultaneous confidence bands that are not only asymptotically valid but also have good small sample performance. They can, however, not easily be adjusted to be applicable to time series with missing data, as the autoregressive sieve bootstrap requires imputation of the missing values through for instance the Kalman filter. While this is certainly possible, it complicates implmentation. One can also argue about how accurate imputation methods are when a majority of the data are missing, as we face in our climatological application. Instead, we provide a much simpler alternative that requires no adjustments at all in the presence of missing data.
To illustrate our methodology, we study a time series of atmospheric ethane emissions for which almost 70% of the data points are missing. When weather conditions are unfavorable -- in particular due to cloud cover -- measurements cannot be taken. The series has previously been investigated by Franco. Atmospheric ethane is an indirect greenhouse gas which can be used as an indicator of atmospheric pollution and transport. It is emitted during shale gas extraction and since shale gas has become more and more important as a source of natural gas, nonparametric trend analysis in ethane data provides geophysicists and climatologists with a tool to link trend changes to shale gas extraction activities, as well as study long-term climatological change.
The paper is organized as follows. In Section (ref), our trend model is introduced along with the missing data generating mechanism. Section (ref) describes the estimation procedure and the construction of bootstrap confidence bands. Subsequently, Section (ref) derives the asymptotic properties of our method. Finite sample performance is analyzed in Section (ref) in a simulation study. Trends in atmospheric ethane are studied in Section (ref). Section (ref) concludes. All technical details including proofs are given in Appendix A, while Supplementary Appendices B to D provide further results.
Finally, a word on notation. We denote by $\xrightarrow{d}$ weak convergence and by $\xrightarrow{p}$ convergence in probability. Whenever a quantity has a subscript $^{\ast}$, it denotes a bootstrap quantity, conditional on the original sample. For instance, bootstrap weak convergence in probability is denoted by $\xrightarrow{d^*}_p$ GineZinn. $\left\lfloor x\right\rfloor$ stands for the largest integer smaller than or equal to $x$. For any functions $f(x)$ and $g(x)$, defined on the same domain, $f^{(i)}(x) = \frac{d^i}{d x^i} f(x)$ and $\left[f g \right]^{(i)} (x) = \frac{d^i}{d x^i} f(x) g (x)$.
Consider the following data generating process (DGP):
where $m\left(\cdot \right)$ is a smooth deterministic trend function and $z_t=\sigma_t u_t$ is a weakly dependent stochastic component. $\sigma_t$ captures unconditional heteroskedasticity and $\left\{u_t\right\}$ is a linear process
with autocovariance function $R_U(k)= \operatorname{\mathbb{E}} u_t u_{t+k}$ and long-run variance
Not all observations $y_1, \ldots, y_n$ are observed in practice. For this purpose, define the process $\{D_t\}$ as an indicator for whether the observations at each time are observed:
Assumptions (ref) to (ref) contain the formal conditions that $\{y_t\}$ and $\{D_t\}$ satisfy.
Assumption (ref) postulates that the trend $m(\cdot)$ is sufficiently smooth, which is the fundamental assumption for the estimation method to work. While it rules out abrupt structural breaks, this does not appear to be particularly restrictive for climatological applications, as many climatological processes tend to be such that change occurs gradually. In particular, as many series are measured daily or even multiple times a day, only instantaneous breaks, which are extremely unlikely in atmospheric processes, would not be covered by the smooth trend model.
Assumption (ref) allows for a wide array of unconditional heteroskedasticity. While excluding abrupt breaks, these can be allowed for by generalizing the function $\sigma(\cdot)$ to be piecewise Lipschitz as in SU. However, given the limited relevance of abrupt breaks for our climatological focus, we do not pursue this in the current paper.
Assumption (ref) is a standard linear process assumption that ensures that sufficient moments of $\left\{u_t\right\}$ exist and $\{u_t\}$ is weakly dependent and strictly stationary. These assumptions are satisfied by a large class of processes including, but not limited to, all finite order stationary ARMA models. The assumption also implies that $\Omega_U = \sigma_{\varepsilon}^2 \sum_{i=-\infty}^\infty \sum_{j=0}^\infty \psi_j \psi_{j + \left\lverti\right\rvert} < \infty$ (cf. Lemma (ref)). While our current assumption does not allow for conditional heteroskedasticity, this could be relaxed at the expense of increasing the complexity of the theoretical arguments, by allowing $\epsilon_t$ to be a martingale difference sequence. Similarly, alternative dependence concepts such as mixing, which is considered in the same bootstrap context by SU, could be used as well. However, as conditional heteroskedasticity is not the focus of our paper, we do not consider these extensions.
Assumption (ref) allows the missing data generating mechanism to be weakly dependent and non-stationary. The mixingale assumption (i) along with the summability condition on $\zeta_i$ assures weak dependence and summable autocovariances (Lemma (ref)). For technical reasons, we need to put the mixingale assumption on the product $D_s D_t$, but it directly implies that $D_t$ is a mixingale, as well, by setting $s=t$. By (ii) and (iii), the first two moments of the missing data process are allowed to vary smoothly over time, which thereby allows for instance for smooth periodical changing probabilities (e.g. due to seasonal variation), or long-term changes related to climate change. Smoothness is required as our estimator performs an implicit nonparametric estimate of the missing probability, and thus must behave similarly smooth as the trend function. Assumptions (i)-(iii) are met by a large class of generating processes, including many Markov chains with smoothly varying transition probabilities.
Assumption (iv) can be interpreted as an exogeneity assumption on the missing data generating mechanism, which for instance is satisfied if $\{D_t\}$ is independent of $\{u_t\}$. While this assumption could be argued to be restrictive, it does not appear to be problematic for our focus. Though inconclusive, some recent research has found evidence of a relation between greenhouse gases and the occurrence of cloud cover through climate change, see e.g. Norris. As cloud cover may cause missing observations, our assumption might appear restrictive, but this kind of long-run dependence can be accommodated through the slowly varying trend affecting both $\{y_t\}$ and $\{D_t\}$. As such, the exogeneity assumption mostly rules out short-run effects of ethane on cloud cover and vice versa, which we argue is reasonable.
Our goal is to conduct inference on the trend function $m(\cdot)$ defined in Section (ref). We first describe point estimation of $m(\cdot)$, followed by our bootstrap method, and finally treat the construction of the confidence bands.
We consider local polynomial estimation which is common in the nonparametric regression literature. In particular, we focus on the local constant or Nadaraya-Watson estimator Nadaraya, Watson, defined as
where $K(\cdot)$ is a kernel function and $h>0$ is a bandwidth, which should satisfy Assumptions (ref) and (ref) given below. Note that by construction of the $\{D_t\}$ series, the formulation in (ref) implies that the estimator only depends on the actually observed data.
Assumption (ref) is a standard assumption in the nonparametric kernel smoother literature, and is satisfied by many commonly used kernels. Assumption (ref) provides the range in convergence rates allowed for $h$ to ensure consistency and asymptotic normality of the kernel estimator.
The bandwidth, or smoothing parameter $h$ plays an important role. Large bandwidths produce a very smooth estimate, while small bandwidth produce a rough, wiggly trend estimate. Although Assumption (ref) gives some guidance, it does not provide us with a practical bandwidth choice. Data-driven bandwidth selection is therefore important for implementation. Leave-one-out cross-validation is the most popular data-based method for bandwidth selection, but it is designed for independent observations and therefore inappropriate for time series data. ChuMarron show that in the presence of positive correlation, this criterion systematically selects very small bandwidths, producing estimates which are too wiggly. With negative correlation, bandwidths will be large and the estimate too smooth. Therefore, ChuMarron propose to use a time series version of this criterion, called modified cross-validation (MCV). It is based on minimizing the criterion function $\frac{1}{n}\sum_{t=1}^n D_t \left(\hat{m}_{k,h}\left(\frac{t}{n}\right)-y_t\right)^2$ with respect to $h$, where
is a leave-$(2k+1)$-out version of the leave-one-out estimator of ordinary cross-validation, which leaves out the observation receiving the highest weight. Next to formal selection methods, visual inspection of the estimated trend function for a range of different bandwidths can help determining an appropriate bandwidth.
To construct confidence bands around the trend estimate, we modify the wild bootstrap, originally designed to handle heteroskedastic data DF, to account fo serial dependence. The wild bootstrap generates bootstrap errors as $ z_t^* = \xi_t^* \hat{z}_t$, where $\hat{z}_t$ are residuals of the nonparametric trend regression. In the standard wild bootstrap, the random variables $\left\{\xi_t^*\right\}$ are i.i.d. and thus, any dependence present in the data gets removed in the bootstrap errors. To overcome this drawback, Shao proposed the dependent wild bootstrap (DWB) in which $\left\{\xi_t^* \right\}$ are generated as $\ell$-dependent random variables with $\operatorname{\mathbb{C}ov}(\xi_s^*, \xi_t^*) = K_{DWB} \left(\frac{s - t}{\ell} \right)$, where $K_{DWB}(\cdot)$ is a kernel function. As the tuning parameter $\ell$, it has to be selected by the user.
Building on this idea, SU propose the autoregressive wild bootstrap (AWB) where $\left\{\xi_t^* \right\}$ is generated as an AR(1) process with parameter $\gamma = \gamma(n)$. The AWB has as advantage over the DWB that it is easier to implement and has a more intuitive interpretation. Moreover, as $\left\{\xi_t^* \right\}$ is not $\ell$-dependent, the AWB has the potential to capture more serial correlation and to be less sensitive to the choice of tuning parameter $\gamma$. In the context of unit root testing, SU show that the AWB generally has a superior finite sample performance compared to the DWB. For these reasons, we mainly focus on the AWB in the following, although we consider the DWB in our simulation study as well. The AWB algorithm can be described as follows.
Note that in Step 3 we only have to draw bootstrap observation for the dates where $D_t = 1$ and we observed the realization $y_t$; in the description the missing ones are artificially set to zero, but they are actually not used anywhere. In Step 2, we generate $\left\{\xi_t^* \right\}$ for all $t=1,\ldots,n$, although subsequently we only use the subset that corresponds to the actually observed data points. The missing data structure is preserved in the bootstrap sample, while the correlation between consecutive non-missing observations is determined only by their distance, which ensures a coherent bootstrap sample. In this way, the missing data structure is automatically taken into account in the bootstrap without any need for modifications.
Although we suggest to generate $\{\nu_t^*\}$ as a sequence of normally distributed random variables, inspection of the proofs shows normality is not needed; all one needs is a sequence of i.i.d. random variables with $\operatorname{\mathbb{E}}^* \nu_t^* = 0$, $\operatorname{\mathbb{E}}^* \nu_t^{*2} = 1 - \gamma^2$ and $\operatorname{\mathbb{E}}^{*4} \nu_t^{*4} < \infty$. Normality is simply a convenient option that is easy to implement; alternatively one could implement a variant of the Rademacher distribution (corrected for the right variance), which for the independent wild bootstrap has good properties DF.
For the tuning parameter $\gamma$, we follow SU and let $\gamma=\theta^{1/\ell}$ where $\ell$ is the “block length” parameter also found in the DWB and $0 < \theta < 1$ is a fixed parameter. This specification has the advantage that $\ell$ can be interpreted in a similar way as the block length parameter in a block bootstrap; its choice constitutes a trade-off between capturing more of the dependence with a large value of the tuning parameter, and allowing for more variation in the bootstrap samples with a smaller value for $\ell$. Additionally, it provides a convenient framework for studying the theoretical properties of our method. Specifically, we need $\ell \rightarrow \infty$ as $n \rightarrow \infty$, such that $\gamma \rightarrow 1$. This is analogous to the block (and dependent wild) bootstrap, where the block size must increase to capture more dependence when the sample size increases. Assumption (ref) postulates the formal conditions that $\ell$ needs to satisfy. They imply that $\gamma \rightarrow 1$, but not too fast.
Note that we propose to use a different bandwidth $\tilde{h}$ in Step 1 of the algorithm. This is a common feature in the literature on bootstrap methods for nonparametric regression. By either selecting a larger (oversmoothing) or smaller bandwidth (undersmoothing) than used for the estimator, one can account for the asymptotic bias that is present in the local polynomial estimation, see HaHo for an extensive literature review. While undersmoothing, such as used in the related paper by NP, aims at making the bias asymptotically negligible, oversmoothing aims at producing a consistent estimator of the (non-negligible) bias. Both have advantages and disadvantages, see the extensive discussion in HaHo. We follow Buhlmann and consider a solution based on oversmoothing, which we find to work well in practice; also see HM. After presenting our theoretical results in Section (ref), Remark (ref) provides an intuition of why oversmoothing allows to consistently estimate the asymptotic bias.\footnote{HaHo propose an alternative bootstrap approach that requires neither under- nor oversmoothing, however their approach only delivers pointwise intervals, and is therefore not considered in this paper.} We now state the formal conditions that $\tilde{h}$ must satisfy in Assumption (ref); one is that $h / \tilde{h} \rightarrow 0$ as $n \rightarrow \infty$, which ensures the oversmoothing.
Pointwise bootstrap confidence intervals with a confidence level of $\left(1 - \alpha \right)$ for $m(\tau)$, are denoted by $I_{n,\alpha}^{(p)} (\tau)$ and constructed with the objective that
Using our bootstrap algorithm, we can construct such pointwise intervals as
As these intervals are constructed separately for each $\tau$, links over time cannot be established with these intervals. Therefore, we next consider how to construct simultaneous confidence bands. Let $I_{n,\alpha}^G (\tau)$, for $\tau \in G$, denote a confidence band that is simultaneous over the set $G$. Formally, we seek to construct $I_{n,\alpha}^G (\tau)$ such that
Our practical implementation follows the three-step procedure proposed by Buhlmann:
In the second step, a pointwise error $\alpha_s$ is found for which a fraction of approximately $\left(1-\alpha\right)$ of all centered bootstrap estimates falls within the resulting confidence intervals, for all points of the set $G$. As such, the confidence intervals with pointwise coverage $\left(1-\alpha_s\right)$ become simultaneous confidence bands with coverage $\left(1-\alpha\right)$.\footnote{We provide R code to implement the estimator and bootstrap confidence bands on \href{http://www.stephansmeekes.nl}{www.stephansmeekes.nl}.}
We first provide the pointwise limiting normal distribution of the local constant estimator $\hat{m}(\cdot)$. Although the result is similar to the non-bootstrap part of Theorem 3.1 in Buhlmann, we extend the asymptotic theory for the local constant estimator to allow for the presence of nonstationary volatility and missing data. As we feel this is a noteworthy result in its own right, we present it in Theorem (ref).
The term $B_{as}(\tau)$ reflects the familiar asymptotic bias generally found in local polynomial estimators, although the exact form is different due the presence of the missing data parameter $p(\tau)$. The asymptotic variance $\sigma^2_{as}(\tau)$ is not only affected by $p(\tau)$, but also by the volatility process $\sigma^2 (\tau)$. If one were to use these distributions directly for inference, one would need to plug in consistent estimators of these nuisance parameters. However, as we show next, in the bootstrap these are automatically consistently estimated, and we have consistency of the autoregressive wild bootstrap method for the local constant estimator.
The pointwise validity of the bootstrap confidence intervals in the sense of (ref) follows directly from this pointwise convergence result. Note that, as the bias term $B_{as}(\tau)$ is the same in both theorems, it is consistently estimated by the bootstrap. As such, we do not need the bias to disappear, which happens when undersmoothing if $n h^{5} \rightarrow 0$, or to be $O(1)$, when $n h^{5} \rightarrow c$. Even if $n h \rightarrow \infty$, and the asymptotic bias dominates the stochastic variation, the bootstrap correctly mimics this and can be used for asymptotically valid inference. As such, we can relax the assumption in Buhlmann that $h \sim C n^{-1/5}$ to allow for a wider range of bandwidths. In practice, this means that the bootstrap provides additional protection against a misspecified bandwidth, by letting the widths of confidence bands automatically adapt.
Next, to study the validity of simultaneous confidence bands as in (ref), we consider $h$-neighborhoods around time points $\tau$. We do so because estimates $\hat{m}(\tau_1)$ and $\hat{m}(\tau_2)$ are asymptotically independent for $\tau_1 \neq \tau_2$ being two fixed distinct time points. When the distance between $\tau_1$ and $\tau_2$ is of order $h$, the estimators show a non-zero correlation. Therefore a major benefit of “zooming in” on local $h$-neighborhoods is that we can study how the bootstrap mimics the correlation between close points, a feature which is lost when considering simultaneity globally.
Theorem (ref) establishes the uniform validity of the bootstrap within an $h$-neighborhood around any point $0<\tau_0<1$, where, since $h=o(1)$, we assume without loss of generality that $m(\tau_0+\tau h)$ is always defined. Note that the interval $[-1,1]$ is mainly chosen out of convenience, and the results can trivially be shown to hold over any interval $[\tau_0 - a h, \tau_0 + b h]$ with $0<a,b < \infty$. Moreover, it follows directly from Theorem (ref) that the bootstrap will be valid uniformly on sets that contain a union of any finite number of such $h$-neighborhoods, see e.g. Buhlmann. While, in finite samples, one can always take $h$ and the intervals such that the full sample is covered in $G$, this kind of “too large” simultaneity should be considered with caution, as this is not what the asymptotic analysis covers. Although simultaneity over such local sets might appear less attractive than simultaneity over the whole sample, it can nevertheless be of great interest in applications. For example, constructing confidence bands with simultaneous coverage over two time periods - one located early in the sample and the other one at the end - is useful when judging if there was an upward (or downward) movement of the trend at the end of the time period when compared to the beginning. This allows the empirical researcher to draw conclusions about developments spanning time stretches, which is not possible with pointwise confidence intervals.
For the simulation exercise, we simulate time series with a trending behavior in both mean and variance, inspired by patterns observed in climatological time series, and allow for similar patterns of missing data. We will first describe the setting and then present and discuss the results.
We consider the following smooth transition model:
where for $\lambda>0$,
The error term $\left\{u_t\right\}$ follows an ARMA$(1,1)$ model
where we vary the parameters $\phi$ and $\psi$ to investigate the impact of serial correlation on our method. The variance of $\epsilon_t$ is normalized such that the signal to noise ratio does not depend on the specific choice of the AR and MA parameter. Furthermore, we introduce heteroskedasticity with the process $\left\{\sigma_t\right\}$. We consider two scenarios, where $\sigma_t$ is constant over time or $\sigma_t = \sigma(t/n)$, with the volatility process $\sigma(\tau)$ given by
Equation (ref) is a shifting mean model as considered by GT, and can be seen as a smooth transition version of a broken trend model with one break. The function $G(\tau,\lambda,c)$ as given in (ref) is the transition function with time as transition variable. Its inputs apart from time are the location of the shift -- the parameter $c$ -- as well as the smoothness of the shift, determined by $\lambda$. For large values of $\lambda$ the shift happens almost instantaneous, while it is smoother for smaller values of this parameter. In our simulations, we fix $\lambda=10$. The other parameters of our DGP will be chosen in such a way that the time series experiences a downward trend during the first three quarters which turns into a steeper upward trend in the last quarter.
This mimics the general pattern which is expected to occur in atmospheric ethane time series and therefore fits our application well. More specifically, this means we set the location of the shift to occur at $c=0.9$. The slope of the trend gradually changes from $\beta_1=-1$ before the shift to $\beta_2=2.5$ after the shift. This is illustrated in Figure (ref). For the variance process, inspired by the series considered in the empirical application, we consider a cyclical component with trend. We have to choose four parameters in (ref): the start and end point of the trend -- $\sigma_0$ and $\sigma_{\ast}$ -- as well as the specifics of the cyclical component. The parameter $a$ fixes the amplitude of the cycle, while $k$ determines how many cycles there are. We set $\sigma_0=1$, $\sigma_{\ast}=2$ and consider different combinations of values for $a$ and $k$. We let $a=0.3,0.5,0.7$ and $k=2,3,4$. An example of this process is displayed in Figure (ref).
We also consider different degrees of dependence by varying the AR and MA parameters. For the AR parameter we take $\phi=0,0.2,0.5,-0.5$, while the MA parameter varies between $\psi=0$, $\psi=0.2$ and $\psi=0.5$. We only look at pure AR or MA processes with these coefficient values. The different specifications will be abbreviated in the tables with self-explanatory names, e.g. we write $AR_{-0.5}$ for $\phi=-0.5$, $\psi=0$ and use $MA_{0.5}$ when $\phi=0$ and $\psi=0.5$.
In addition, we consider cases of missing data for which we generate a missing pattern that is representative for the ethane data. We implement a first-order Markov Chain for $D_t$ with transition probabilities
which are estimated from the ethane time series considered in Section (ref). This transition matrix results in an average fraction of around 70% missing observations.
In the estimation step, we apply the local constant estimator based on the Epanechnikov kernel which is given by the function $K(x)=\frac{3}{4}(1-x^2)\mathbbm{1}_{\left\{|x|\leq 1\right\}}$. For the bandwidth parameter $h$ we use $h=0.02$, $h=0.04$ and $h=0.06$. In the first step of the bootstrap procedure we follow the recommendation of Buhlmann to use $\tilde{h}=Ch^{5/9}$ with $C=2$. In the second step of the bootstrap, we consider different values for the AR parameter $\gamma$; next to $\gamma=0$, which reduces the AWB to a standard wild bootstrap (WB), we also consider $\gamma=0.2,0.4,0.6$.
For each specification, we run 5000 Monte Carlo simulations. We report average pointwise as well as simultaneous coverage for a sample size of $n=200$, based on $B=999$ bootstrap replications. For ease of comparison, we choose the sample size in cases with missing data such that we have approximately 200 data points remaining. Given the large fraction of missings, an expected effective sample size of 200 translates into a original sample size for our Markov Chain of $n=666$.
The nominal coverage in all cases is 95%. We also report the average median length of the confidence intervals in parenthesis underneath the respective coverage. For simultaneous coverage, the trend curve has to lie within the confidence bands for all points of the considered set $G$, for which we take the two sets $G_{sub}$ and $G$ considered by Buhlmann, where $G_{sub}=U_1(h)\cup U_4(h)$ and $G=\bigcup_{i=1}^4U_i(h)$, with $U_i(h)=\left\{(i/5)-h+j/100;\; j=0,...,[200h]\right\}$.
To compare the performance of our AWB method to related bootstrap methods, we also implement the dependent wild bootstrap (DWB) and the sieve wild bootstrap (SWB), with a standard normal distribution for generation of the wild bootstrap errors. Since the SWB cannot easily be adapted to work with missing data, we provide results for this method only for the cases with no missing data. For the DWB, we convert the $\gamma$ parameter into the corresponding value for the tuning parameter $\ell$, using the formula $\gamma=\theta^{1/\ell}$. SU found in their simulation study that $\theta=0.01$ provides a sensible conversion between the AWB and DWB, in the sense of yielding comparable performance of the two methods, therefore we use $\theta=0.01$ as well to convert the AWB parameter into the DWB parameter. This tuning parameter does not exist in the case of the SWB; instead, the lag length has to be selected, for which we use AIC.
First, we report results for equally spaced data with no missing observations and a variance process with the same specifications ($k=4$ and $a=0.5$) as displayed in Figure (ref). The results on pointwise coverage are given in Table (ref), while Tables (ref) and (ref) show simultaneous coverage probabilities for the two sets $G.sub$ and $G$, respectively. The tables consist of three main blocks, one for each bandwidth. Within each block, the individual rows contain results for different combinations of AR and MA parameters. Results for other choices of the variance parameters $k$ and $a$, including the homoskedastic case, are available in Supplementary Appendix C. Qualitatively, these settings yield the same conclusions as the ones considered here.
Table (ref) shows that the autoregressive wild bootstrap provides confidence intervals with good pointwise coverage in the presence of heteroskedasticity and mild autocorrelation. For the independent case and the cases with negative or small positive correlation, the coverage probabilities are close to the nominal level. The only specifications for which the coverage lies below the nominal level are when $\phi=0.5$ and $\psi=0.5$. In these cases, the data deviate from the trend line in clusters due to the strong positive correlation. This causes the nonparametric estimate to go through these clusters and thus, to deviate significantly from the true trend. The confidence bands are in these situations not wide enough to cover the true trend, resulting in too low coverage. Interestingly, all methods/tuning parameters are similarly affected.
Concerning the autoregressive parameter of the wild bootstrap, we can observe that whenever the data are serially correlated, the autoregressive wild bootstrap ($\gamma\neq 0$) provides better coverage than the standard wild bootstrap ($\gamma=0$). In addition, with stronger correlation, a larger value for $\gamma$ should be preferred, except that the case $\gamma=0.6$ provides consistently lower coverage, indicating that simply going for a very large value of $\gamma$ is not sensible in practice. However, even if we do see these patterns, in general, the coverage probabilities do not vary substantially with the autoregressive parameter and therefore, the sensitivity to this parameter appears to be fairly limited.
When we look at the different blocks of Table (ref), the bootstrap shows a similar overall performance regardless of the value we select for the bandwidth parameter. Since the bandwidth plays such an important role in nonparametric estimation, yet there are no fully satisfactory ways to select optimal bandwidth from the data in most applications, robustness to bandwidth “misspecifcation” is an important finding. This implies that the bootstrap can correct for poorly chosen bandwidths.
We observe similar patterns in Tables (ref) and (ref), while overall coverage is lower for the set $G$ than for $G_{sub}$. This is not surprising, since the former set covers twice as many points as the latter. Interestingly, the confidence bands are consistently more narrow with $G$ than they are with $G_{sub}$. This appears counterintuitive at first, as $G$ is twice as large as $G_{sub}$. However, while $G_{sub}$ is made up entirely of points relatively close to the boundaries, for which estimation is more variable, $G$ additionally contains “stable” regions closer to the center. It may be that the stability of these regions has an offsetting effect compared with the boundary regions, reducing the size of the intervals. As a side effect, overall coverage is also reduced. As such, if one is only interested in coverage near the beginning and the end of the sample, it may be wiser to only attempt to achieve uninformity over these regions, rather than over the full sample.
Similar to the pointwise coverage results, the simultaneous coverage is close to the nominal level for the two cases with negative correlation as well as the independent case. Weak positive correlation can also be handled decently. The cases $\phi=0.5$ and $\psi=0.5$ are more problematic, as coverage drops to around 60% for $G$ and 70% for $G_{sub}$. In these cases, the smallest bandwidth $h=0.02$ seems to be preferred.
Comparing the AWB to the other two bootstrap methods, we can see that in almost all cases, AWB and DWB show similar results and they outperform the sieve version. Often, the AWB results in slightly higher coverage with shorter intervals. Only when $\gamma=0.6$, the DWB displays better coverage. Further increases of this parameter did not lead to improvements. An exception is the DGP with negative correlation, where the DWB displays coverage that is too high, independent of the choice of $\gamma$. In such cases, the AWB is often more accurate for $\gamma=0.4$. In all other cases, we see that for both methods, the best performance is similar in magnitude but obtained at a different value of the tuning parameter. Since the tuning parameters do not have exactly the same meaning in both methods, there is no reason to expect identical variation. We chose $\theta=0.01$ to link the two methods; changing this value will likely change the relation between the methods as well.
Next, we consider the setting with missing data. Given the previous findings, we restrict ourselves to one bandwidth ($h=0.06$) but consider all AR and MA models. The results for pointwise as well as simultaneous coverage probabilities are given in Table (ref); further results, with similar conclusions, are available in Supplementary Appendix C.
The first block presents pointwise coverage, while the second and third blocks show results for the sets $G_{sub}$ and $G$. The AWB performs well even if a significant proportion of the data are missing, as both pointwise and simultaneous coverage is close to the nominal level for almost all models. There is a significant increase in pointwise coverage for the cases with strong positive correlation, which is now close to 90%. The same increase is visible for both $G_{sub}$ and $G$. This phenomenon does not come as a surprise, as the missing data points create space between consecutive observations, thus effectively reducing the serial dependence between observed points. Comparing the AWB with the DWB, we can see that the coverage is slightly closer to 95% for the AWB in many cases. As before, the best performance is obtained at different values of $\gamma$ for the AWB and the DWB. The general pattern, however, is as in the previous setting. The DWB outperforms the AWB for higher values of $\gamma$, while the AWB obtains the most accurate coverage for smaller values of this parameter. A notable exception are again the cases with negative autocorrelation. Results for other bandwidths are similar, and presented in Supplementary Appendix C.
Overall, this simulation study indicates that the autoregressive wild bootstrap performs well in most of our considered scenarios. The pointwise confidence intervals show coverage close to the nominal level in the presence of heteroskedasticity and serial correlation. In addition, the method still performs well in the presence of missing data. The AWB provides simultaneous confidence bands with good coverage as long as the correlation does not become too strong. It outperforms the sieve wild bootstrap whenever we could compare results. In comparison to the dependent wild bootstrap, we saw that both methods provide very similar coverage probabilities, while the DWB produces slightly wider intervals.
We use our methodology to investigate the trending behavior of a time series of atmospheric ethane emissions which is derived from observations performed at the Jungfraujoch station in the Swiss Alps. This station can be found on the saddle between the Jungfrau and the M\"onch, located at 46.55$^{\circ}$ N, 7.98$^{\circ}$ E, 3580 m altitude. Ethane is the most abundant hydrocarbon gas in the atmosphere after methane and it is used as a measure of atmospheric pollution. It contributes to the formation of ground-level ozone and it influences the lifetime of methane which classifies it as an indirect greenhouse gas. This series, which has been studied by Franco, is available from the Network for the Detection of Atmospheric Composition Change website at \href{ftp://ftp.cpc.ncep.noaa.gov/ndacc/station/jungfrau/hdf/ftir/}{ftp://ftp.cpc.ncep.noaa.gov/ndacc/station/jungfrau/hdf/ftir/}. It is argued in Franco that the measurement conditions are very favorable at the Jungfraujoch location due to high dryness and low local pollution. Further details on the ground-based station and on how the measurements are obtained can be found in the aforementioned reference. It is a time series consisting of daily ethane columns (i.e. the number of molecules integrated between the ground and the top of the atmosphere) recorded under clear-sky conditions between September 1994 and August 2014 with a total of 2260 data points. Whenever more than one measurement is taken on one day, a daily mean is considered.
The average number of data points per year is 112.6 - giving an indication of the severity of the missing data problem present in this series. This shows that, in line with the above discussion, it is of major importance to use a bootstrap method which can replicate the missing data pattern correctly. The estimated transition probabilities of a first order Markov Chain reported in (ref) already indicated the presence of (weak) serial dependence in the missing data generating mechanism. As a further exploration of this mechanism, note that the local constant estimator implicitly provides an estimator for the smoothly varying proportion of non-missing data $p(\tau)$. We can write (ref) as
and $\hat{p}(\tau)$ can be seen as an estimator of $p(\tau)$.\footnote{Lemma (ref) establishes the consistency of this estimator.} In Figure (ref), we plot $\hat{p}(\cdot)$ as a diagnostic tool to investigate how data availability evolves over time. It fluctuates around the average proportion of 0.3, with a maximum of almost 0.4 and a minimum of 0.2. Although no overall trend appears to be present, the fluctuations are serious enough to cast doubt on the stationarity of the missing data generating mechanism; however, our method can handle this without problems.
In addition, the data exhibit strong seasonality, as ethane degrades faster in summer than it does in winter, causing the series to displays peaks during winter and troughs during summer. Franco take care of this seasonality with the help of Fourier terms by fitting the following model to ethane measurements $x_t$:
They continue their analysis with the residuals from this estimation, where $M=3$ is selected by inspecting the residual variance. To investigate the sensitivity of the choice of $M$, we perform a frequency domain analysis to give more insight about the form of the periodic pattern present in our data. Due to the missing data, we use the Lomb-Scargle periodogram, which is suitable in this situation (see Lomb76,Scargle82). Figure (ref) plots the periodogram of the Jungfraujoch series with the frequency, rescaled to years, on the horizontal axis.
The peaks around zero are the smooth long-run trend we model with our trend estimator. The present seasonality induces additional peaks at higher frequencies. There is a large peak at 1, representing an annual periodicity, which is so pronounced that it obscures peaks at other frequencies. Therefore, the right panel displays the same spectrum as the left panel, but with a smaller vertical axis such that other peaks are observed more clearly. Moreover, we can observe that there are peaks at 2 and 3. They are, however, not as clear-cut as the peak at 1 and might not contribute as much to the seasonality, yet provide further justification for the choice $M=3$. In Supplementary Appendix D, we consider the periodograms of the residuals of the regression on 1 up to 4 Fourier terms. These show that while inclusion of one term is clearly needed, including more terms does indeed further reduce periodicity, although for increasing $M$, the effect becomes less pronounced.
To corroborate our results and to be able to compare our findings to Franco, we additionally look at the Akaike and Bayesian information criteria as well as the residual variance (MSE in Franco) from the above regression for different values of $M$. The results are summarized in Table (ref). While the Akaike criterion (AIC) is indifferent between adding 4 or 5 Fourier terms, BIC selects 2. The residual variance (MSE) decreases by only small increments when more than 3 terms are included. Based on our analysis, it is not clear how many Fourier terms we should include; any value of $M$ between 1 and 4 seems reasonable. In the following, we report results for applying the trend estimation on the residuals of the regression with $M=3$ Fourier terms, as in Franco. In Supplementary Appendix D we perform the same analysis with $M$ varying between 1 and 4, with hardly any difference in the results.
We next investigate bandwidth selection. As suggested in Section (ref), we determine a possible bandwidth using modified cross-validation. In line with the discussion in ChuMarron, for our series, the ordinary leave-one-out cross-validation criterion selects a bandwidth which is too small ($h_{cv}=0.0006$). This value for the bandwidth parameter gives almost no smoothing of the data and the resulting trend curve is too wiggly. Leaving out $k=5$ observations on each side of any point, the modified criterion yields a value of $h_{mcv}=0.03$. Albeit a still small bandwidth, this value gives a much more reasonable picture of the trend estimate. The resulting nonparametric estimate as well as 95% simultaneous confidence bands are depicted in Figure (ref). The confidence bands are simultaneous over the whole sample. Although the validity has not been established, the algorithm works when we cover the whole sample and the results are easier to interpret. The bands are obtained using $B=999$ replications of the bootstrap procedure and an autoregressive parameter of $\gamma=0.5$.
As a robustness check, we also applied the trend estimation to the data without explicitly modeling the seasonality using Fourier terms. The nonparametric kernel estimator can be interpreted as a low pass filter which suppresses high frequency oscillations. A sufficiently large bandwidth should introduce enough smoothing to provide a trend curve which is not driven by the seasonal component. The bandwidth selected by MCV is too small for this effect to appear in our data. When we increase the bandwidth to $h=0.06$, the resulting trend shows the same pattern as the one in Figure (ref). This analysis shows that when strong seasonality is present in the data, the nonparametric kernel estimator can be used to filter out the seasonality if the bandwidth is large enough. In this case, however, bandwidth selection becomes a critical issue and the proposed MCV criterion should be applied with care. Further details can be found in Supplementary Appendix D.
We observe a slight downward trend of the ethane time series until around 2009, with local peaks in 1998 and 2002-2003, and an upward trend thereafter. This general development of the trend supports the findings in Franco who estimate a linear trend model with a break at the beginning of 2009. They find a negative slope of the trend line before the break and a positive slope after the break. As mentioned by Franco, the initial downward trend can be explained by a general emission reduction since the mid 1980's, of the fossil fuel sources in the Northern Hemisphere. This has also been reported by Simpson. The upward trend seems to be a more recent phenomenon. Studies attribute it to the recent growth in the exploitation of shale gas and tight oil reservoirs, taking place in North America, see e.g. Vinci and Franco2. Since previous studies have mainly used methods based on linear trends, the two local peaks have to our knowledge not yet been analyzed. They can potentially be explained by boreal forest fires which were taking place mainly in Russia during both periods. Geophysical studies have investigated these events in association with anomalies in carbon monoxide emissions Yurganov1, Yurganov2. In such fires, carbon monoxide is co-emitted with ethane, such that these events are likely explanations for the peaks we observe.
As a final step, we look at the standard deviation of the residuals. When estimating it with a nonparametric kernel smoother, we see a cyclical pattern with upward trend, similar to the process we generate in our simulations. We plot the estimated standard deviation in Figure (ref). This clearly shows that the residuals are heteroskedastic which further underlines the importance of a flexible bootstrap method.
In this paper we have proposed a dependent version of the wild bootstrap, the autoregressive wild bootstrap, to construct confidence intervals around a nonparametrically estimated trend. Consistency of the bootstrap has been established such that it can be used to construct pointwise and simultaneous confidence bands. While the pointwise intervals always show good coverage in finite samples, simulation results for the simultaneous bands indicate that strong positive autocorrelation leads to a drop in coverage whenever simultaneous confidence bands are considered. However, other bootstrap methods such as the dependent wild bootstrap and the sieve bootstrap are equally affected, and overall the autoregressive wild bootstrap performs at least on par with these other methods, and often outperforms them.
One major advantage of the proposed approach is its broad applicability as it can be used under general forms of serial dependence and heteroskedasticity. Furthermore, it can be applied without further adjustments when data points are missing. This feature of the autoregressive wild bootstrap is particularly relevant in economic and climatological applications where the problem of missing data is often encountered. In addition to simulation results, we provide a rigorous asymptotic analysis where asymptotically pervasive missing data are allowed for. While the missing data generating mechanism affects the asymptotic distribution of the estimator, our bootstrap method correctly mimics this and is therefore valid in the presence of general forms of missing data patterns.
An application to atmospheric ethane measurements from Switzerland demonstrates our methodology. An upward trend in this time series is an indication of increasing atmospheric pollution and it has been visible in the data for the last quarter. This finding is in line with previous studies in geophysics and provides further evidence that an increased activity in shale gas extraction might have caused an increase in the ethane burden measured over the Jungfraujoch. In addition, we find two local peaks in the ethane series, which can be explained by boreal forest fires. Natural limitations of linear trend estimation have prevented these peaks from being discovered in previous research. This underlines the flexibility of our approach compared to parametric methods.
An open end to our analysis is the choice of the autoregressive parameter in the autoregressive wild bootstrap. Although our simulation results suggest that a range of values for this parameter perform adequately, its selection in practice remains an open issue. Theoretical results on the choice of this parameter are not trivial; moreover, such theoretical results do not translate directly into selection methods with good properties in small samples. This issue therefore merits deeper study and is left as an exercise for future research.
We would like to thank guest editor Tommaso Proietti, two referees, Eric Beutner, David Hendry, Franz Palm and Hanno Reuvers as well as participants at the conference Econometric Models of Climate Change in Aarhus, the Maastricht Workshop on Advances in Quantitative Economics II and a seminar at Oxford University, for helpful discussions and valuable comments. We further thank Whitney Bader, Bruno Franco, Bernard Lejeune and Emmanuel Mahieu for helpful discussions regarding the ethane application. The second author thanks the Netherlands Organisation for Scientific Research (NWO) for financial support.
\numberwithin{equation}{section} \numberwithin{lemma}{section} \numberwithin{assumption}{section} {0.7\baselineskip} {0cm}