Theory of Evolutionary Spectra for Heteroskedasticity and Autocorrelation Robust Inference in Possibly Misspecified and Nonstationary Models

Theory of Evolutionary Spectra for Heteroskedasticity and Autocorrelation Robust Inference in Possibly Misspecified and Nonstationary Models

\setcounter{page}{0}

\raggedbottom

{\bf{JEL Classification}}: C12, C13, C18, C22, C32, C51\\ {\bf{Keywords}}: Fixed-$b$, HAC standard errors, HAR, Long-run variance, Nonstationarity, Misspecification, Outliers, Segmented locally stationary.

\onehalfspacing \thispagestyle{empty}

Introduction

The literature on heteroskedasticity and autocorrelation robust (HAR) inference is extensive and quite mature by now. For concreteness, consider the linear model where $x_{t}$ is a vector of regressors and $e_{t}$ is an unobservable disturbance, which can be serially correlated. It is now common practice to use OLS and correct the standard errors. This entails the estimation of the covariance matrix (referred to as the long-run variance, LRV) of $V_{t}=x_{t}e_{t}$ or ($2\pi$ times) the spectral density of $V_{t}$ at frequency zero when the latter is stationary (of course, in general, the relevant process $V_{t}$ can be generated from a more complex model; e.g., a moment condition in a GMM context). Early important contributions in econometrics are newey/west:87 (newey/west:87; newey/west:94) and andrews:91 who proposed heteroskedasticity and autocorrelation consistent (HAC) estimators with some optimal properties. This approach aims at devising good estimate of the LRV of $V_{t}$. An alternative method foregoes that aim and concentrates on having a test with a pivotal non-normal limit distribution that is obtained through an inconsistent estimate of the LRV of $V_{t}$ that keeps the bandwidth at a fixed fraction of the sample size. This is the so-called fixed-$b$ HAR inference initiated by Kiefer/vogelsang/bunzel:00 and Kiefer/vogelsang:02 (Kiefer/vogelsang:02; kiefer/vogelsang:05). The drawback of this approach is that the limit distribution changes depending on the context and critical values are to be obtained numerically on a case by case basis. The literature since then has focused on providing various refinements, mostly to have tests having exact sizes closer to the nominal level.\footnote{The fixed-$b$ or post-HAC literature is vast; see dou:18, lazarus/lewis/stock:17, Lazarus et al. lazarus/lewis/stock/watson:18, \nocite{Goncalves/vogelsang:11} dejong/davidson:00, ibragimov/muller:10, jansson:04, muller:07 (muller:07; mueller:14), phillips:05, politis:11, preinerstorfer/potscher:16, potscher/preinerstorfer:18, rho/vogelsang:2020, robinson:98, sun:13 sun:13,sun:14,sun:14a, sun/phillips/jin:08 and zhang/shao:13.}

Most of this literature relies on stationarity with exception of the consistency results in newey/west:87 and of some results in andrews:91 which, however, do not provide accurate approximations. Yet, another strand of the literature has argued convincingly that the processes governing economic data $\{x_{t}\}$ and the errors in the relevant regressions $\left\{ e_{t}\right\} $ are nonstationary.\footnote{By nonstationary we mean non-constant moments. As in the literature, we consider processes whose sum of absolute autocovariances is finite. That is, we rule out processes with unbounded second moments (e.g., unit root).} This can occur for several reasons: changes in the moments of $x_{t}$ induced by changes in the model parameters that govern the data {[}cf. perron:89, stock/watson:96 and the surveys of ng/wright:13 and giacomini/rossi:15{]}; changes in the moments of $e_{t}$ (think about the Great Moderation with the decline in variance for many macroeconomic variables or the effects of the COVID-19 pandemic); smooth changes in the distributions governing either processes that arise from transitory dynamics; and so on. All these induces nonstationarity in $\{V_{t}\}$, which then makes $\mathbb{E}(V_{t}V'_{t-k})$ depend on both $k$ and $t$. Furthermore, even if the data and primitive shocks $\{e_{t}\}$ are stationary, many HAR testing problems are such that the relevant process $\{V_{t}\}$ is stationary under the null hypothesis but is affected by changes in means (or other forms of nonstationarity) under the alternative. This occurs, for instance, when using tests involving structural breaks based on estimating the model under the null hypothesis; e.g., popular tests for forecast evaluation {[}e.g., diebold/mariano:95{]}, tests for forecast instability {[}cf. casini_CR_Test_Inst_Forecast, giacomini/rossi:09 and perron/yamamoto:18{]}, tests for structural change {[}cf. casini/perron_Oxford_Survey and perron:06{]}. When standardized by classical HAC estimators such tests may suffer from issues such as non-monotonic power, i.e., power that goes to zero as the alternative gets farther away from the null value. Various forms of misspecication and/or nonstationarity generate low frequency contamination and make the series or residuals appear much more persistent. As a consequence, HAC standard errors are too large and when used as normalizing factors of test statistics, the tests lose power {[}see casini/perron_Low_Frequency_Contam_Nonstat:2020 for formal details{]}.\footnote{A partial list of works that present evidence of power issues with HAR inference is altissimo/corradi:2003, casini/perron_Oxford_Survey (casini/perron_Oxford_Survey, casini/perron_Lap_CR_Single_Inf, casini/perron_SC_BP_Lap), chan:2020, crainiceanu/vogelsang:07, juhl/xiao:09, kim/perron:09, martins/perron:16, perron/yamamoto:18, shao/zhang:2010, vogeslang:99, xu:2013, zhang/lavitas:2018.} This applies even more forcefully to the fixed-$b$ type methods and to the recent refinements by lazarus/lewis/stock:17 and Lazarus et al. lazarus/lewis/stock/watson:18, since they involve more lagged autocovariances (or long bandwidths) and, hence, larger contaminations.

This points to the importance of extending the methods for HAR inference so that they have the correct size and good power even under nonstationarity. This is the aim of the paper. We first develop a theoretical framework under which to analyze the statistical properties of our suggested estimate. We introduce a class of nonstationary processes which possess a spectrum that varies both over frequencies and time, thereby encompassing the nonstationary models used in applied work. We work in an infill asymptotic setting akin to the one used in nonparametric regression {[}cf. robinson:89{]}. For a process $V_{t}$, its spectrum at frequency $\omega$ and time $u=t/T$, denoted by $f\left(u,\,\omega\right)$, is allowed to change slowly yet continuously as well as to change abruptly in $u$ at a finite number of time points; the latter allows for structural breaks in the spectrum of $V_{t}$. We label this class as Segmented Locally Stationary (SLS). It is related to the locally stationary processes introduced by dahlhaus:96 that have the characterizing property of behaving as a stationary process in a small neighborhood of $u$. This is achieved via smoothness of $f\left(u,\,\omega\right)$ in $u$. By allowing discontinuities across some segments, we can deal with relevant features such as structural change, regime switching-type and threshold models {[}cf. bai/perron:98, casini/perron_CR_Single_Break (casini/perron_Oxford_Survey, casini/perron_Lap_CR_Single_Inf, casini/perron_SC_BP_Lap and casini/perron_CR_Single_Break), hamilton:89 and hansen:00ecma{]}. The SLS class extends some of the analysis of dahlhaus:96 to processes having a more general time-varying spectrum.\footnote{A few authors used a notion of local stationarity that allows for breaks {[}see, e.g., dahlhaus:2009 and last/shumway:08{]}. However, none of these works was concerned with HAR inference. dahlhaus:2009 presented some results for local spectral density estimation and required smoothness (see Example 4.2 there). last/shumway:08 considered testing for change-points in a locally stationary series which under the alternative hypothesis results in a piecewise locally stationary series. Furthermore, our notion of SLS processes and related framework are more general from a theoretical standpoint since we provide precise definitions and establish theoretical results about the identification of the local spectral density when there are discontinuities. } Our framework is of independent interest and can be useful in many contexts in econometrics if one is interested in deriving the properties of estimators or inference under nonstationarity.

Under this framework, we introduce a double kernel HAC (DK-HAC) estimator in order to flexibly account for nonstationarity and we show that it is robust to low frequency contamination and other misspecifications. This entails an extension of the classical HAC estimators since in addition to the usual smoothing procedure over lagged autocovariances, it applies a second smoothing over time for each lagged autocovariance, involving a second kernel and bandwidth. If $\{V_{t}\}$ is Segmented Locally Stationary, $\mathbb{E}(V_{t}V'_{t-k})$ changes smoothly in $t$, as long as $t$ is away from the change-points in the spectrum $f\left(t/T,\,\omega\right)$. Thus, the smoothing over time yields good estimates for the time path of $\mathbb{E}(V_{t}V'_{t-k})$ for all $k$. We determine the optimal kernels and optimal values for both bandwidth sequences under a mean-squared error (MSE) criterion. We establish new MSE bounds that show how nonstationarity affects the bias-variance trade-off and are more informative than previously established MSE bounds. We use them to construct data-dependent bandwidths relying on the “plug-in” approach. Unlike andrews:91, our candidate parametric models have time-varying parameters which can be estimated by applying standard methods to local data, akin to using rolling regressions. The procedure depends on three elements: the bandwidths for the smoothing over autocovariances and over time, and a block size to separate the regimes. In this paper, we consider a sequential bandwidth selection procedure by first deriving the optimal bandwidth for smoothing over time, then conditioning on this to obtain the optimal bandwidth for smoothing over autocovariances.

The DK-HAC estimators can result in HAR tests that are oversized when there is high temporal dependence in the data, a well-known problem for all methods, though for ours these distortions are relatively minor compared to, e.g., the methods of newey/west:87 and andrews:91. Still, in order to improve the size control of HAR tests, casini/perron_PrewhitedHAC, using the theory of this paper, propose a nonparametric nonlinear VAR prewhitened DK-HAC estimators. This form of prewhitening differs from those discussed previously {[}e.g., andrews/monahan:92 and rho/shao:13{]} in that it accounts explicitly for nonstationarity. HAR tests based on prewhitened DK-HAC estimators have size control competitive to fixed-$b$ HAR tests when the latter work well (i.e., under stationarity). Notably, non-prewhitened and prewhitened DK-HAC have excellent power properties even when existing HAR tests have serious issues with power.

Comparison to Existing Literature

There are two main approaches to HAR inference differing on whether the LRV estimator is consistent or not. The classical approach relies on consistency, which results in HAC estimators {[}cf. newey/west:87 newey/west:87,newey/west:94, andrews:91 and hansen:92ecma{]}, and on bandwidths chosen via MSE criterion. Inference is standard because HAR tests follow asymptotically standard distributions. The researcher then uses corrected standard errors and asymptotic critical values. It was shown early that classical HAC standard errors can result in oversized tests when there is substantial temporal dependence. This stimulated a second approach based on inconsistent LRV estimators that keep the bandwidth at some fixed fraction of the sample size {[}cf. Kiefer/vogelsang/bunzel:00{]}. Because of the inconsistency, inference is nonstandard and HAR tests do not asymptotically follow standard distributions. Critical values are to be obtained numerically. Long bandwidths/fixed-$b$ methods require stationarity and reduce the oversize problem of HAR tests. The bandwidth choice is often based on testing-oriented criteria {[}e.g., sun/phillips/jin:08{]}. Our approach falls in the first category; we propose HAC estimators and standard HAR inference.

We now compare in detail our approach to the existing literature. We believe that a fair comparison has to consider the following four criteria: (1) applicability to general HAR inference tests; (2) size of HAR tests; (3) power of HAR tests; (4) theoretical validity under stationarity/nonstationarity. In terms of (1), it is clear that the HAC and DK-HAC estimators are generally and immediately applicable to any HAR inference test and that they are simple to use in practice. This also explains why the classical HAC estimators have become the standard practice in econometrics. Methods that rely on long bandwidths/fixed-$b$ do not share the same property. They are not generally applicable to HAR inference tests because a researcher would first need to derive a new asymptotic non-standard fixed-$b$ distribution. This can be unfeasible in non-standard testing problems {[}e.g., tests for parameter instability, etc.{]}. Turning to (2), all existing HAR inference tests are known to be oversized when there is strong serial dependence. However, fixed-$b$ HAR tests (or versions thereof) are less oversized than other tests based on the classical HAC estimators. The stronger is the temporal dependence the larger is the difference in size between the two approaches. Our prewhitened DK-HAC estimators are competitive with fixed-$b$ HAR tests in controlling the size. Moving to (3), prewhitened and non-prewhitened DK-HAC estimators have excellent power under either stationarity or nonstationarity whereas existing methods have serious problems with power under nonstationarity or under nonstationary alternative hypotheses. These problems result in non-monotonic power and little or no power in relevant circumstances especially in HAR tests outside the stable linear regression model. Fixed-$b$ or long bandwidths methods suffer most from these problems. Finally, turning to (4), our method like the classical HAC approach is valid under nonstationarity whereas methods using long bandwidths/fixed-$b$ are only valid under stationarity {[}cf. casini_fixed_b_erp{]}. It should be mentioned that the fixed-$b$ approach is shown to achieve (pointwise) higher-order refinements under stationarity while the MSE-based optimality of the HAC or DK-HAC estimators pertains only to the first-order but it holds under nonstationarity.\footnote{preinerstorfer/potscher:16 pointed out some limitations of these approaches because optimality does not hold uniformly over DGPs. They showed negative non-asymptotic results about size and power of HAR tests when one allows for all correlation structures corresponding to stationary Gaussian AR(1) processes. }

Recently, lazarus/lewis/stock:17 and Lazarus et al. lazarus/lewis/stock/watson:18 made some progress on the applicability of fixed-$b$ methods. They showed that the $t$-test using a LRV estimator based on equally-weighted cosine (EWC) under fixed-$b$ asymptotics can achieve a $t$-distribution with the degrees of freedom depending on the bandwidth choice.\footnote{hwang/sun:2017 proposed a modification to the trinity of test statistics in the two-step GMM setting and showed that the modified test statistics are asymptotically $F$ distributed under fixed-$b$ asymptotics and stationarity. However, following the logic of lazarus/lewis/stock:17, procedures such as e.g., bootstrap-based autocorrelation robust tests, modification of the test statistic, etc., are not as simple as HAC-based inference and therefore may find less traction than HAC-based inference in empirical work. } However, casini/perron_PrewhitedHAC showed that EWC is oversized relative to the original fixed-$b$ of Kiefer/vogelsang/bunzel:00 and to the prewhitened DK-HAC when there is strong dependence. lazarus/lewis/stock:17 relied on the Neyman-Pearson Lemma or simply “apple-to-apple comparison” to compare HAR tests. This is certainly a reasonable criterion. The lemma suggests to compare the power of tests that have an empirical size no greater than the significance level. However, the Neyman-Pearson Lemma alone does not suffice to find the “best” test in this context because all tests are oversized when there is strong dependence. Indeed, it does not even apply in this context. We face a trade-off between size and power. Our proposed method is competitive with the existing methods which is least oversized {[}i.e., original fixed-$b$ of Kiefer/vogelsang/bunzel:00{]} and has excellent power even when existing HAR tests do not have any. We believe that our method strikes a good balance with respect to criteria (1)-(4).

Our approach is different from methods based on subsampling of $t$-statistics {[}see, e.g., ibragimov/muller:10{]}. The latter rely on splitting the sample in subsamples and estimating the model within each subsample. Under the assumption that the estimates from the subsamples are asymptotically independent, the test statistic based on an average of estimates across subsamples follows asymptotically a $t$-distribution. In terms of point (1) above, this approach is not general enough compared to the HAC/DK-HAC approach because this changes the test statistic and its asymptotic distribution. Also, subsampling test statistics which are not $t$-tests can be challenging/unfeasible and would require at best extra work in general HAR inference contexts to derive the new distribution. Finally, our simulation experience (not reported) suggests that this method suffers from the same finite-sample issues about size and power as the fixed-$b$ methods.

Related Work

This paper is part of a set of papers on HAR inference by the author and collegues. The current paper provides the core theoretical and empirical elements that are used in all other papers, which can be viewed as providing extensions or refinements. casini/perron_PrewhitedHAC used our theoretical framework to derive minimax MSE bounds for LRV estimation that are sharper than previously established and extended some of our theoretical results to general nonstationarity. As a finite-sample refinement, they also developed a new prewhitening procedure robust to nonstationarity for DK-HAC estimators. Even though the latter procedure is included in our simulations, we established the corresponding theoretical results in a separate paper because of the extent of the work needed in the analysis. casini/perron_Low_Frequency_Contam_Nonstat:2020 showed analytically that the poor finite-sample performance of existing LRV-based HAR tests under nonstationarity and misspecification is induced by low frequency contamination. Belotti et al. belotti/casini/catania/grassi/perron_HAC_Sim_Bandws used our theoretical framework to propose alternative data-dependent bandwidths for DK-HAC estimators that are optimal under a global MSE criterion. casini:change-point-spectra considered change-point detection in time series with evolutionary spectra. Initially, it was intended to be used in the method suggested in this paper to improve the finite-sample size and power properties, given that we work with segmented locally stationary processes. However, it turns out that the current method to select the blocks (see Section (ref)) is able to handle even abrupt structural change.

The remainder of the paper is organized as follows. Section (ref) introduces the statistical setting and the new HAC estimator. Section (ref) presents consistency, rates of convergence and the asymptotic MSE results for the DK-HAC estimators. Asymptotically optimal kernels and bandwidths are derived in Section (ref). A data-dependent method for choosing the bandwidths and its asymptotic properties are discussed in Section (ref). Section (ref) presents a Monte Carlo study. Section (ref) concludes the paper. The supplemental materials {[}cf. casini_hac_supp and an additional supplement not for publication{]} contain some implementation details and all mathematical proofs. The code to implement our methods is provided in $\mathrm{\mathsf{Matlab}}$, $\mathrm{\mathrm{\mathsf{R}}}$ and $\mathrm{\mathrm{\mathsf{Stata}}}$ languages through a $\mathrm{\mathsf{Github}}$ repository.

The Statistical Environment

To motivate our approach, consider the linear regression model estimated by least-squares (LS): $y_{t}=x'_{t}\beta_{0}+e_{t}$ $(t=1,\ldots,\,T)$, where $\beta_{0}\in\Theta\subset\mathbb{R}^{p}$, $y_{t}$ is an observation on the dependent variable, $x_{t}$ is a $p$-vector of regressors and $e_{t}$ is an unobserved disturbance. The LS estimator is given by $\widehat{\beta}=(X'X)^{-1}X'Y$, where $Y=(y_{1},\ldots,\,y_{T})'$ and $X=(x_{1},\ldots,\,x_{T})'$. Classical inference about $\beta_{0}$ requires estimation of $\mathrm{Var}(\sqrt{T}(\widehat{\beta}-\beta_{0}))$ where

where “$\triangleq$” is used for definitional equivalence. Consistent estimation of $\mathrm{Var}(\sqrt{T}(\widehat{\beta}-\beta_{0}))$ relies on consistent estimation of $\mathrm{lim}_{T\rightarrow\infty}T^{-1}\sum_{s=1}^{T}\sum_{t=1}^{T}$ $\mathbb{E}(e_{s}x_{s}(e_{t}x{}_{t})')$. More generally, one needs a consistent estimate of $J\triangleq\mathrm{lim}_{T\rightarrow\infty}J_{T}$ where $J_{T}=T^{-1}\sum_{s=1}^{T}\sum_{t=1}^{T}\mathbb{E}(V_{s}(\beta_{0})$ $V_{t}(\beta_{0})')$ with $V_{t}(\beta)$ being a random $p$-vector for each $\beta\in\Theta$. For the linear regression model, $V_{t}(\beta)=(y_{t}-x'_{t}\beta)x_{t}$. For HAR tests outside the regression model $V_{t}(\beta)$ takes different forms.\footnote{If one suspects that $\beta_{0}$ may not be constant, one can use appropriate tests for parameter instability. However, our discussions and methods still apply because these tests are HAR inference tests and one needs a LRV estimate of $J$ based on the appropriate $V_{t}.$ } Hence, our problem is to estimate $J$ when $\{V_{t}\}$ is a Segmented Locally Stationary process, as defined in Section (ref).\footnote{casini/perron_PrewhitedHAC extended the results to the case where $\left\{ V_{t}\right\} $ is generally nonstationary (i.e., $\left\{ V_{t}\right\} $ is a sequence of unconditionally heteroskedastic random variables). } Such estimate can then be used to conduct HAR inference in the usual way using the theory developed below. By a change of variables, $J_{T}$ can be rewritten as

and $V_{t}=V_{t}\left(\beta_{0}\right)$. The rest of this section is structured as follows. In Section (ref) we introduce a new class of nonstationary time series that we use as the underlying framework for our theoretical analysis. Section (ref) presents the DK-HAC estimator. We adopt the following notational conventions. The $j$th element of a vector $x$ is indicated by $x^{\left(j\right)}$ while the $\left(j,\,l\right)$th element of a matrix $X$ is indicated by $X^{\left(j,\,l\right)}$. $\mathrm{tr}(\cdot)$ denotes the trace and $\otimes$ denotes the tensor product. The $p^{2}\times p^{2}$ matrix $C_{pp}$ is a commutation matrix that transforms $\mathrm{vec}\left(A\right)$ into $\mathrm{vec}\left(A'\right)$, i.e., $C_{pp}=\sum_{j=1}^{p}\sum_{l=1}^{p}\iota_{j}\iota_{l}'\otimes\iota_{l}\iota_{j}'$, where $\iota_{j}$ is the $j$th elementary $p$-vector. $\lambda_{\max}\left(A\right)$ denotes the largest eigenvalue of $A$. $W$ and $\widetilde{W}$ are used for $p^{2}\times p^{2}$ weight matrices. $\mathbb{C}$ is used for the set of complex numbers. $\overline{A}$ is used for the complex conjugate of $A\in\mathbb{C}$. Let $0=\lambda_{0}<\ldots<\lambda_{m+1}=1$. A function $G\left(\cdot,\,\cdot\right):\,\left[0,\,1\right]\times\mathbb{R}\rightarrow\mathbb{C}$ is said to be piecewise (Lipschitz) continuous with $m+1$ segments if it is (Lipschitz) continuous within each segment (e.g., it is piecewise Lipschitz continuous if for each $j=1,\ldots,\,m+1$ it satisfies $\sup_{u\neq v}\left|G\left(u,\,\omega\right)-G\left(v,\,\omega\right)\right|\leq K\left|u-v\right|$ for any $\omega\in\mathbb{R}$ with $\lambda_{j-1}<u,\,v\leq\lambda_{j}$ for some $K<\infty$). We define $G_{j}\left(u,\,\omega\right)=G\left(u,\,\omega\right)$ for $\lambda_{j-1}<u\leq\lambda_{j}$. If we say piecewise Lipschitz continuous with index $\vartheta>0$, then the above inequality is replaced by $\sup_{u\neq v}\left|G\left(u,\,\omega\right)-G\left(v,\,\omega\right)\right|\leq K\left|u-v\right|^{\vartheta}$. A function $G\left(\cdot,\,\cdot\right):\,\left[0,\,1\right]\times\mathbb{R}\rightarrow\mathbb{C}$ is said to be left-differentiable at $u_{0}$ if $\partial G\left(u_{0},\omega\right)/\partial_{-}u\triangleq\lim_{u\rightarrow u_{0}^{-}}\left(G\left(u_{0},\,\omega\right)-G\left(u,\,\omega\right)\right)/\left(u_{0}-u\right)$ exists $\forall\omega\in\mathbb{R}$.

Segmented Locally Stationary Processes

Suppose $\left\{ V_{t}\right\} _{t=1}^{T}$ is defined on an abstract probability space $\left(\Omega,\,\mathscr{F},\,\mathbb{P}\right)$, where $\Omega$ is the sample space, $\mathscr{F}$ is the $\sigma$-algebra and $\mathbb{P}$ is a probability measure. In order to introduce a framework to analyze time series models with a time-varying spectrum it is necessary to introduce an infill asymptotic setting whereby we rescale the original discrete time horizon $\left[1,\,T\right]$ by dividing each $t$ by $T.$ Letting $u=t/T$ and $T\rightarrow\infty,$ this defines a new time scale $u\in\left[0,\,1\right]$ which we interpret as saying that as $T\rightarrow\infty$ we observe more and more realizations of $V_{t}$ close to time $t$, i.e., we observe the rescaled process $V_{Tu}$ on the interval $\left[u-\varepsilon,\,u+\varepsilon\right]$, where $\varepsilon>0$ is a small number.

In order to define a general class of nonstationary processes, we shall start from processes that have a time-varying spectral representation specified by:

where $i\triangleq\sqrt{-1}$, $\mu\left(t/T\right)$ is the trend function, $A\left(t/T,\,\omega\right)$ is the transfer function and $\xi\left(\omega\right)$ is some stochastic process whose properties are specified below. Observe that this representation is similar to the spectral representation of stationary processes {[}see anderson:71, brillinger:75, hannan:70 and priestley:85{]}. We shall see that the main difference is that $A\left(t/T,\,\omega\right)$ and $\mu\left(t/T\right)$ are not constant in $t$.\footnote{In HAR inference, a minimal assumption on $V_{t}$ under the null hypothesis is that it has zero mean (i.e., $\mu\left(t/T\right)=0$ for all $t$). However, in this subsection we allow for arbitrary $\mu\left(t/T\right)$ so as to introduce a general framework, also applicable under various alternative hypotheses in both within and outside the regression model. } dahlhaus:96 used the time-varying spectral representation to define the so-called locally stationary processes which are characterized, broadly speaking, by smoothness conditions on $\mu\left(\cdot\right)$ and $A\left(\cdot,\,\cdot\right)$. Locally stationary processes have been used widely in both statistics and economics, though in the latter field they are best known as time-varying parameter processes {[}see, e.g., cai:07 and chen/hong:12{]}. The smoothness restrictions exclude many prominent models that account for time variation in the parameters. For example, structural change and regime switching-type models do not belong to this class because parameter changes occur suddenly at a particular point in time. We propose a class of nonstationarity processes which allow both continuous and discontinuous changes in the parameters. Stationarity and local stationarity are recovered as special cases.

The smoothness properties of $A$ in $u$ guarantees that $V_{t,T}$ has a piecewise locally stationary behavior. Later we will require additional smoothness properties for $A$.

If there is only a single regime (i.e., $m_{0}=0$) then $V_{t,T}$ is locally stationary {[}cf. dahlhaus:96{]}. If $\mu$ and $A^{0}$ do not depend on $t$, then $V_{t,T}$ is stationary and the spectral representation of stationary processes applies. However, $m_{0}=0$ rules out structural change and regime switching models. With $m_{0}>0$, we propose a framework where parameter variation can occur either smoothly or abruptly, both being relevant for economic data.\footnote{Some authors have used alternative notions of local stationarity that allow for discontinuities (i.e., piecewise locally stationary) and have established some results in other contexts which are not related to HAR inference {[}see, e.g., dahlhaus:2009, last/shumway:08 and zhou:2013{]}. In particular, our framework is more general because we also define (and work with) the covariance between observations belonging to different regimes whereas previous works considered only the covariance between observations belonging to the same regime thereby using smoothness which restricts the framework substantially.}

Let $\left\lfloor \cdot\right\rfloor $ denote the largest smaller integer function and let $\mathcal{T}\triangleq\{T_{1}^{0},\,\ldots,\,T_{m_{0}}^{0}\}$. We define the spectrum of $V_{t,T}$ in (ref) (for fixed $T$) as

with $A_{1,t,T}^{0}\left(\omega\right)=A_{1}\left(0,\,\omega\right)$ for $t<1$ and $A_{m_{0}+1,t,T}^{0}\left(\omega\right)=A_{m_{0}+1}\left(1,\,\omega\right)$ for $t>T$. Our definition coincides with the Wigner-Ville spectrum {[}cf. martin/flandrin:85{]} when there are no change-points (i.e., $m_{0}=0$). Below we show that $f_{j,T}\left(u,\,\omega\right)$ tends in mean-squared to $f_{j}\left(u,\,\omega\right)\triangleq\left|A_{j}\left(u,\,\omega\right)\right|^{2}$ for $T_{j-1}^{0}/T<u=t/T\leq T_{j}^{0}/T$ which is the spectrum that corresponds to the spectral representation. Therefore, we call $f_{j}\left(u,\,\omega\right)$ the time-varying spectral density matrix of the process.

Let $f\left(u,\,\omega\right)=f_{j}\left(u,\,\omega\right)$ if $Tu\in(T_{j-1}^{0},\,T_{j}^{0}]$ so as to suppress the subscript $j$ from $f$. It is well-known that, even when $m_{0}=0,$ the spectral representation (ref) is not unique {[}cf. priestley:85, Chapter 11.1{]}. A consequence of Theorem (ref) is that $\{f_{j}\left(u,\,\omega\right)=\left|A_{j}\left(u,\,\omega\right)\right|^{2},\,j=1,\ldots,\,m_{0}+1\}$ is uniquely determined from the whole triangular array $\{V_{t,T}\}$.

For $Tu\notin\mathcal{T}$ with $T_{j-1}^{0}/T<u=t/T<T_{j}^{0}/T$, only the realizations of $V_{t,T}$ in the time interval $u\in\left[u-n/T,\,u+n/T\right]$ with $n\rightarrow\infty$ contribute to $f_{j}\left(u,\,\omega\right)$. Since this interval is fully contained in a segment $j$ where $A_{j}\left(u,\,\omega\right)$ is smooth, and given that the length of this interval tends to zero, $V_{t,T}$ becomes “asymptotically stationary” on this interval. The length of the interval in which $V_{t,T}$ can be considered stationary is given by $n\ln n/T^{\vartheta}\rightarrow0$ . For $Tu\in\mathcal{T}$, the arguments are different. Suppose $Tu=T_{j}^{0}$. The spectrum $f_{j,T}\left(u,\,\omega\right)$ is defined in such a way that only observations prior to $T_{j}^{0}$ are used in order to construct an approximation to $f_{j}\left(u,\,\omega\right)$. Since the length of this interval tends to zero and $A_{j}\left(u,\,\omega\right)$ is left-Lipschitz continuous, then those observations become “asymptotically stationary” and thus provide the same information about $f_{j}\left(u,\,\omega\right)$.

Given $f\left(u,\,\omega\right),$ we can define the local covariance of $V_{t,T}$ at rescaled time $u$ with $Tu\notin\mathcal{T}$ and lag $k\in\mathbb{Z}$ as \[ c\left(u,\,k\right)\triangleq\int_{-\pi}^{\pi}e^{i\omega k}f\left(u,\,\omega\right)d\omega. \] The same definition is also used when $Tu\in\mathcal{T}$ and $k\geq0$. For $Tu\in\mathcal{T}$ and $k<0$ it is defined as $c\left(u,\,k\right)\triangleq\int_{-\pi}^{\pi}e^{i\omega k}A\left(u,\,\omega\right)A\left(u-k/T,\,-\omega\right)d\omega$.

DK-HAC Estimation

In model (ref), if $m_{0}=0$ and $A^{0}$ is constant in its first argument, then $\left\{ V_{t,T}\right\} $ is second-order stationary. Its spectral density matrix is then equal to $f\left(\omega\right)\triangleq\left(2\pi\right)^{-1}\sum_{k=-\infty}^{\infty}\Gamma\left(k\right)e^{-i\omega k}$ where $\Gamma\left(k\right)\triangleq\mathbb{E}(V_{t,T}V'_{t-k,T})$. When evaluated at frequency $\omega=0$ it plays a prominent role because $\lim_{T\rightarrow\infty}J_{T}=2\pi f\left(0\right).$ Nonstationarity implies that the spectral density is time-varying since $\mathbb{E}(V_{t}V'_{t-k})$ now depends on $k$ as well as on $t$. The SLS processes introduced above accommodate this property because they have a time-varying spectrum $f\left(u,\,\omega\right)$. Accordingly, we introduce the notation $\Gamma_{u}\left(k\right)\triangleq\mathbb{E}(V_{Tu,T}V'_{Tu-k,T})$ where $u=t/T$. We show below that $\Gamma_{u}\left(k\right)=c\left(u,\,k\right)+O\left(T^{-1}\right)$ uniformly in $1\leq j\leq m+1$, $Tu\leq T_{j}^{0}$ and $k\in\mathbb{Z}.$ Under the rescaling $u=t/T,\,u\in\left[0,\,1\right]$, the limit of $J_{T}$ for SLS processes is given by,

Using the definition of $f\left(u,\,\omega\right)$ it can be shown that $J=2\pi\int_{0}^{1}f\left(u,\,0\right)du$. dahlhaus:2009 discussed how to estimate $f\left(u,\,\omega\right)$ for the scalar case under smoothness in both arguments using the smoothed local periodogram. Our goals are to estimate $J$ using a time-domain method and to relax the smoothness assumption in $u$. This is different from Dahlhaus' work that considered local problems (i.e., estimation of $f\left(u,\,\omega\right)$ under smoothness) and not full-sample problems (i.e., estimation of $J$). The class of estimators of $J$ relies on double kernel smoothing over lags and time,

where $K_{1}\left(\cdot\right)$ is a real-valued kernel in the class $\boldsymbol{K}_{1}$ defined below, $b_{1,T}$ is a bandwidth sequence discussed below, $n_{T}\rightarrow\infty$ satisfying the conditions given below, and

with $K_{2}^{*}$ being a real-valued kernel and $b_{2,T}$ is a bandwidth sequence discussed below. $\widehat{c}_{T}\left(u,\,k\right)$ is an estimate of the local autocovariance $c\left(u,\,k\right)$ of lag $k$ at time $u=rn_{T}/T$. Estimation of $c\left(u,\,k\right)$ for locally stationary processes was considered by dahlhaus:12. For positive semi-definiteness, it is necessary that $K_{2}^{*}$ takes the following form: \[ K_{2}^{*}\left(\frac{\left(r+1\right)n_{T}-\left(s-k/2\right)}{Tb_{2,T}}\right)=\left(K_{2}\left(\frac{\left(r+1\right)n_{T}-s}{Tb_{2,T}}\right)K_{2}\left(\frac{\left(r+1\right)n_{T}-\left(s-k\right)}{Tb_{2,T}}\right)\right)^{1/2}\quad\mathrm{for}\,k\geq0, \] \[ K_{2}^{*}\left(\frac{\left(r+1\right)n_{T}-\left(s+k/2\right)}{Tb_{2,T}}\right)=\left(K_{2}\left(\frac{\left(r+1\right)n_{T}-s}{Tb_{2,T}}\right)K_{2}\left(\frac{\left(r+1\right)n_{T}-\left(s+k\right)}{Tb_{2,T}}\right)\right)^{1/2}\quad\mathrm{for}\,k<0. \] Setting $K_{2}\left(x\right)=(\int_{0}^{1}h\left(x\right)^{2}dx)^{-1}h\left(x+1/2\right)^{2}$ and $N_{T}=Tb_{2,T}$, we see that positive semi-definiteness requires the use of a data taper $h\left(\cdot\right)$ with length $N_{T}$. This follows because we need each $\widehat{V}_{t}$ $\left(t=1,\ldots,\,T\right)$ to be assigned the same weight across different $k$ for any given $r$. Then, letting $\widehat{V}_{t}^{\ifmmode\begingroup\defbold{bold} \text{\ifx\math@versionbold\bfseries\fi\textdegree}\endgroup\else\textdegree\fi}=\ensuremath{\left(K_{2}\left(\left(\left(r+1\right)n_{T}-t\right)/Tb_{2,T}\right)\right)^{1/2}\widehat{V}_{t}}$ we can use the same arguments as in andrews:91 applied now to $\widehat{V}_{t}^{\ifmmode\begingroup\defbold{bold} \text{\ifx\math@versionbold\bfseries\fi\textdegree}\endgroup\else\textdegree\fi}$ to show that $J_{T}$ is positive semi-definite for the appropriate choice of $K_{1}$.

The estimator $\widehat{J}_{T}$ involves two kernels: $K_{1}$ smooths the lagged sample autocovariances, akin to the classical HAC estimators, while $K_{2}$ applies smoothing over time. The factor $T/\left(T-p\right)$ is an optional small-sample degrees of freedom adjustment. In Section (ref)-(ref), we consider estimators $\widehat{J}_{T}$ for which $b_{1,T}$ and $b_{2,T}$ are given sequences. In Section (ref), we consider adaptive estimators $\widehat{J}_{T}$ for which $b_{1,T}$ and $b_{2,T}$ are data-dependent. Observe that the optimal $b_{2,T}$ actually depends on the properties of $\left\{ V_{t,T}\right\} $ in any given block. Since the order of $b_{2,T}\left(\cdot\right)$ is the same across blocks, we omit this notation for the developments of the asymptotic results. However, when we determine the data-dependent estimate of $b_{2,T}\left(\cdot\right)$, we will estimate $b_{2,T}\left(rn_{T}/T\right)$ for each $r$. We consider the following class of kernels {[}cf. andrews:91{]},

Examples of kernels in $\boldsymbol{K}_{1}$ include the Truncated, Bartlett, Parzen, Quadratic Spectral (QS) and Tukey-Hanning kernel. We shall show below that the QS kernel has certain optimality properties:

HAC Estimation with Predetermined Bandwidths

In Section (ref) we present some asymptotic properties of $\widehat{c}\left(\cdot,\,\cdot\right)$. We use them in Section (ref) in order to establish consistency, rate of convergence and MSE properties of predetermined bandwidths HAC estimators. Let $\widetilde{J}_{T}$ denote the pseudo-estimator identical to $\widehat{J}_{T}$ but based on $\{V_{t,T}\}=\{V_{t,T}(\beta_{0})\}$ rather than on $\{\widehat{V}_{t,T}\}=\{V_{t,T}(\widehat{\beta})\}$. We first require some smoothness of $A\left(u,\,\cdot\right)$ in $u$.

We also need to impose conditions on the temporal dependence of $V_{t}=V_{t,T}$. Let

where $\left\{ V_{\mathscr{N},t}\right\} $ is a Gaussian sequence with the same mean and covariance structure as $\left\{ V_{t}\right\} $. $\kappa_{V,t}^{\left(a,b,c,d\right)}\left(u,\,v,\,w\right)$ is the time-$t$ fourth-order cumulant of $(V_{t}^{\left(a\right)},\,V_{t+u}^{\left(b\right)},\,V_{t+v}^{\left(c\right)},$ $\,V_{t+w}^{\left(d\right)})$ while $\kappa_{\mathscr{N}}^{\left(a,b,c,d\right)}(t,\,t+u,$ $\,t+v,\,t+w)$ is the time-$t$ centered fourth moment of $V_{t}$ if $V_{t}$ were Gaussian.

If $\left\{ V_{t,T}\right\} $ is stationary then the cumulant condition of Assumption (ref)-(i) reduces to the standard one used in the time series literature {[}see also Assumption A in andrews:91{]}. We do not require fourth-order stationarity but only that the time-$t=Tu$ fourth order cumulant is locally constant in a neighborhood of $u$. One can show that $\alpha$-mixing and moment conditions imply that the cumulant condition of Assumption (ref) holds.

Estimation of the Local Covariance

Let $\widetilde{c}_{T}\left(u,\,k\right)$ denote the estimator that uses $\left\{ V_{t,T}\right\} $. We consider the following class of kernels:

Results on DK-HAC Estimation with Predetermined Bandwidths

Following parzen:57, we define $K_{1,q}\triangleq\lim_{x\downarrow0}\left(1-K_{1}\left(x\right)\right)/\left|x\right|^{q}$ for $q\in[0,\,\infty);$ $q$ increases with the smoothness of $K_{1}\left(\cdot\right)$ with the largest value being such that $K_{1,q}<\infty$. When $q$ is an even integer, $K_{1,q}=-\left(d^{q}K_{1}\left(x\right)/dx^{q}\right)|_{x=0}/q!$ and $K_{1,q}<\infty$ if and only if $K_{1}\left(x\right)$ is $q$ times differentiable at zero. We define the index of smoothness of $f\left(u,\,\omega\right)$ at $\omega=0$ by $f^{\left(q\right)}\left(u,\,0\right)\triangleq\left(2\pi\right)^{-1}\sum_{k=-\infty}^{\infty}\left|k\right|^{q}c\left(u,\,k\right)$, for $q\in[0,\,\infty)$. If $q$ is even, then $f^{\left(q\right)}\left(u,\,0\right)=\left(-1\right)^{q/2}\left(d^{q}f\left(u,\,\omega\right)/d\omega^{q}\right)|_{\omega=0}$. Further, $||f^{\left(q\right)}\left(u,\,0\right)||<\infty$ if and only if $f\left(u,\,\omega\right)$ is $q$ times differentiable at $\omega=0$. We define

If $b_{2,T}^{2}/b_{1,T}^{q}\rightarrow\nu<\infty$ replaces $b_{2,T}^{2}/b_{1,T}^{q}\rightarrow0$ in part (ii), then the asymptotic bias for the case of locally stationary processes becomes

For the general case of SLS processes the term involving $\left(\partial^{2}/\partial^{2}u\right)c\left(u,\,k\right)$ is different. The second summand on the right-hand side of (ref) cancels when $\int_{0}^{1}\left(\partial^{2}/\partial^{2}u\right)c\left(u,\,k\right)du=0$. The latter occurs when the process is stationary. dahlhaus:12 presented MSE results for a pointwise estimate of $f\left(u,\,\omega\right)$ under continuity in both components by applying smoothing over $u$ and $\omega.$ His results depends on the local behavior of $f\left(u,\,\omega\right)$ at time $u$ and frequency $\omega$ whereas in our problem the MSE results depend on properties of the full time path of $f\left(u,\,0\right)$. The theorem suggests that the optimal choice of $b_{1,T}$ hinges on the degree of nonstationary in the data, a feature that does not appear from the corresponding results in the literature. The results are derived as $n_{T}\rightarrow\infty$. It is possible and indeed easier to keep $n_{T}$ fixed, in which case the results are unchanged. However, the case with $n_{T}$ fixed can have some disadvantages when the spectrum is discontinuous because then the estimator would be often dealing with observations from different regimes, which as explained above might lead to low frequency contamination. We now move to the results concerning $\widehat{J}_{T}$.

Assumption (ref)-(i,iii) is the same as Assumption B in andrews:91. As remarked above, we interpret $\beta_{0}$ as the pseudo-true parameter $\beta^{*}$ when the model is misspecified. Part (iv) of the assumption is satisfied by most commonly used kernels. In order to obtain rate of convergence results we replace Assumption (ref) with the following assumptions.

The consistency result of $\widehat{J}_{T}$ in part (i) applies to kernels $K_{1}\left(\cdot\right)$ with unbounded support and to bandwidths $b_{1,T}$ and $b_{2,T}$ such that $1/b_{1,T}b_{2,T}$ grows at rate $o(\sqrt{T/b_{2,T}})$. Part (ii) yields the consistency of $\widehat{J}_{T}$ with $b_{1,T}$ only required to be $o\left(Tb_{2,T}\right)$. This rate is slower than the corresponding rate $o\left(T\right)$ of the classical kernel HAC estimators as shown by andrews:91 in his Theorem 1-(b). However, this property is of little practical import because optimal growth rates typically are less than $T^{1/2}$; for the QS kernel HAC estimator the optimal growth rate is $T^{1/5}$ while it is $T^{1/3}$ for the Newey-West HAC estimator. Part (ii) of the theorem presents the rate of convergence of $\widehat{J}_{T}$ which is $\sqrt{Tb_{2,T}b_{1,T}}$. In Section (ref), we compare the rate of convergence of $\widehat{J}_{T}$ with that of the classical HAC estimators when the respective optimal bandwidths are used.

Optimal Kernels, Bandwidths and Choice of $n_{T}$

In this section, we show the optimality of quadratic-type kernels under MSE criterion.\footnote{Besides andrews:91 and newey/west:87 in the context of LRV estimation, the MSE-optimality criterion was also used more recently by whilelm:2015 in a GMM context to determine the optimal bandwidth of the nonparametric estimator of the optimal weighting matrix. } For $K_{1},$ the result states that the QS kernel minimizes the asymptotic MSE for any $K_{2}\left(\cdot\right)$. Let

where $\widetilde{W}_{T}$ is some $p\times p$ positive semidefinite matrix. The optimal bandwidths $b_{1,T}^{\mathrm{opt}}$ and $b_{2,T}^{\mathrm{opt}}$ satisfy the following sequential MSE criterion:

The first inequality above has to hold as $T\rightarrow\infty$. The above criterion determines the globally optimal $b_{1,T}^{\mathrm{opt}}$ given the integrated locally optimal $b_{2,T}^{\mathrm{opt}}\left(u\right)$. Thus, $b_{1,T}^{\mathrm{opt}}$ and $\overline{b}_{2,T}^{\mathrm{opt}}$ need not be the same as the bandwidths $(\widetilde{b}_{1,T}^{\mathrm{opt}},\,\widetilde{b}_{2,T}^{\mathrm{opt}})$ that jointly minimize the global asymptotic MSE,

Theorem (ref)-(ii) states that, under the condition $b_{2,T}^{2}/b_{1,T}^{q}\rightarrow0$, the bias only depends on the smoothing over lagged autocovariances but not on $b_{2,T}$. Then, the global solution $\widetilde{b}_{2,T}^{\mathrm{opt}}$ would be trivial: $b_{2,T}$ affects the MSE only through the variance term and optimality requires to set the bandwidth as large as possible. In contrast, the MSE criterion (ref) based on the MSE given in Theorem (ref)-(iii) leads to a unique solution which can be obtained analytically. Under the condition $b_{2,T}^{2}/b_{1,T}^{q}\rightarrow\nu<\infty$, Belotti et al. belotti/casini/catania/grassi/perron_HAC_Sim_Bandws\nocite{belotti/casini/catania/grassi/perron:2020} determined the bandwidths $(\widetilde{b}_{1,T}^{\mathrm{opt}},\,\widetilde{b}_{2,T}^{\mathrm{opt}})$ that jointly minimize (ref). They showed that $\widetilde{b}_{1,T}^{\mathrm{opt}},\,\widetilde{b}_{2,T}^{\mathrm{opt}}=O(T^{-1/6})$ while the optimal bandwidths $(b_{1,T}^{\mathrm{opt}},\,\overline{b}_{2,T}^{\mathrm{opt}})$ from (ref) satisfy $b_{1,T}^{\mathrm{opt}}=O(T^{-4/25})$ and $\overline{b}_{2,T}^{\mathrm{opt}}=O(T^{-1/5})$. Thus, the criterion (ref) leads to a slightly shorter block length relative to the global criterion (ref) (i.e., $T\overline{b}_{2,T}^{\mathrm{opt}}<T\widetilde{b}_{2,T}^{\mathrm{opt}}$). A shorter bock length is beneficial if there is substantial nonstationarity and implies less sensitivity to low frequency contamination from not properly accounting for nonstationarity {[}cf. Casini et al. casini/perron_Low_Frequency_Contam_Nonstat:2020{]}. For a throughout comparison between the two criteria see Belotti et al. belotti/casini/catania/grassi/perron_HAC_Sim_Bandws.\nocite{belotti/casini/catania/grassi/perron:2020}

Optimal $K_{2}\left(\cdot\right)$ and $b_{2,T}$

Let $F\left(K_{2}\right)\triangleq\int_{0}^{1}K_{2}^{2}\left(x\right)dx$, $H\left(K_{2}\right)=(\int_{0}^{1}x^{2}K_{2}\left(x\right)dx)^{2}$, and for any $k\in\mathbb{Z}$,

The optimal kernel $K_{2}^{\mathrm{opt}}\left(x\right)$ is a transformation of the Epanechnikov kernel. Optimality of quadratic kernels under a MSE criterion has been shown in many contexts {[}cf. epanechnikov:69 and priestley:85{]}. The optimal bandwidth sequence decreases at rate $T^{-1/5}$ which is the same optimal rate derived in the context of parameter estimation of locally stationary processes {[}see e.g., Dahlhaus/Giraitis:98{]}. The term $D_{1}\left(u_{0}\right)$ is due to nonstationary, while the term $D_{2}\left(u_{0}\right)$ measures the variability of $\widehat{c}_{T}\left(u_{0},\,k\right)$. The bandwidth $b_{2,T}^{\mathrm{opt}}$ converges to zero at a slower rate as the process becomes closer to stationary (i.e., as the square root of $D_{1}\left(u_{0}\right)$ decreases).

Optimal $K_{1}\left(\cdot\right)$

We next determine the optimal kernel $K_{1}$ and the optimal bandwidth sequence $b_{1,T}$ given any $K_{2}$ and any $b_{2,T}$ of order $O(T^{-1/5})$, i.e., the same order of $b_{2,T}^{\mathrm{opt}}\left(u\right)$ for any $u\in\left[0,\,1\right]$. Let $\widehat{J}_{T}^{\mathrm{QS}}$ denote $\widehat{J}_{T}$ when the latter is based on the QS kernel. For some results below, we consider a subset of $\boldsymbol{K}_{1}$. Let $\boldsymbol{\widetilde{K}}_{1}=\{K_{1}\left(\cdot\right)\in\boldsymbol{K}_{1}|\,\widetilde{K}\left(\omega\right)\geq0\,\forall\,\omega\in\mathbb{R}\}\bigr\}$ where $\widetilde{K}\left(\omega\right)=\left(2\pi\right)^{-1}\int_{-\infty}^{\infty}K_{1}\left(x\right)e^{-ix\omega}dx.$ The function $\widetilde{K}\left(\omega\right)$ is referred to as the spectral window generator. The set $\widetilde{\boldsymbol{K}}_{1}$ contains all kernels $K_{1}$ that necessarily generate positive semidefinite estimators in finite samples.

We adopt the notation $\widehat{J}_{T}\left(b_{1,T}\right)=\widehat{J}_{T}\left(b_{1,T},\,b_{2,T},\,K_{2}\right)$ to denote the estimator $\widehat{J}_{T}$ that uses $b_{1,T},\,b_{2,T}=\overline{b}_{2,T}^{\mathrm{opt}}+o\left(T^{-1/5}\right)$ and $K_{2}^{\mathrm{}}\left(\cdot\right)$. We then compare two kernels $K_{1}$ using comparable\textcolor{red}{ }bandwidths $b_{1,T}$ which are defined as follows. Given $K_{1}\left(\cdot\right)\in\widetilde{\boldsymbol{K}}_{1}$, the QS kernel $K_{1}^{\mathrm{QS}}\left(\cdot\right)$, and a bandwidth sequence $\left\{ b_{1,T}\right\} $ to be used with the QS kernel, define a comparable bandwidth sequence $\left\{ b_{1,T,K_{1}}\right\} $ for use with $K_{1}\left(\cdot\right)$ such that both kernel/bandwidth combinations have the same asymptotic variance when scaled by the same factor $Tb_{1,T}b_{2,T}$. This means that

This definition yields $b_{1,T,K_{1}}=b_{1,T}/(\int K_{1}^{2}\left(x\right)dx)$ and $b_{1,T,\mathrm{QS}}=b_{1,T}$ since $\int\left(K_{1}^{\mathrm{QS}}\right)^{2}\left(x\right)dx=1$.

The requirement $\int_{0}^{1}||f^{\left(2\right)}\left(u,\,0\right)||du<\infty$ is not stringent and it reduces to the one used by andrews:91 when $\left\{ V_{t,T}\right\} $ is stationary. If $\int_{0}^{1}||f^{\left(q\right)}\left(u,\,0\right)||du<\infty$ only for some $1\leq q<2$, one can show that any kernel with $K_{1,q}=0$ has smaller asymptotic MSE than a kernel with $K_{1,q}>0$. The QS, Parzen, and Tukey-Hanning kernels have $K_{1,q}=0$ for $1\leq q<2$, whereas the Bartlett has $K_{1,q}>0$ for $1\leq q<2$. Thus, the asymptotic superiority of the former kernels over the Bartlett kernel holds even if $\int_{0}^{1}||f^{\left(q\right)}\left(u,\,0\right)||du<\infty$ only for $1\leq q<2$.

Optimal Predetermined Bandwidth Sequence $b_{1,T}$

We now present the predetermined bandwidth sequence that minimizes the asymptotic MSE given $b_{2,T}=O(b_{2,T}^{\mathrm{opt}})$ and $K_{2}=K_{2}^{\mathrm{opt}}.$ This optimality result applies to each kernel $K_{1}\left(\cdot\right)\in\boldsymbol{K}_{1}$ for which $K_{1,q}\in\left(0,\,\infty\right)$ for some $q\in\left(0,\,\infty\right)$. Thus, most commonly used kernels are allowed with the exception of the truncated kernel. Let

The optimal bandwidth is $b_{1,T}^{\mathrm{opt}}=(2qK_{1,q}^{2}\phi\left(q\right)Tb_{2,T}^{\mathrm{\mathrm{opt}}}/(\int K_{1}^{2}\left(y\right)dy\int_{0}^{1}K_{2}^{2}\left(x\right)dx))^{-1/\left(2q+1\right)}$, where $\phi\left(q\right)$ is a function of the\textcolor{red}{ }unknown spectral density $f\left(\cdot,\,\cdot\right)$. Hence, the optimal bandwidth $b_{1,T}^{\mathrm{opt}}$ is unknown in practice, and we consider data-dependent estimates of $\phi\left(q\right)$ in Section (ref).

In Corollary (ref), $q=2$ for the QS kernel and so $b_{1,T}^{\mathrm{opt}}=0.6584(\phi\left(2\right)Tb_{2,T}^{\mathrm{\mathrm{opt}}})^{-1/5}(\int_{0}^{1}K_{2}^{2}\left(y\right)dy)^{1/5}.$ For $K_{2}\left(y\right)=K_{2}^{\mathrm{opt}}\left(y\right),$ the latter reduces to,

The optimal bandwidth is of order $T^{-4/25}.$ Thus, the optimal bandwidth sequence decreases to zero at a slower rate than the optimal bandwidth sequence for the QS kernel-based HAC estimator of andrews:91, for which the rate is of order $T^{-1/5}$. The slower rate is due to the fact that our estimator smooths the spectrum over time through $K_{2}\left(\cdot\right)$ and this restricts the smoothing of $K_{1}\left(\cdot\right)$. In particular, the optimal choice of $b_{1,T}$ hinges on the degree of nonstationary through $b_{2,T}^{\mathrm{\mathrm{opt}}}$. The more nonstationary are the data, the smaller is $b_{2,T}^{\mathrm{\mathrm{opt}}}$ and the large is $b_{1,T}^{\mathrm{\mathrm{opt}}}$ which means that less weight is given to $\widehat{\Gamma}\left(k\right)$ for $k\neq0$. In contrast, the optimal choice of $b_{1,T}$ for the methods proposed in the literature is independent of the degree of nonstationarity. When $b_{1,T}$ and $b_{2,T}$ are chosen optimally, the convergence rate from Theorem (ref) reduces to $T^{8/25}$. Thus, the rate is slower than the corresponding one for the QS kernel HAC estimator considered in andrews:91. However, it is misleading to compare our DK-HAC estimator with the classical HAC estimators only on the basis of the rate of convergence. In fact, the DK-HAC estimators account flexibly for nonstationarity and are robust to low frequency contamination induced by nonstationarity/misspecification whereas the classical HAC estimators are not in general {[}cf. Casini et al. casini/perron_Low_Frequency_Contam_Nonstat:2020{]}.

Choice of $n_{T}$

Our MSE analysis does not indicate an optimal value for $n_{T}$. It only suggests growth rate bounds. When $K_{1}^{\mathrm{QS}}$ is used, $n_{T}$ cannot grow faster than $T^{2/3}$. We set $n_{T}=T^{0.66}$ for the $\mathrm{QS}$ kernel. That is, we choose $n_{T}$ to be the largest possible value allowed by the condition. Our sensitivity analysis (not reported) suggests that choosing a smaller $n_{T}$ might result in excessive overlapping of regimes when the process is SLS (i.e., $m_{0}>0$). See Belotti et al. belotti/casini/catania/grassi/perron_HAC_Sim_Bandws for more details.

Data-Dependent Bandwidths

In this section we consider estimators $\widehat{J}_{T}$ that use bandwidths $b_{1,T}$ and $b_{2,T}$ whose values are determined via data-dependent methods. We use the “plug-in” method which is characterized by plugging-in estimates of unknown quantities into an asymptotic formula for an optimal bandwidth parameter (i.e., the expressions for $b_{1,T}^{\mathrm{opt}}$ and $b_{2,T}^{\mathrm{opt}}$ from Section (ref)). Section (ref) explains how to construct the automatic bandwidths while Section (ref) presents the corresponding theoretical results.

Implementation

Let us begin with $b_{1,T}^{\mathrm{opt}}$ and then move to $b_{2,T}^{\mathrm{opt}}$. The first step for the construction of data-dependent bandwidth parameters is to specify $p$ univariate parametric models for the elements of $V_{t}=(V_{t}^{\left(1\right)},\ldots,\,V_{t}^{\left(p\right)})'$. The second step involves the estimation of the parameters. In our context, the logical estimation methods to use are local (weighted) least-squares (LS) (i.e., LS method applied to rolling windows) and nonparametric kernel methods. In a third step, we replace the unknown parameters in $\phi\left(q\right)$ with corresponding estimates. Such estimate $\widehat{\phi}\left(q\right)$ is then substituted into the expression for $b_{1,T}^{\mathrm{opt}}$ to yield the data-dependent bandwidth $\widehat{b}_{1,T}$:

where $\widehat{\overline{b}}_{2,T}=\left(n_{T}/T\right)\sum_{r=1}^{\left\lfloor T/n_{T}\right\rfloor -1}\widehat{b}_{2,T}\left(rn_{T}/T\right)$. $\widehat{\overline{b}}_{2,T}$ is an average of the estimates $\widehat{b}_{2,T}\left(\cdot\right)$. Since $b_{2,T}$ depends on $u$, it is more efficient to estimate it for each block as its optimal value can change over time. In practice, a reasonable candidate for an approximating parametric model is the class of first order autoregressive {[}AR(1){]} models for $\{V_{t}^{\left(r\right)}\},\,r=1,\ldots,\,p$ (with different parameters for each $r$) or a first order vector autoregressive {[}VAR(l){]} model for $\{V_{t}\}$. These classes were also used by andrews:91. However, in our context it is reasonable to allow the parameters of the AR(1) model to be time-varying. For parsimony, we consider a time-varying AR(1) with no breaks in $f\left(u,\,\omega\right)$, i.e., $V_{t}^{\left(r\right)}=a_{1}\left(t/T\right)V_{t-1}^{\left(r\right)}+u_{t}^{\left(r\right)},$ where the $u_{t}^{\left(r\right)}$ need not be independent across $r$.

The use of $p$ univariate parametric models requires a simple form for the weight matrix $W$. In particular, $W$ has to be a diagonal matrix which in turn implies that $\phi\left(q\right)$ reduces to

The usual choice is $W^{\left(r,r\right)}=1$ for all $r$ except that which corresponds to an intercept for which it is set to zero. An estimate of $f^{\left(r,r\right)}\left(u,\,0\right)$ is $\widehat{f}^{\left(r,r\right)}\left(u,\,0\right)=\left(2\pi\right)^{-1}(\widehat{\sigma}^{\left(r\right)}\left(u\right))^{2}(1-\widehat{a}_{1}^{\left(r\right)}\left(u\right))^{-2}$ while $f^{\left(2\right)\left(r,r\right)}\left(u,\,0\right)$ can be estimated by $\widehat{f}^{\left(2\right)\left(r,r\right)}\left(u,\,0\right)=3\pi^{-1}$ $((\widehat{\sigma}^{\left(r\right)}\left(u\right))^{2}\widehat{a}_{1}^{\left(r\right)}\left(u\right))(1-\widehat{a}_{1}^{\left(r\right)}\left(u\right))^{-4}$ where $\widehat{a}_{1}^{\left(r\right)}\left(u\right)$ and $\widehat{\sigma}^{\left(r\right)}\left(u\right)$ are the LS estimates computed using local data to the left of $u=t/T$:

where $n_{2,T}\rightarrow\infty$. Then, for the QS kernel $K_{1}$,

For most of the results below we can take $n_{3,T}=n_{2,T}=n_{T}.$ After plugging-in $\widehat{\phi}\left(2\right)$ into the formula (ref), we have $\widehat{b}_{1,T}=0.6828(\widehat{\phi}\left(2\right)T\widehat{\overline{b}}_{2,T})^{-1/5}.$

We now propose a data-dependent procedure for $b_{2,T}\left(u_{r}\right),$ where $u_{r}=rn_{T}/T$ for $r=1,\ldots,\,\left\lfloor \left(T-n_{T}\right)/n_{T}\right\rfloor $. We assume that the parameters of the approximating time-varying AR(1) models change slowly such that the smoothness of $f\left(\cdot,\,\omega\right)$ and thus of $c\left(\cdot,\,k\right)$ is the same as the one that would arise if $a_{1}\left(u\right)=0.8\left(\cos1.5+\cos4\pi u\right)$ and $\sigma\left(u\right)=\sigma=1$ for all $u\in\left[0,\,1\right]$ {[}cf. dahlhaus:12{]}. The reason for imposing this condition is that it is otherwise difficult to estimate $\left(\partial^{2}/\partial u^{2}\right)c\left(u,\,k\right)$, which enters $D_{1}\left(u\right)$, from the data. Under the above specification, the exact expression of $D_{1}\left(u\right)$ can be computed analytically:

An estimate of $D_{1}\left(u\right)$ is given by

where $\left[S_{\omega}\right]$ is the cardinality of $S_{\omega}$ and $\omega_{s+1}>\omega_{s}$ with $\omega_{1}=-\pi,\,\omega_{\left[S_{\omega}\right]}=\pi.$ In our simulations we use $S_{\omega}=\left\{ -\pi,\,-3,\,-2,\,-1,\,0,\,1,\,2,\,3,\,\pi\right\} $. Note that we have computed $\widehat{D}_{1}\left(u\right)$ for $k=0$ because it makes the computation simpler. Further, this is consistent with our sequential MSE criterion because $k=0$ is the only lag for which $K_{1}(0)=1$ for all $K_{1}$ so that the choice of $K_{1}$ does not influence $b_{2}^{\mathrm{opt}}\left(\cdot\right)$. It remains to derive an estimate of $D_{2}\left(u\right)$ since $F\left(K_{2}\right)$ and $H\left(K_{2}\right)$ can be computed for a given $K_{2}\left(\cdot\right)$. We assume that the innovations of the approximating time-varying AR(1) model satisfy $\mathbb{E}(u_{t}^{\left(r\right)})=0$, $\mathbb{E}((u_{t}^{\left(r\right)})^{2})=\sigma^{2}$ and $\mathbb{E}((u_{t}^{\left(r\right)})^{4})=3\sigma^{4}$ so that $D_{3}\left(u\right)=0$ for all $u\in\left(0,\,1\right)$. That is, the term involving the cumulant drops from $b_{2}^{\mathrm{opt}}\left(u\right)$. In practice this is convenient because it is complex to deal with consistent estimation of cumulant terms. Note also that $D_{3}\left(u\right)=0$ if $u_{t}$ is Gaussian. Since $c\left(u,\,k\right)$ can be consistently estimated by $\widehat{c}_{T}\left(u,\,k\right)$, an estimate of $D_{2}\left(u\right)$ is given by \[ \widehat{D}_{2}\left(u_{0}\right)\triangleq p^{-1}\sum_{r=1}^{p}\sum_{l=-\left\lfloor T^{4/25}\right\rfloor }^{\left\lfloor T^{4/25}\right\rfloor }\widehat{c}_{T}^{\left(r,r\right)}\left(u_{0},\,l\right)\left[2\widehat{c}_{T}^{\left(r,r\right)}\left(u_{0},\,l\right)\right], \] where the number of summands grows at the same rate as $(b_{1,T}^{\mathrm{opt}})^{-1}$; a different choice is allowed as long as it grows at a slower rate than $T^{2/5}$. Hence, the estimate of the optimal $b_{2,T}$ is given by

Theoretical Results

Next, we establish consistency, rate of convergence and asymptotic MSE results for the estimator $\widehat{J}_{T}(\widehat{b}_{1,T},\,\widehat{\overline{b}}_{2,T})$ that uses the data-dependent bandwidths $\widehat{b}_{1,T}$ and $\widehat{\overline{b}}_{2,T}$. As in andrews:91, we need to restrict the class of admissible kernels to the following class:

Let $\widehat{\theta}$ denote the estimator of the parameter of the approximate (time-varying) parametric model(s) introduced above. For example, with univariate AR(1) approximating parametric models, $\widehat{\theta}=(\int_{0}^{1}\widehat{a}_{1}^{\left(1\right)}\left(u\right)du,\,\int_{0}^{1}(\widehat{\sigma}^{\left(1\right)}\left(u\right))^{2}du,\ldots,\,\int_{0}^{1}\widehat{a}_{1}^{\left(p\right)}\left(u\right)du,\,\int_{0}^{1}(\widehat{\sigma}^{\left(p\right)}\left(u\right))^{2}du)'$. Let $\theta^{*}$ denote the probability limit of $\widehat{\theta}$. $\widehat{\phi}\left(q\right)$ is the value of $\phi\left(q\right)$ with $\widehat{\theta}$ instead of $\theta$. Its probability limit is denoted by $\phi_{\theta^{*}}$.

Parts (i)-(ii) of Assumption (ref) are the nonparametric analogue to Assumption E and F, respectively, in andrews:91. Part (iii) is satisfied if $\left\{ V_{t}\right\} $ is strong mixing with mixing numbers that are less stringent than those sufficient for the cumulant condition in Assumption (ref)-(i). Part (iv) and (vi) extend (i)-(ii) to $\widehat{D}_{1}$ and $\widehat{D}_{2}$. Part (v) is needed to apply the convergence of Riemann sums. Part (vi) follows from the asymptotic results about $\widehat{c}_{T}\left(u,\,k\right)$. Part (vii) requires $K_{2}$ to satisfy Lipschitz continuity. Note that $\phi_{\theta^{*}}$ coincides with the optimal value $\phi\left(q\right)$ only when the approximate parametric model indexed by $\theta^{*}$ corresponds to the true data-generating mechanism.

Let $b_{\theta_{1},T}=(2qK_{1,q}^{2}\phi_{\theta^{*}}T\overline{b}_{\theta_{2},T}/\int K_{1}^{2}\left(y\right)dy\int_{0}^{1}K_{2}^{2}\left(x\right)dx)^{-1/\left(2q+1\right)},$ where $\overline{b}_{\theta_{2},T}\triangleq\int_{0}^{1}b_{2,T}^{\mathrm{opt}}\left(u\right)du$. The asymptotic properties of $\widehat{J}_{T}(\widehat{b}_{1,T},\,\widehat{\overline{b}}_{2,T})$ are shown to be equivalent to those of $\widehat{J}_{T}(b_{\theta_{1},T},\,b_{\theta_{2},T})$.

When the chosen parametric model indexed by $\theta$ is correct, it follows that $\phi_{\theta^{*}}=\phi\left(q\right)$ and $\widehat{\phi}\left(q\right)\overset{\mathbb{P}}{\rightarrow}\phi\left(q\right)$. The theorem then implies that $\widehat{J}_{T}(\widehat{b}_{1,T},\,\widehat{\overline{b}}_{2,T})$ exhibits the same optimality properties presented in Theorem (ref) and Corollary (ref). We omit the details.

Small-Sample Evaluations

We conduct a Monte Carlo analysis to evaluate the properties of HAR inference based on the HAC estimator $\widehat{J}_{T}$. We consider HAR tests in the linear regression model as well as HAR tests used in the forecast evaluation literature, namely the diebold/mariano:95's diebold/mariano:95 test\nocite{diebold/mariano:95} and the forecast breakdown test of giacomini/rossi:09. The linear regression models have an intercept and a stochastic regressor. We focus on the $t$-statistics $t_{r}=\sqrt{T}(\widehat{\beta}^{\left(r\right)}-\beta_{0}^{\left(r\right)})/\sqrt{\widehat{J}_{X,T}^{\left(r,r\right)}}$ where

is a consistent estimate of the limit of $\mathrm{Var}(\sqrt{T}(\widehat{\beta}-\beta_{0}))$ and $r=1,\,2$. $t_{1}$ is the $t$-statistic for the parameter associated with the intercept while $t_{2}$ is associated with the stochastic regressor $x_{t}$. Results for the $F$-test are qualitatively similar {[}see casini_diss{]}. Six basic regression models are considered. We run a $t$-test on the intercept in model M1 and M5 whereas a $t$-test on the coefficient of $x_{t}$ is run in model M2-M4 and M6. The models are based on,

for the $t$-test on the intercept (i.e., $t_{1}$) and

for the $t$-test on $\beta_{0}^{\left(2\right)}$ (i.e., $t_{2}$) where $\delta=0$ under the null hypotheses. In Model M1 $e_{t}=0.5e_{t-1}+u_{t},\,u_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(0,\,0.5\right),$ $x_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(1,\,1\right)$, $\beta_{0}^{\left(1\right)}=0$ and $\beta_{0}^{\left(2\right)}=1.$\footnote{For the results with AR coefficient 0.9 see Table 1 in casini/perron_PrewhitedHAC and footnote (ref) below.} Model M2 involves $e_{t}=0.8e_{t-1}+u_{t},\,u_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(0,\,1\right),$ $x_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(1,\,1\right)$, and $\beta_{0}^{\left(1\right)}=\beta_{0}^{\left(2\right)}=0.$ In Model M3 we have segmented locally stationary errors $e_{t}=\rho_{t}e_{t-1}+u_{t},\,u_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(0,\,1\right),\,\rho_{t}=\max\left\{ 0,\,-1\left(\cos\left(1.5-\cos\left(5t/T\right)\right)\right)\right\} $ for $t<4T/5$ and $e_{t}=0.9e_{t-1}+u_{t},\,u_{t}\sim\mathrm{\mathrm{i.i.d.\,}}\mathscr{N}\left(0,\,1\right)$ for $t\geq4T/5,$ and $x_{t}=0.4x_{t-1}+u_{X,t},\,u_{X,t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(0,\,1\right)$. Note that $\rho_{t}$ varies smoothly between 0 and 0.8071. Model M4 involves some misspecification that induces nonstationarity in the errors,

where $e_{t}=\rho_{t}e_{t-1}+u_{t},\,u_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(0,\,1\right),\,\rho_{t}$ as in M3, $x_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(1,\,1\right)$, and $w_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(2,\,1\right)$ independent from $x_{t}$. Model M5 involves misspecification under $H_{1}$ via a smooth change in the coefficient $\beta_{0}^{\left(2\right)}$ toward the end of the sample. This situation is very common in practice and it is motivated by the model for the variable “cay” from bianchi/lettau/ludvigson:18 (cf. Figure 3 in their paper). The model is given by

where $d_{t}=1.5\delta\left(t-4.5T/5\right)/T$, $e_{t}=\rho_{t}e_{t-1}+u_{t},\,u_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(0,\,1\right),\,\rho_{t}=0.8(\cos(1.5-\cos(t/$ $2T)))$ for $t\in\left\{ 1,\ldots,\,T/2-1\right\} \cup\left\{ T/2+T/4+1,\ldots,\,T\right\} $ and $e_{t}=0.2e_{t-1}+2u_{t},\,u_{t}\sim\mathrm{\mathrm{i.i.d.\,}}\mathscr{N}\left(0,\,1\right)$ for $T/2\leq t\leq T/2+T/4$, and $x_{t}=2+0.5x_{t-1}+u_{X,t},\,u_{X,t}\sim\mathrm{i.i.d.}\,\mathscr{N}\left(0,\,1\right)$. That is, $\rho_{t}$ varies smoothly between 0 and 0.7021. Model M6 is given by (ref) where $e_{t}=\rho_{t}e_{t-1}+u_{t},\,u_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(0,\,1\right),\,\rho_{t}=\max\left\{ 0,\,0.3\left(\cos\left(1.5-\cos\left(t/5T\right)\right)\right)\right\} $ for $t\in\left\{ 1,\ldots,\,T/2-1\right\} \cup\{T/2+4,\ldots,$ $T-16\}$ and $e_{t}=0.99e_{t-1}+2u_{t},\,u_{t}\sim\mathrm{\mathrm{i.i.d.\,}}\mathscr{N}\left(0,\,1\right)$ for $T/2\leq t\leq T/2+3$ and $e_{t}=0.9e_{t-1}+2u_{t},\,u_{t}\sim\mathrm{\mathrm{i.i.d.\,}}\mathscr{N}\left(0,\,1\right)$ for $T-15\leq t\leq T$, and $x_{t}\sim\mathrm{i.i.d.}\,\mathscr{N}\left(1,\,1\right)$. Note that $\rho_{t}\in\left[0,\,0.2633\right]$.

Next, we move to the forecast evaluation tests. The Diebold-Mariano test statistic is defined as $t_{\mathrm{DM}}\triangleq\sqrt{T_{n}}\,\overline{d}_{L}/\sqrt{\widehat{J}_{d_{L},T}}$, where $\overline{d}_{L}$ is the average of the loss differentials between two competing forecast models, $\widehat{J}_{d_{L},T}$ is an estimate of the asymptotic variance of the the loss differential series and $T_{n}$ is the number of observations in the out-of-sample. Throughout we use the quadratic loss. In model M7, we consider an out-of-sample forecasting exercise with a fixed scheme where, given a sample of $T$ observations, $0.5T$ observations are used for the in-sample and the remaining half is used for prediction. The true model for the target variable is given by $y_{t}=\beta_{0}^{\left(1\right)}+\beta_{0}^{\left(2\right)}x_{t-1}^{(0)}+e_{t}$ where $x_{t-1}^{(0)}\sim\mathrm{i.i.d.}\,\mathscr{N}\left(1,\,1\right)$, $e_{t}=0.3e_{t-1}+u_{t}$ with $u_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(0,\,1\right)$ and we set $\beta_{0}^{\left(1\right)}=\beta_{0}^{\left(2\right)}=1.$ The two competing models both involve an intercept but differ on the predictor used in place of $x_{t}^{(0)}$. The first forecast model uses $x_{t}^{(1)}$ while the second uses $x_{t}^{(2)}$ where $x_{t}^{(1)}$ and $x_{t}^{(2)}$ are independent $\mathrm{i.i.d.}\,\mathscr{N}\left(1,\,1\right)$ sequences, both independent from $x_{t}^{(0)}$. Each forecast model generates a sequence of $\tau\left(=1\right)$-step ahead out-of-sample losses $L_{t}^{(i)}$ $\left(i=1,\,2\right)$ for $t=T/2+1,\ldots,\,T-\tau.$ Then $d_{t}\triangleq L_{t}^{(2)}-L_{t}^{(1)}$ denotes the loss differential at time $t$. The Diebold-Mariano test rejects the null of equal predictive ability when (after normalization) $\overline{d}$ is sufficiently far from zero.

Finally, we consider model M8 which we use to investigate the performance of a $t$-test for forecast breakdown {[}cf. giacomini/rossi:09{]}. Suppose we want to forecast a variable $y_{t}$ following the equation $y_{t}=\beta_{0}^{\left(1\right)}+\beta_{0}^{\left(2\right)}x_{t-1}+e_{t}$ where $x_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(1,\,1.5\right)$ and $e_{t}=0.3e_{t-1}+u_{t}$ with $u_{t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(0,\,1\right)$. For a given forecast model and forecasting scheme, the test of giacomini/rossi:09 (GR) detects a forecast breakdown when the average of the out-of-sample losses differs significantly from the average of the in-sample losses. The in-sample is used to obtain estimates of $\beta_{0}^{\left(1\right)}$ and $\beta_{0}^{\left(2\right)}$ which are in turn used to construct out-of-sample forecasts $\widehat{y}_{t}=\widehat{\beta}_{0}^{\left(1\right)}+\widehat{\beta}_{0}^{\left(2\right)}x_{t-1}$. We set $\beta_{0}^{\left(1\right)}=\beta_{0}^{\left(2\right)}=1.$ We consider a fixed forecasting scheme and one-step ahead forecasts. The GR's (2009) test statistic is defined as $t^{\mathrm{GR}}\triangleq\sqrt{T_{n}}\overline{SL}/\sqrt{\widehat{J}_{SL}}$ where $\overline{SL}\triangleq T_{n}^{-1}\sum_{t=T_{m}}^{T-\tau}SL_{t+\tau}$, $SL_{t+\tau}$ is the surprise loss at time $t+\tau$, i.e., the difference between the time $t+\tau$ out-of-sample loss and in-sample loss, $SL_{t+\tau}=L_{t+\tau}-\overline{L}_{t+\tau}$, $T_{n}$ is the sample size in the out-of-sample, $T_{m}$ is the sample size in the in-sample and $\widehat{J}_{SL}$ is an HAC estimator. We restrict attention to $\tau=1.$

Throughout our study we consider the following LRV estimators: $\widehat{J}_{T}$ with automatic bandwidths; $\widehat{J}_{T}$ with automatic bandwidths and the prewhitening of casini/perron_PrewhitedHAC; Andrews's (1991) HAC estimator with automatic bandwidth; Andrews's (1991) HAC estimator with automatic bandwidth and the prewhitening procedure of andrews/monahan:92; Newey and West's (1987) HAC estimator with the automatic bandwidth as proposed in newey/west:94; Newey and West's (1987) HAC estimator with the automatic bandwidth as proposed in newey/west:94 and the prewhitening procedure; Newey-West with the fixed-$b$ method of Kiefer/vogelsang/bunzel:00.\footnote{To save space, we do not report results for the Empirical Weighted Periodogram (EWP) or Empirical Weighted Cosine (EWC) of Lazarus et al. lazarus/lewis/stock:17 and Lazarus et al. lazarus/lewis/stock/watson:18, respectively. Their performance is similar to the method of Kiefer/vogelsang/bunzel:00. The LRV estimator of Kiefer/vogelsang/bunzel:00 leads to HAR tests that have better size control. casini/perron_PrewhitedHAC showed that EWC leads to oversized tests when there is strong dependence in the data relative to the fixed-$b$ method of Kiefer/vogelsang/bunzel:00 and to the prewhitened DK-HAC. The power properties of tests normalized by the EWP and EWC are similar to those using the method of Kiefer/vogelsang/bunzel:00.} casini/perron_PrewhitedHAC proposed three forms of prewhitening: (1) $\widehat{J}_{T,\mathrm{pw},1}$ uses a stationary model to whiten the data; (2) $\widehat{J}_{T,\mathrm{pw},\mathrm{SLS}}$ uses a nonstationary model to whiten the data; (3) $\widehat{J}_{T,\mathrm{pw},\mathrm{SLS},\mu}$ is the same as $\widehat{J}_{T,\mathrm{pw},\mathrm{SLS}}$ but it adds a time-varying intercept in the VAR to whiten the data. For model M7 we also report results using $\widehat{J}_{T}$ and $\widehat{J}_{T,\mathrm{pw},\mathrm{SLS}}$ with the pre-test for breaks in the spectrum as developed in casini:change-point-spectra. We do not report the results for the pre-test for model M1-M6 and M8 because they are equivalent to those without the pre-test.

For all versions of $\widehat{J}_{T}$ we use $K_{1}^{\mathrm{opt}}$ and $K_{2}^{\mathrm{opt}}$. We set $n_{T}=T^{0.66}$ as explained in Section (ref) and $n_{2,T}=n_{3,T}=n_{T}.$ We consider the following sample sizes for M1-M6: $T=200,\,400$. Simulation results for additional data-generating processes involving ARMA, ARCH and heteroskedastic errors are not discussed here because the results are qualitatively equivalent {[}see, e.g., casini_diss and casini/perron_PrewhitedHAC{]}. The significance level is $\alpha=0.05$ throughout the study.

Empirical Sizes of HAR Inference Tests

Table (ref)-(ref) report the rejection rates for model M1-M8. We begin with the $t$-test in the linear regression models. As a general pattern, we confirm previous evidence that the Newey-West's (1987) and Andrews' (1991) HAC estimators lead to $t$-tests that are oversized when the data are stationary {[}cf. model M1-M2{]}. The same problem occurs for the Newey-West (1987) HAC estimator using the usual “rule” to determine the number of lags (not reported). For extreme temporal dependence, simulations in casini/perron_PrewhitedHAC showed that the size distortions can be even larger especially for the $t$-test on the intercept. Prewhitening is often effective in helping the HAC estimators to better control the size under stationarity. However, the simulation results in casini/perron_PrewhitedHAC and in the literature show that the prewhitened HAC estimators can lead to oversized tests when there is high serial dependence. The rejection rates of tests normalized by the Newey-West estimator with fixed-$b$ are the most accurate in model M1-M2 for $T=200$. Overall, the results in the literature along with those in casini_diss and casini/perron_PrewhitedHAC showed that under stationarity the original fixed-$b$ method of KVB is in general the least oversized across different degrees of dependence among all existing methods. Table (ref) shows that for the $t$-test on the intercept the non-prewhitened DK-HAC leads to HAR tests that are oversized while they are accurate for the $t$-test on the coefficient on the stochastic regressor. The table also shows that the prewhitened DK-HAC estimators are competitive with the KVB's fixed-$b$ in controlling the size. $\widehat{J}_{T,\mathrm{pw},1}$ is the most accurate among the DK-HAC estimators. Since $\widehat{J}_{T,\mathrm{pw},1}$ uses a stationarity VAR model to whiten the data, it works better than $\widehat{J}_{T,\mathrm{pw},\mathrm{SLS}}$ and $\widehat{J}_{T,\mathrm{pw},\mathrm{SLS},\mu}$ when stationarity actually holds which is consistent with the results of Table (ref).

Turning to nonstationarity, Table (ref) casts concerns about the performance of existing methods in this context. For both model M3 and M4, existing LRV estimators lead to HAR tests that have either size equal or close to zero. The methods that use long bandwidths (i.e., many lags) such as KVB's fixed-$b$ suffer most from this problem relative to the classical HAC estimators. This is consistent with the argument in casini/perron_Low_Frequency_Contam_Nonstat:2020 who showed analytically that nonstationarity induces a positive bias for each sample autocovariance. That bias is constant across different lags. Since existing LRV estimators are weighted sum of sample autocovariances, the larger the bandwidth (i.e., the more lagged autocovariances are included) the larger the positive bias. Thus, LRV estimators are inflated and HAR tests have rejection rates lower than the significance level. This mechanism has consequences for power as well, as we show below that traditional HAR tests have low power. In model M3-M4 the non-prewhitened DK-HAC and the prewhitened DK-HAC (except $\widehat{J}_{T,\mathrm{pw},1}$) perform well. $\widehat{J}_{T,\mathrm{pw},1}$ suffers from the same problem as the existing estimators because it uses stationarity and when this is violated its performance is affected. In model M5, the classical HAC estimators yield HAR tests that are oversized. Also the non-prewhitened DK-HAC is oversized. In contrast, the KVB's fixed-$b$ and the prewhitened DK-HAC have rejection rates close to the significance level. In model M6, the KVB's fixed-$b$ HAR tests tend to be undersized whereas the HAC and DK-HAC estimators lead to tests that control the size more accurately.

Turning to the HAR tests for forecast evaluations, Table (ref) shows that for model M7 the KVB's fixed-$b$ HAR test has size essentially equal to zero while the classical HAC estimators yield HAR tests that are somewhat oversized. In contrast, the tests normalized by the prewhitened DK-HAC estimators have most accurate rejection rates. In model M8, the KVB's fixed-$b$ HAR tests are well-sized whereas the classical HAC estimators lead to tests that are severely undersized. The DK-HAC estimators control the size reasonably well.

In summary, the prewhitened DK-HAC estimators yield $t$-tests in regression models with rejection rates that are relatively close to the nominal size. The non-prewhitened DK-HAC can lead to oversized tests for the $t$-tests on the intercept if there is high dependence. Our results confirm the oversize problem of the HAR tests normalized by the classical HAC estimators documented in the literature under stationarity. The Fixed-$b$ HAR tests control the size well when the data are stationary but can show severe undersized issues under nonstationarity, a problem also affecting the tests normalized by the classical HAC estimators. Thus, with regards to size control, the prewhitened DK-HAC estimators are competitive with fixed-$b$ methods under stationarity and they also perform well when the data are nonstationary.

Empirical Power of HAR Inference Tests

For model M1-M6 we report the values of the power in Table (ref)-(ref). The sample size is $T=200$. Power functions for the Diebold-Maraino and for the forecast breakdown test are presented next. For model M1, the non-prewhitened HAC and DK-HAC lead to tests that have the highest power but they were more oversized than the other methods. The KVB's fixed-$b$ LRV leads to $t$-tests that sacrify some power relative to the prewhitened HAC and DK-HAC estimators. In model M2, a similar picture arises. HAR tests normalized by either classical HAC or DK-HAC estimators have similarly good power while HAR tests based on KVB's fixed-$b$ have relatively less power. In model M3, the prewhitening HAC estimators and $\widehat{J}_{T,\mathrm{pw},1}$ (which uses a stationary model to whiten the data) have low power. The best power is achieved by tests normalized by Andrews' (1991) HAC estimator and $\widehat{J}_{T}$, followed by $\widehat{J}_{T,\mathrm{pw},\mathrm{SLS}}$ and $\widehat{J}_{T,\mathrm{pw},\mathrm{SLS},\mu}$. The KVB's fixed-$b$ leads to relatively less power than the latter methods. The Newey-West's (1987) estimator leads to tests that have good power but they were shown to be oversized. Similar comments apply to model M4. Here Andrews' (1991) HAC estimator leads to tests that have better power for small to medium breaks while tests based on $\widehat{J}_{T}$ have better power for large breaks. In model M5, prewhitening HAC estimators and KVB's fixed-$b$ lead to HAR tests that have non-monotonic power and reach zero as $\delta$ increases. This does not occur for the classical HAC estimators which, however, were oversized. HAR tests based on $\widehat{J}_{T,\mathrm{QS}}$, $\widehat{J}_{T,\mathrm{pw},\mathrm{SLS}}$ and $\widehat{J}_{T,\mathrm{pw},\mathrm{SLS},\mu}$ perform best for this model. $\widehat{J}_{T,\mathrm{pw},1}$ results in HAR tests that have lower power relative to the tests based on the other DK-HAC because stationarity is violated. In model M6, all HAR tests enjoy monotonic power with small differences.

Next, let us move to the evaluation of the power properties of the $t$-tests used in the forecasting literature. We begin with the Diebold-Mariano test. For this test, the separation between the null and alternative hypotheses does not depend on the value of a single parameter. Thus, the data-generating mechanism is different from the one under the null. The two competing forecast models are as follows: the first model uses the actual true data-generating process while the second model differs in that in place of $x_{t-1}^{(0)}$ it uses $x_{t-1}^{(2)}=x_{t-1}^{(0)}+u_{X_{2},t}$ for $t\leq3T/4$ and $x_{t-1}^{(2)}=\delta+x_{t-1}^{(0)}+u_{X_{2},t}$ for $t>3T/4$ with $u_{X_{2},t}\sim\mathrm{i.i.d.\,}\mathscr{N}\left(0,\,1\right)$. Evidently, the null hypothesis of equal predictive ability should be rejected by the Diebold-Mariano test whenever $\delta>0$. Table (ref) reports the power for several values of $\delta.$ The HAR tests based on existing estimators have lower power relative to the $\widehat{J}_{T}$ DK-HAC estimators for small values of $\delta$. When we raise $\delta$ the tests based on the HAC estimators of andrews:91 and Newey and West (1987), and KVB's fixed-$b$ method display non-monotonic power gradually converging to zero. In contrast, the DK-HAC estimators lead to tests that have monotonic power that reach and maintain unit power. The only exception is the test based on $\widehat{J}_{T,\mathrm{pw},1}$ that has lower power because stationarity is violated. The table also reports $\widehat{J}_{T}$ and $\widehat{J}_{T,\mathrm{pw},\mathrm{SLS}}$ with the pre-test for breaks in the spectrum {[}cf. casini:change-point-spectra{]} that is used for choosing more efficiently how to split the sample in blocks to compute $\widehat{\Gamma}\left(k\right)$. The pre-test yields HAR tests with higher power while having the same size as the corresponding HAR tests with no pre-test. We have not reported the results with the pre-test for model M1-M6 and M8 because they are the same as with no pre-test.

Finally, we move to the $t$-test of giacomini/rossi:09. The data-generating process under $H_{1}:\,\mathbb{E}\left(\overline{SL}\right)\neq0$ is given by $y_{t}=1+x_{t-1}+\delta x_{t-1}\mathbf{1}\left\{ t>T_{1}^{0}\right\} +e_{t}$, where $x_{t-1}\sim\mathrm{i.i.d.}\,\mathscr{N}\left(1.5,\,1\right)$, $e_{t}=0.3e_{t-1}+u_{t}$, $u_{t}\sim\mathrm{i.i.d.}\mathscr{N}\left(0,\,1\right)$ and $T_{1}^{0}=T\lambda_{1}^{0}$ with $\lambda_{1}^{0}=0.8$. Under this specification there is a break in the coefficient associated to the predictor $x_{t-1}$. Thus, there is a forecast failure and the test of giacomini/rossi:09 should reject $H_{0}$. From Table (ref) it appears that all versions of the classical HAC estimators of andrews:91 and newey/west:87, and KVB's fixed-$b$ lead to $t$-tests that have, essentially, zero power for all $\delta$. The only exception is Andrews' (1991) HAC estimator with prewhitening that shows some power but it is not monotonic. In contrast, the $t$-test based on the DK-HAC estimators have good power. The failure of existing LRV estimators cannot be attributed to the sample size because as we raise the sample size to 400 or 800, the tests still display no power {[}see casini_diss{]}.

The failure of the HAR tests based on the existing LRV estimators occurring in some of the data-generating mechanisms reported here can be simply reconciled with the fact that in such models the spectrum of $V_{t}$ is not constant. In other words, the autocovariance of $V_{t}$ depends not only on the lag order but also on $t$. Existing LRV estimators estimate an average of a time-varying spectrum. Because of this instability in the spectrum, they overestimate the extent of the dependence or variation in $V_{t}$. This is explained analytically in Casini et al. casini/perron_Low_Frequency_Contam_Nonstat:2020 who showed in a general setting that nonstationarity/misspecification alters the low frequency components of a time series making the latter appear as more persistent. Since traditional LRV estimators are a weighted sum of a large number of low frequency periodogram ordinates, these estimates turn to be inflated. Similarly, LRV estimators using long bandwidths (i.e., fixed-$b$) are weighted sum of a large number of sample autocovariances. Each sample autocovariance is biased upward so that the latter estimates are even more inflated than the classical HAC estimators that use a smaller number of sample autocovariances. This explains why KVB's fixed-$b$ HAR tests are subject to more power problems, even though the classical HAC estimators are also largely affected.

The introduction of the smoothing over time in the DK-HAC estimators avoids the low frequency contamination because observations belonging to different regimes are not mixed up when computing sample autocovariances. This guarantees good power properties also under nonstationarity/misspecification or under nonstationary alternative hypotheses (e.g., HAR tests for forecast evaluation discussed above). Casini et al. casini/perron_Low_Frequency_Contam_Nonstat:2020 reconciled this issue with some results in the unit root and long memory literature. Tests for a unit root are known to struggle to reject the unit root hypotheses if a process is second-order stationary (i.e., no unit root) but it is contaminated by breaks in the mean or trend {[}cf. perron:89 (perron:89, perron:90){]}. Similarly, a short memory sequence contaminated by structural breaks can approximate a long memory series in that the autocorrelation function has the same properties as that of a long memory series {[}cf. diebold/inoue:01, hillebrand:05, mccloskey/hill:2017 and mikosh/starica:04{]}.

Conclusions

Economic time series are highly nonstationary. Methods constructed under the assumption of stationarity might then have undesirable properties. This paper developed a theoretical framework for inference in settings where the data may be nonstationary. A new class of double kernel heteroskedasticity and autocorrelation consistent (DK-HAC) estimators was presented. In addition to the usual smoothing procedure over lagged autocovariances, the estimator applies smoothing over time. This is important in order to account flexibly for the variation over time of the structural properties of the economic time series. Optimality results under MSE criterion concerning bandwidths and kernels have been established. A data-dependent method based on the “plug-in” approach has been proposed. There are empirical relevant circumstances where HAR tests, either in linear regression models or other contexts, standardized by existing LRV estimators perform poorly. These may result in size distortions as well as significant power losses, even when the sample size is large. In contrast, when the proposed DK-HAC estimator is used the same HAR tests do not suffer from those issues. DK-HAC estimators lead to HAR tests that have competitive size control relative to fixed-$b$ HAR tests, when the latter work well, and have good power, irrespective of whether there is weak or strong dependence in the data.

\addcontentsline{toc}{section}{References}

\pagenumbering{arabic}