Econometrics and Statistics 27 (2023) 36–61

Contents lists available at ScienceDirect

Econometrics and Statistics
journal homepage: www.elsevier.com/locate/ecosta

Robust Covariance Matrix Estimation in Time Series: A Review
Masayuki Hirukawa
Faculty of Economics, Ryukoku University, 67 Tsukamoto-cho, Fukakusa, Fushimi-ku, Kyoto 612-8577, Japana r t i c l e

i n f oArticle history:
Received 26 February 2021
Revised 1 December 2021
Accepted 1 December 2021
Available online 16 December 2021
JEL classiﬁcation:
C12
C13
C22
C32

a b s t r a c t
In the analysis of economic, ﬁnancial and other time series, long-run variance estimatorsplay an important role in estimating model parameters more eﬃciently and drawing moreaccurate statistical inference on the parameters. A non-technical review of long-run variance estimation is provided. Both parametric and nonparametric estimators are discussed.
Kernel methods are dominant among all estimation procedures, and therefore recent developments in kernel-smoothed estimators and related inference are presented. The information given can help practitioners decide on a suitable long-run variance estimator.
© 2021 EcoSta Econometrics and Statistics. Published by Elsevier B.V. All rights reserved.

Keywords:
Bandwidth
Generalized method of moments
Heteroskedasticity and autocorrelation
robust inference
KernelLong-run variancePositive semi-deﬁnite

1. Introduction
Many applications of the analysis of economic, ﬁnancial and other time series involve estimation of a covariance matrix ofa zero-mean stationary vector process. An example of such a process is the cross-product of instruments and the regression
error. The process typically exhibits conditional heteroskedasticity and serial correlation. Heteroskedasticity and temporaldependence may be of unknown form.
To estimate model parameters more eﬃciently and draw more accurate statistical inference on the parameters in thiscircumstance, we must estimate the long-run variance (LRV) matrix of the process. Heteroskedasticity and autocorrelation
robust (HAR) covariance matrix estimators smoothed by kernel functions are frequently employed in LRV estimation. The
estimators of this class have long been called heteroskedasticity and autocorrelation consistent (HAC) covariance matrix
estimators due to the inﬂuential work by Andrews (1991). However, inconsistent kernel-smoothed estimators can yield welldeﬁned test statistics. Therefore, the abbreviation HAR is adopted throughout, unless otherwise noted.
This article provides an overview of LRV estimation. LRV estimation forms an integral part of the generalized method of
moments (GMM) estimation by Hansen (1982). With the development of GMM estimation, many excellent surveys on LRVor HAR estimation have been published. Comprehensive surveys on this subject include Den Haan and Levin (1997), Cushingand McGarvey (1999), and West (2008). Hall (2005, Section 3.5) focuses on LRV estimation within the framework of GMM

E-mail address: hirukawa@econ.ryukoku.ac.jp

https://doi.org/10.1016/j.ecosta.2021.12.001
2452-3062/© 2021 EcoSta Econometrics and Statistics. Published by Elsevier B.V. All rights reserved.

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

estimation, emphasizing that the inverse of an LRV estimator serves as the optimal weighting matrix of the eﬃcient GMM
criterion function. To the best of our knowledge, the most recent survey is given by Wang and Wu (2012). The primarytarget of that survey, however, is users of Stata, one of the most popular statistical packages among economists.
There are two objectives in this article. First, nearly a decade has passed since the latest survey was published. So, thisarticle not only presents a variety of LRV estimators, but it also discusses recent developments in HAR estimation and inference that have not been covered before. Second, the article emphasizes the role of LRV estimators as a tool for eﬃcientestimation of model parameters. On the one hand, development of LRV estimation has been paralleled by improvement in
HAR inference. Indeed, HAR inference procedures have been reﬁned since the seminal work on test statistics using inconsistent kernel-smoothed HAR estimators by Kiefer et al. (20 0 0) and Kiefer and Vogelsang (20 02, 20 05). On the other hand, little
is known about the implementation of HAR estimators in estimating parameters in GMM and cointegrating regressions efﬁciently. When employing a kernel HAR estimator for these estimation problems, we typically rely on a so-called automatic
bandwidth, which is an estimate of the bandwidth that minimizes the asymptotic mean squared error (AMSE) of the HAR
estimator. This bandwidth does not necessarily minimize the AMSE of the parameter estimator of interest. To the best of ourknowledge, the only available result is Wilhelm (2015), who derives the bandwidth that optimizes a quadratic loss based ona higher-order approximation to the eﬃcient GMM estimator. It is hoped that the Monte Carlo study in Section 6 can give
guidance for improving estimation-optimal bandwidth selection procedures one step further.
Furthermore, preference is given to practical rather than theoretical aspects. In particular, we pay attention to whether
each LRV estimator necessarily generates positive semi-deﬁnite (PSD) estimates. This property is highly desirable, because
the estimator will be used to construct the criterion function for eﬃcient GMM estimation and to compute standard errors
and test statistics.
Before proceeding, it is worth remarking on the relation of LRV estimation to spectral density estimation. The LRV matrix
of a zero-mean stationary vector process equals 2π times its spectral density matrix evaluated at frequency zero. Therefore,
LRV estimation dates back to pioneering works of spectral density estimation in statistical time series analysis such as
Jowett (1955), Hannan (1957) and Parzen (1957, 1961). There are also book-length treatments of spectral estimation available.
Examples include Hannan (1970), Anderson (1971), Brillinger (1975), and Priestley (1981). Kernel methods are dominant insuch earlier works. As GMM estimation became a standard tool in economics, spectral matrix estimation was revived in thename of HAR estimation. Indeed, kernel-smoothed HAR estimators studied by Newey and West (1987), Gallant (1987) andAndrews (1991), for example, inherit their form and their asymptotic properties from the literature on spectral estimation.
In addition, studies on HAR estimation and inference can be found among more than 200 articles in Econometrics andStatistics. Wenger and Leschinski (2021) improve HAR inference in the context of structural break testing, whereas Kawka(2020) investigates statistical properties of HAR estimators for locally stationary processes. Details will be described in
Sections 5.4 and 5.6.1, respectively.
The remainder of this article is organized as follows. Section 2 provides an overview of LRV estimation. Various applications of LRV estimates in estimation and testing problems are also documented. Procedures for LRV estimation may beclassiﬁed into two broad categories. Sections 3 focuses on parametric LRV estimators. Section 4 discusses nonparametric LRV
estimators with emphasis on sample autocovariance-based ones. While much of the material in Sections 3 and 4 appears
in previous surveys, Section 5 presents several recent developments in HAR estimation and inference that have not been
covered before. Section 6 reports the results of two Monte Carlo experiments. The simulation study focuses on appraising the inﬂuence of different HAR estimation procedures on parameter estimates of GMM and cointegrating regressions.
Section 7 summarizes and concludes the article.
√
The article adopts the following notational conventions: 0 p×q is the p × q zero matrix; i = −1;  is used to represent
the tensor (or Kronecker) product; I p is the p-dimensional identity matrix; L is the lag operator, i.e., L j Yt = Yt− j for j =1, 2, . . .; a ∨ b = max {a, b} and a ∧ b = min {a, b} for a, b ∈ R; · denotes the integer part; ∗ signiﬁes conjugate transposeof a complex-valued matrix A, i.e., A∗ = A ; K pp is the p2 × p2 commutation matrix that transforms vec (A ) into vec A ,p pi.e., K pp = j=1 k=1 e j ek  ek e j , where e j is the jth elementary p-vector; and χ 2 (ν ) and t (υ ) signify the chi-squared
distribution with ν degrees of freedom and the Student-t distribution with υ degrees of freedom, respectively. Lastly, limits
and orders of magnitude are taken as T → ∞ unless otherwise noted.
2. LRV Estimation: An Overview
2.1. Deﬁnition of the LRVLet gt be an  × 1 stationary vector process with E (gt ) = 0×1 . Throughout gt is assumed to have an  ×  spectral matrix.
Denote the jth-order autocovariance of gt for j = 0, 1, . . . as  j = E (gt gt− j ), where − j =  j . The  ×  LRV matrix of theprocess gt is in the form of

=

∞

j=−∞

 j = 0 +

∞

j=1

j +

∞


 j = 0 +  +  .

(1)

j=1

Most of the LRV estimators in this article can be obtained in the time domain. However, some are computed in the frequency
domain (see, e.g., Sections 5.1.1 and 5.6.2). For this purpose, the spectral density matrix of gt at frequency ω ∈ (−π , π ) is37

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

also considered. It is given by

fgg (ω ) =∞1 
 j e−i jω .
2π

(2)

j=−∞

Observe that  = 2π fgg (0 ) holds.
2.2. Roles of the LRV in the GMM Framework
There are two major reasons why we estimate the LRV matrix (1). These are (i) eﬃcient estimation and (ii) statisticalinference. To clarify the reasons, we further specify the process gt as gt = gt (θ0 ) = g(wt ; θ0 ) for a stationary vector process
wt and a known function g(·; · ) up to the parameter of interest θ0 ∈  ⊆ R p ( p ≤  ). Also assume thatE {gt (θ )} = 0×1

(3)

holds if and only if θ = θ 0 .
As a concrete example, consider a system of r equations

yt = Xt θ0 + ut

(4)

for the r × 1 vector of dependent variables yt , the p × r matrix of regressors Xt and the r × 1 vector of regression errors
ut . If Xt is endogenous and the  × r matrix of instruments Zt is chosen to estimate θ0 , then the orthogonality condition(3) becomes E (Zt ut ) = E Zt yt − Xt θ0= 0×1 .

(5)

Alternatively, as in Hansen and Singleton (1982), the h × 1 vector of instruments Rt may be orthogonal to each ele-ment of ut . In this scenario, we simply put Zt = Rt  Ir and  = hr. Observe that wt = yt , vec (Xt ) , vec (Zt )









and



wt = yt , vec (Xt ) , vec (Rt )in the ﬁrst and second cases, respectively, and that the moment function gt = Zt yt − Xt θ0
in both cases.
T
Given T observations {wt }t=1
, the orthogonality condition (3) motivates us to estimate θ0 by GMM so that θ˜ =Targ minθ ∈ G(θ ) WG(θ ), where G(θ ) = (1/T ) t=1
gt (θ ) is the sample counterpart of E {gt (θ )}, and W is some  ×  symmetric PSD weighting matrix. Now the reason (i) becomes clear. In eﬃcient GMM estimation, the inverse of a consistentˆ . Then, for the eﬃcient GMM estimator
estimator of (1) is chosen as W. Denote the estimator as 

ˆ −1 G(θ ),θˆ = arg min G(θ ) θ ∈

it holds that

√ 



(6)









−1
d
T θˆ − θ 0 → N 0 p×1 , V = N 0 p×1 , D −1 D, where D = E ∂ gt (θ )/∂ θθ = θ0.

ˆ −1 Dˆ = Dˆ ˆThe reason (ii) appears at this stage. The asymptotic variance V can be consistently estimated as V−1,

ˆ = ∂ G(θ )/∂ θˆ can be used to compute standard
where Dis the sample counterpart of D. The variance estimator Vθ =θˆerrors for θˆ and to construct test statistics for θ0 . Another example of inference is the J-test of overidentifying restrictions,

where J = T G θˆd
ˆ −1 G θˆ →

χ 2 ( − p) for  > p, under the null hypothesis that (3) holds true.

2.3. Applications of LRV Estimators
Other than GMM, there are many examples of applying LRV estimators for eﬃcient estimation and inference. In reference
to eﬃcient estimation, LRV estimators are used for eﬃciency gain via generalized least squares-type estimation methods
applied to seemingly unrelated regression models (e.g., Creel and Farell, 1996) and two-pass cross-sectional regressions ofasset returns (e.g., Ahn et al., 2012). Moreover, in some class of cointegrating regression models (e.g., Phillips and Hansen,
1990; Phillips, 1991; Park, 1992), an estimate of the one-sided LRV matrix = 0 +  = 0 +

∞


 j,

(7)

j=1

as well as that of the two-sided one (1), is employed to correct the so-called second-order effect that arises when the
regressor is endogenous.
LRV estimators are also applied for a variety of testing problems other than inference on parameters in regression andmoment-condition models. Examples include unit-root and stationarity tests (e.g., Phillips, 1987; Phillips and Perron, 1988;
Kwiatkowski et al., 1992), cointegration tests (e.g., Phillips and Ouliaris, 1990), structural break tests (e.g., Andrews, 1993),
and forecast comparison tests (e.g., Diebold and Mariano, 1995; West, 1996; McCracken, 20 0 0; Rossi and Inoue, 2012), toname a few.
38

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

2.4. Classiﬁcations of LRV Estimation Procedures
Our review of LRV estimation starts from the next section. Particular attention is paid to whether the LRV estimatorin discussion necessarily generates PSD estimates. This property is highly desirable, because the estimator will be used toconstruct the optimal weighting matrix for the eﬃcient GMM criterion function and to compute standard errors and test
statistics.
LRV estimation procedures can be categorized into two broad classes. The ﬁrst class is parametric procedures coveredin Section 3. In the parametric procedures, a linear time series model such as a vector autoregressive (VAR) or a vector
moving-average (VMA) model is ﬁtted to the vector process gt . Then, (2π times) the spectral density matrix at frequency
zero implied by the model is computed. Parametric LRV estimators become PSD by construction.
The second class is nonparametric procedures. The most popular estimator is in the form of a weighted sum of sample
autocovariances, where the weights are determined by a kernel and bandwidth. Section 4 focuses on this class of estimatorsand refers to a few kernels that necessarily generate PSD estimates.

3. Parametric LRV Estimation
3.1. VMA-Based Estimators
This section begins with the parametric LRV estimator implied by specifying gt as a VMA model of known order. Finiteorder serial dependence typically arises when gt is implied by a multi-period ahead predictive regression or comes from an
Euler equation for a rational expectations model.
Suppose that gt is known to obey a VMA model of order n a priori, as is the case with an (n + 1 )-period ahead predictive
regression.
it can be expressed as gt = εt + 1 εt−1 + · · · + n εt−n = (L )εt , where εt is the  × 1 innovation in gt with
 Then,ε = E εt εt , and (L ) = I + 1 L + · · · + n Ln is the lag polynomial with 1 , . . . , n being  ×  coeﬃcient matrices. It
follows that the LRV matrix of gt becomes  = (1 )ε (1 ) , and a consistent estimator of  can be obtained, in principle,
by replacing 1 , . . . , n with their consistent estimates. This estimator is PSD by construction and T 1/2 -consistent. However,it has not been applied in empirical works, because computational cost of estimating VMA coeﬃcient matrices 1 , . . . , n
is yet to be inexpensive. The estimation suffers from slow convergence of likelihood functions even for small-dimensional
systems.
Instead of directly estimating VMA coeﬃcients, Cumby et al. (1983) and Eichenbaum et al. (1988) ﬁrst estimate a ﬁniteorder (approximate) VAR model for gt and then invert the model to recover the VMA coeﬃcients. While the resulting LRV
estimate is PSD, its convergence rate may be no longer the parametric one; actually, the rate is determined by the expansion
rate of the lag order of the VAR model.
Hodrick (1992) and West (1997) take a different estimation strategy. Their estimators are PSD and T 1/2 -consistent but
can be computed with less computational burden. Suppose that the regression error ut in (4) is known to obey a VMA(n)model a priori. The focus here is on the estimator proposed by West (1997), because Hodrick’s (1992) estimator may beviewed as a special case of West’s (1997) estimator in which ut is a scalar and all MA coeﬃcients are set equal to unity.
Since r = dim (ut ) is usually (much) smaller than  = dim (gt ) (e.g., 1 = r ≤ p ≤  in a single-equation case), it is easier toˆ t to obtain u
ˆ t = εˆt +estimate the VMA(n) model for ut than gt . Then, West (1997) ﬁts a VMA(n) process to the residuals uψˆ 1 εˆt−1 + · · · + ψˆ n εˆt−n , where the εˆt are ﬁtted innovations with εˆt ≡ 0r×1 for t ≤ 0, and ψˆ1 , . . . , ψˆ n are consistent estimatesˆ t+n = Zt + Zt+1 ψˆ 1 + · · · + Zt+n ψˆ n εˆt forof r × r coeﬃcient matrices. Given the instrument matrix Zt deﬁned in (5), let d
t = 1, . . . , T − n. Then, the VMA(n)-based LRV estimator by West (1997) isˆ V MA =T −n

1 ˆ ˆ
dt+n dt+n .
T −n
t=1

3.2. VAR-Based Estimator
VAR-based LRV estimators are also available. A beneﬁt of VAR models is low computational cost, because coeﬃcients canbe estimated by ordinary least squares (OLS). Den Haan and Levin (1996, 1997, 20 0 0) propose to ﬁt an inﬁnite-order VAR
model to gt and then compute the LRV estimator as 2π times the spectral density matrix at frequency zero implied by thismodel. The resulting estimator is called the VARHAC estimator because of this particular procedure.
Before proceeding, a few remarks are in order. First, the idea of approximating univariate time series models by inﬁniteorder AR processes dates back to Berk (1974). In this respect, the VARHAC estimator can be viewed as an application of themultivariate extension of Berk’s (1974) method. Second, consistency of the VAR process is essential for consistency of the
VARHAC estimator. This is a sharp contrast to VAR-based prewhitening to be discussed in Section 4.2.3. Therefore, the lagorder must be chosen in a data-driven manner as documented below.
39

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

The VARHAC estimator is implemented as follows. Suppose that gt can be speciﬁed as a possibly inﬁnite-order VAR
model

gt =

∞


 j gt− j + t

(8)

j=1





with  = E t t , where the  j are  ×  coeﬃcient matrices. The VAR(∞) representation is possible if gt obeys an invertible VMA or VARMA model, as well as a ﬁnite-order VAR model.
To approximate (8) by a ﬁnite-order VAR model, let M = MT be the maximum lag order. It is assumed that MT → ∞.
Coeﬃcient estimation and order selection can be done in an element-by-element manner. To see this, suppose that an
AR(m) model is ﬁtted to ga,t , the ath element of gt , for m = 1, . . . , M and a = 1, . . . , . The ﬁtted model can be expressed as

ga,t =m φˆ a,b, j (m )gb,t− j + ˆa,t (m ), t = M + 1, . . . , T ,(9)

j=1 b=1

where the φˆ a,b, j (m ) are OLS estimates of coeﬃcients and ˆa,t (m ) is the residual. In this model, the number of lagged dependent variables ga,t− j and that of other predetermined elements gb,t− j for b = a are assumed to be the same. Alternatively, it
is possible to differ the former from the latter; see Den Haan and Levin (20 0 0) for more details. The optimal lag order forˆa = m
ˆ aT ∈ {1, . . . , M}, is determined as the minimizer of either Akaike’s (1973) information criterion (AIC)
ga,t , denoted as m



T
1  2
ˆa,t (m )T

AICa (m ) = log

+

t=M+1

2
m,
T

(10)

or Schwarz’ (1978) Bayesian information criterion (BIC)



T
1  2
ˆa,t (m )T

BICa (m ) = log

t=M+1

+

log T
m.
T

(11)

 

ˆj m
ˆ =m
ˆ T = max1≤a≤ m
ˆ a . Also let ˆ be the  ×  estimated coefDenote the maximum optimal lag order as m
 
 
 
ˆ
ˆ
ˆ , where the (a, b) element of  j m
ˆ is φa,b, j mˆ with φˆ a,b, j m
ˆ ≡ 0 for j > mˆ a.
ﬁcient matrix for lag j = 1, . . . , m
  j
 
mˆ
ˆ mˆ (L ) is given by ˆ mˆ (L ) = I −ˆ ˆ L . For the vector of residuals ˆt m
ˆ =Then, the estimated lag polynomial 
j=1  j m











ˆ
ˆ 1 , . . . , ˆ,t mˆ
ˆ1,t m, t = M + 1, . . . , T , the instantaneous variance of the innovation can be estimated by 
ˆ) =
ˆ (m
   −1 Tˆ ˆt mˆ . In the end, the VARHAC estimator is deﬁned as(T − M )t=M+1 ˆt m−1ˆ VAR = ˆ mˆ (1 )−1 ˆˆ ˆ (1 ).
ˆ ) mˆ (m(12)

The VARHAC
 estimator is PSD by construction. In addition, Den Haan and Levin (1996, Theorem 3) demonstrate that if
ˆ T is determined via BIC (11) and some other regularity conditions hold, then the
MT = O T 1/3 and the optimal lag order mestimator is nearly T 1/2 -consistent.
3.3. Discussion
Parametric LRV estimators have a few advantages over nonparametric, kernel-smoothed ones. The former necessarily
generates PSD estimates and attains (nearly) T 1/2 -consistency. Simulation results by West (1997) and Den Haan and Levin(1997) indicate that test statistics using their parametric LRV estimators often outperform those using nonparametric kernel
estimators. On the other hand, whether a kernel-smoothed LRV estimator is PSD depends on the choice of a kernel, and its
convergence rate is usually slower than the parametric one.
Nevertheless, parametric LRV estimators are less popular than nonparametric ones in practice. This can be attributed
to diﬃculty in estimating model parameters precisely. First, estimating a VMA model is still computationally challenging.
Therefore, a kernel-smoothed estimator with a suitable choice of the bandwidth is preferred even when an economic orother theory suggests that gt has ﬁnite-order dependence. Second, estimating many parameters in the VAR model is aserious concern for VARHAC. In (9) we use (T − M ) observations to estimate m parameters. When  is large relative toT (e.g., for the case of many instruments in (5)), we may resort to regularized VAR estimation (e.g., Basu and Michailidis,2015). However, VARHAC in this direction is yet to be developed, to the best of our knowledge. Moreover, the process gt
computed from economic and ﬁnancial data often exhibits high persistence. Then, OLS estimates of VAR coeﬃcients haveconsiderable ﬁnite-sample bias, which is likely to deteriorate ﬁnite-sample properties of the VARHAC estimator, as suggested
by Sul et al. (2005).
40

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

4. Nonparametric LRV Estimation
4.1. Truncated Estimator
Among all nonparametric LRV estimators, sample autocovariance-based ones are dominant. Therefore, below we focus
exclusively on this class. Suppose that an economic or other theory suggests that gt has zero autocovariances after lag n,

as seen in Section 3.1. The LRV matrix reduces to  = nj=−n  j in this case. To estimate , deﬁne the jth-order sam
ˆ j = (1/T ) (T + j )∧T gt g . A natural estimator of  is its sample
ple autocovariance of gt for j = 0, ±1, . . . , ±(T − 1 ) as 
t=1∨ (1+ j )analogˆ TR =


n


t− j

ˆ j.


j=−n

This is called the truncated estimator, and it is studied in the early literature on GMM including Hansen (1982), Hansen andSingleton (1982), and White and Domowitz (1984). Hansen and Hodrick (1980) consider a similar but more restrictive LRV
estimator, but its nature is the same.
The truncated estimator is consistent for more general cases. White and Domowitz (1984, Theorem 3.5) establish that
ˆ T R becomes consistent if the lag truncation point n = nT satisﬁes 1/nT + nT /T 1/3 → 0 andeven when  j = 0× for all j, some other regularity conditions hold.
A disadvantage of the truncated estimator is that it does not necessarily generate PSD estimates. West (1997) andˆ T R quite frequently yields negative variance estimates in the presHirukawa (2010) report in their simulation studies that 
ence of negative serial dependence. Therefore, this estimator is not much used in practice.
4.2. Weighted Autocovariance Estimators
4.2.1. Deﬁnition and Kernel Choicesˆ = T −1ˆThe most popular LRV estimators take the form of a weighted sum of sample autocovariances 
j=−(T −1 ) w j  j . Following the literature on spectral density estimation, Andrews (1991) explores a general framework for the analysis on thisclass of estimators. In the framework, the weights w j are determined by a kernel k(· ) and a sequence of bandwidth ST ∈ R+
so that w j = k( j/ST ) and
T −1


ˆJ =


j=−(T −1 )



kJ

j
ST



ˆ j,


(13)

where the subscript J signiﬁes that the kernel kJ (· ) is chosen. Such LRV estimators have long been called HAC estimators
due to the pioneered work by Andrews (1991). However, the estimator (13) using the bandwidth ST = T is inconsistent but
can construct well-deﬁned test statistics, as in Section 5.3. Therefore, we refer to what was called HAC estimators as HAR
estimators, whenever there is no confusion.
Some conditions on a kernel k(· ) and bandwidth ST are required for consistency of HAR estimators. A typical set of∞
regularity conditions on k(· ) are summarized as Assumption 1 below. The ﬁnal condition 0 k̄(x )dx < ∞ ensures that certain
Riemann-type sums deﬁned in terms of the kernel converge to their integral representation counterparts; see Jansson (2002,
p.1451) for discussion.
Assumption 1. k : R → [−1, 1], k(0 ) = 1, k(x ) = k(−x ), for all x ∈ R, k(· ) is continuous at 0 and almost everywhere, and∞
0 k̄ (x )dx < ∞ where k̄ (x ) = supy≥x |k (y )|.
Quite a few kernels satisfy Assumption 1, and the bandwidth selection methods in the next section are also designed forthese kernels. From the practical perspective, however, preference is given to the following most commonly used kernels.



Truncated:

kT R ( x ) =Bartlett:

kBT (x ) =Parzen:Quadratic Spectral:

1
0

for |x| ≤ 1otherwise

1 − |x|0

for |x| ≤ 1otherwise⎧⎨1 − 6x2 + 6|x|3 for |x| ≤ 1/2
kPZ (x ) = 2(1 − |x| )3
for 1/2 < |x| ≤ 1⎩0otherwise 
5)kQS (x ) = 12π252 x2 sin6(π6πx/x/− cos 6π5 x
5The truncated estimator discussed in the previous section can be viewed as the HAR estimator in which the truncatedkernel is used and the bandwidth ST is set equal to n. The Bartlett, Parzen and quadratic spectral (QS) kernels are investigated by Newey and West (1987), Gallant (1987, p.533) and Andrews (1991), respectively.
41

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

Table 1
Characteristic Numbers of Kernels Widely Applied in HAR Estimation
Kernelqkq

Bartlett
Parzen
QS

1
2
2

1
6
18π 2 /125

∞

−∞ k

2

(x )dx2/3
151/2801

∞−∞|x|2q k2 (x )dx

1/15
929/295680∞

Whether a kernel necessarily generates PSD estimates is of particular importance and interest. To answer this question,∞deﬁne the spectral window generator of a kernel k(· ) (Andrews, 1991) as K (λ ) = (2π )−1 −∞ k(x )e−ixλ dx. The kernel k(· )necessarily generates PSD estimates if K (λ ) ≥ 0 for all λ ∈ R. The Bartlett, Parzen and QS kernels satisfy this condition, but
the truncated kernel does not. Moreover, Andrews (1991) demonstrate that the QS kernel is optimal in the sense that it
attains the smallest AMSE among all HAR estimators smoothed by kernels yielding PSD estimates. The spectral window
generator of this kernel reduces to the Epanechnikov kernel in the frequency domain, which is known to be optimal among
all nonnegative kernels. Optimality of the QS kernel can be understood analogously.
Now we turn to the regularity condition on the bandwidth ST . For consistency of HAR estimators, let ST satisfy 1/ST +ST /T → 0. This assumption, jointly with examples of the kernels satisfying Assumption 1, implies the following intuitions.
On the one hand, it follows from ST /T → 0 that k( j/ST ) = 0 (or k( j/ST ) → 0 at least) for | j| near T − 1. As a consequence,
sample autocovariances with long lags, which tend to be estimated imprecisely, are downweighted. On the other hand,
ST → ∞ implies that k( j/ST ) → k(0 ) = 1 for each j, which leads to consistency of a HAR estimator. The requirement for the
bandwidth to diverge also draws the important conclusion that we must set ST → ∞ even when gt is known to have at
most n non-zero autocovariances. Putting ST = n is allowed only if the truncated kernel is employed. However, this comesat the cost of non-PSD estimates.
4.2.2. Data-Driven Bandwidth Selection
How to choose the bandwidth ST for a given kernel k(· ) is an important practical problem. Finite-sample performancesof a HAR estimator is largely inﬂuenced by the bandwidth, and the choice of a kernel is thought to be of second importancein practice. This section addresses the issue of selecting ST in a data-driven manner.
To begin with, we refer to measures of smoothness for a kernel and spectral density. A kernel k(· ) is said to have thecharacteristic exponent q (Parzen, 1957) if it has the following properties:
kr = lim

x→0

1 − k (x )
|x|r

 =0
∈ ( 0, ∞ )=∞for r < qfor r = q ,for r > q

where kr and kq are called the rth generalized derivative of the kernel at the origin and the characteristic coeﬃcient (Parzen,
1961), respectively. In addition, the rth generalized derivative of the spectral density (2) at frequency ω ∈ (−π , π ) is deﬁned as(r )r−i jω . Since the truncated kernel has k ≡ 0 for any r > 0 and thus q = ∞, it does not ﬁt withfgg(ω ) = (2π )−1 ∞rj=−∞ | j |  j e
the AMSE-optimal framework below. Furthermore, Parzen (1961, p.179) shows that a kernel that necessarily generates PSD
estimates must satisfy q ≤ 2, and remaining three kernels indeed have q = 1 (Bartlett) or 2 (Parzen, QS). Table 1 presents
values of q and kq , as well as some other characteristic numbers, of these kernels.
We adopt the AMSE as the loss function of a HAR estimator. The AMSE approximates the mean squared error (MSE) of theestimator by AMSE = ABias2 + AV ar, where ABias and AV ar denote dominant terms of the bias and variance of the estimator,
respectively. Andrews (1991) demonstrates that if a HAR estimator is built on a kernel having the characteristic exponent q,−qthen under suﬃcient smoothness of the spectral density at frequency zero, ABias = O(ST ) and AV ar = O(ST /T ). Therefore,−2q+ ST /T ), and thus the estimator attains the best possible convergence rate of T q/(2q+1 ) underpˆ BT →
the bandwidth ST = O(T 1/(2q+1 ) ) that balances ABias2 and AV ar. This means that  at rate T 1/3 if ST = O(T 1/3 ) and
AMSE = ABias2 + AV ar = O(STp

ˆ PZ , ˆ QS →  at rate T 2/5 if ST = O(T 1/5 ).
that 
However, the above result is merely a guidance, and a bandwidth formula of the form ST = CT 1/(2q+1 ) for some constantC > 0 is required for practical purposes. Andrews (1991), Newey and West (1994) and Hirukawa (2010) provide the socalled automatic bandwidth as an estimate of the minimizer of the AMSE of a HAR estimator using a general class ofkernels. There are important differences in the deﬁnition of the MSE of a HAR estimator and the approach to estimating an
unknown quantity in the AMSE-optimal bandwidth. Below we describe the three bandwidth choice rules by illustrating thedifferences.
ˆ as E {  wa (ˆ aa − aa )2 }, where the wa ≥ 0 are
Andrews (1991): Andrews (1991) deﬁnes the MSE of a HAR estimator a=1ˆˆweights, and aa and aa are the ath diagonal elements of  and , respectively. The deﬁnition yields the AMSE-optimal
bandwidth formula as



ST† =

qk2q α (q )

∞

−∞ k

2

(x )dx1/(2q+1)

T 1/(2q+1 ) ,

(14)
42

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

where α (q ) =a=1 wa2(q )fgg,aa/2
a=1 wa f gg,aa , and

(q )

(q )

fgg,aa and fgg,aa are the ath diagonal elements of fgg (0 ) and fgg (0 ),respectively.
(q )To make (14) fully operational, we must prespecify the weights wa and replace the unknown fgg,aa and fgg,aa with theirproxies. For the former, Andrews (1991) set all weights corresponding to slopes equal to unity and the one corresponding
to an intercept equal to zero for regression models. For the latter, an AR(1) model as a reference is ﬁtted to each elementof gt . In the end, using the numbers in Table 1 yields Andrews’ (1991) data-driven (or automatic) bandwidth as⎧⎨1.1447{αˆ (1 )}1/3 T 1/3†
SˆT = 2.6614{αˆ (2 )}1/5 T 1/5⎩
1/5
1.3221{αˆ (2 )} T 1/5

where





a=1 waαˆ (1 ) =for Bartlett
for Parzen ,for QS

4ρˆa2 σˆ a4(1−ρˆa ) (1+ρˆa )


6



a=1 waa=1 wa

2

σˆ a4
4
1
−
( ρˆa )

, αˆ (2 ) =

a=1 wa4ρˆa2 σˆ a4(1−ρˆa )

8



σˆ a4
4
1
−
( ρˆa )

,and (ρˆa , σˆ a ) are OLS estimates of the coeﬃcient and the instantaneous variance of the innovation in the AR(1) model ﬁtted
to the ath element of gt , respectively.
Andrews’ (1991) approach is analogous to the “rule of thumb” for probability density estimation by Silverman (1986,
Section 3.4.2). A potential problem is that the data-driven bandwidth is not consistent for the AMSE-optimal bandwidth
unless the AR(1) model is elementwise correct speciﬁcation of the process gt . Hence, this approach may perform poorly
when gt is not well approximated by an AR(1) model.
ˆ − )w}2 as an
Newey and West (1994). For an  × 1 vector of weights w, Newey and West (1994) consider E {w (ˆ . This deﬁnition leads to the AMSE-optimal bandwidth formulaalternative deﬁnition of the MSE of a HAR estimator 



ST† =



2 1/(2q+1)

qk2q R(q )∞2−∞ k (x )dx

where

R (q ) =



(q )
w fgg
w

w fgg w

T 1/(2q+1 ) ,

∞

(15)

| j|q γ j s(q)= (0 ) ,s
j=−∞ γ j

j=−∞

∞

=



γ j = E ht ht− j is the jth autocovariance of the process ht = w gt , and s(r ) =∞

rj=−∞ | j | γ j . Observe that (14) and (15) co-

incide if gt is a scalar process. Because the quantity R(q ) is the ratio of a density derivative to the density itself, Hirukawa(2010) refers to it as the normalized curvature.
When implementing (15), Newey and West (1994) estimate R(q ) nonparametrically using the truncated kernel to avoid
possible misspeciﬁcation of the process ht . As in Andrews (1991), Newey and West (1994) also set all weights in w corresponding to slopes equal to unity and the one corresponding to an intercept equal to zero for regression models. The
truncated estimator for R(q ) using the bandwidth n isqRˆT( R) =where γˆ j =n| j|q γˆ j,ˆjj=−n γ(16)

t=1+ j ht ht− j /T

is the jth-order sample autocovariance of ht for j = 0, 1, . . . , T − 1 with γˆ− j = γˆ j . The problem

j=−n

n

T

†

is that it is hard to ﬁnd the optimal bandwidth n† = nT for this estimator. Rather, Newey and West (1994) determine n in
an ad hoc manner. The ﬁnal form of their data-driven bandwidth isSˆT† =⎧
 2 1/3
1
⎪
⎪
1
.
1447
RˆT( R)
T 1/3⎪⎪⎨
  1/5
2.6614

2
Rˆ( )

2

T 1/5

TR
⎪
⎪
 2 1/5⎪⎪
⎩1.3221 Rˆ(2)
T 1/5
TR

where

for Bartlett
for Parzen ,for QS⎧  2/9    2/9 
T
T
4 100
, 12 100⎪⎪⎪⎨       
4/25
4/25
T
T
n=
4 100
, 12 100⎪   2/25 ⎪⎪⎩  T 2/25
T
3 100, 4 100

for Bartlett
for Parzen .
for QS43

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

Hirukawa (2010): Yet another approach by Hirukawa (2010) is an analog to the bandwidth choice rule for probability
density estimation by Sheather and Jones (1991). It is built on the AMSE-optimal bandwidth (15) by Newey and West (1994).
Two remarkable differences are: (i) that the normalized curvature R(q ) is estimated nonparametrically using the same kernel†k(· ) but a different bandwidth bT ; and (ii) that the AMSE-optimal bandwidth bT for the normalized curvature estimatorT −1

j=−(T −1 ) k

Rˆ(q ) (bT ) =T −1( j/bT )| j|q γˆ jj=−(T −1 ) k( j/bT )γˆ j

is derived. It is



b†T =

qk2q β 2 (q )∞
(2q + 1 ) −∞ |x|2q k2 (x )dx1/(4q+1)

T 1/(4q+1 ) ,

(17)

2

where β (q ) = s(q ) /s(0 ) − s(2q ) /s(0 ) .
Like Newey and West (1994), it is possible to conduct sequential estimation of the normalized curvature (ﬁrst-stage)††
and the LRV (second-stage) by using estimates of bT and ST . Instead, in the spirit of Sheather and Jones (1991), Hirukawa†(2010) adopts the solve-the-equation plug-in (SEPI) approach as the algorithm for estimating ST . It starts from constructing

a system of two equations with two unknowns


T =

†

†

ST , bT . Solving (15) for T gives(x )dx  † 2q+1.
 2 ST
2

∞
2
−∞ k

qkq R(q )Substituting this into (17) and invoking that R(q ) = s(q ) /s(0 ) , we haveb†T = b†T

 †
ST





∞
δ 2 (q ) −∞k2 (x )dx=∞
(2q + 1 ) −∞ |x|2q k2 (x )dx

1/(4q+1) † (2q+1)/(4q+1)
ST



,

(18)

where δ (q ) = s(q ) /s(0 ) − s(2q ) /s(q ) . Finally, replacing R(q ) in (15) with Rˆ(q ) bT


ST† =



†

 

= Rˆ(q ) bT ST††

yields  2 1/(2q+1)

qk2q Rˆ(q ) b†T ST†∞2−∞ k (x )dx

T 1/(2q+1 ) .

(19)†

The SEPI approach numerically solves the system of equations (18) and (19) for ST . Notice that this does not work for the∞
HAR estimator using the QS kernel, because Table 1 indicates that −∞ |x|2q k2QS (x )dx = ∞.
To operationalize this algorithm, we replace the unknown quantity δ (q ) with its proxy from ﬁtting an AR(1) model toht , as in Andrews (1991). Sheather and Jones (1991) justify ﬁtting a parametric model at this stage by arguing that doing
so is less crucial than ﬁtting it directly to R(q ) . The proxies of δ (q ) for q = 1, 2 using the OLS estimate of the autoregressive





2
†
†
parameter ρˆ are δˆ (1 ) = ρˆ 2 + 1 / ρˆ 2 − 1 and δˆ (2 ) = − ρˆ 2 + 8ρˆ + 1 / ρˆ − 1 . Let Sˆ be the estimate of S . Also denote



T

T

†
†
†
†
the estimate of bT as bˆ T = bˆ T SˆT . Then, the ﬁnal form of the system is given by

⎧

  
⎪
⎨ Sˆ† = 1.1447 Rˆ(1) bˆ † Sˆ†
T

BT

T

2

1 / 3TT 1/3⎪⎩ bˆ † Sˆ†  = 1.2723 δˆ 2 (1 ) 1/5 Sˆ† 3/5
T
T
T
⎧

   2 1 / 5⎪⎨ Sˆ† = 2.6614 Rˆ(2) bˆ † Sˆ†
T 1/5
T

PZ

T

T

⎪
⎩ bˆ † Sˆ†  = 1.4812 δˆ 2 (2 ) 1/9 Sˆ† 5/9
TT

for Bartlett.
for Parzen

T

By construction, the SEPI approach inherits properties of both Andrews’ (1991) parametric and Newey and West’s(1994) nonparametric approaches. While ﬁtting a parametric model to the unknown process ht is required for implementation, the rule itself is built on a nonparametric estimator of the normalized curvature, which limits the inﬂuence of theparametric reference. Nice ﬁnite-sample properties of HAR estimates combined with the SEPI approach are reported in
Hirukawa (2010, 2011). Disadvantages of the approach are excluding the QS kernel and solving the highly nonlinear systemnumerically. For the latter, there may be multiple roots. In such cases, Hirukawa (2010) recommends taking the largest root†as SˆT .
44

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

4.2.3. Prewhitening
The process gt computed from economic and ﬁnancial data is often highly persistent. Essentially, its spectral density has
a sharp peak at frequency zero, and as a result, covariance estimates from HAR estimators tend to be imprecise. It has long
been known in the literature on statistical time series analysis (e.g., Press and Tukey, 1956; Blackman and Tukey, 1958) thatthe spectral density at frequency zero is easier to estimate if it is ﬂat in the vicinity of the origin. This motivates us to ﬁlter
out the original process gt via a simple time series model in order to estimate the spectral density of the ﬁltered process
at frequency zero more easily. Prewhitening studied by Andrews and Monahan (1992) and Lee and Phillips (1994) is built
on the above idea. The procedure of recovering the spectral density of the original process from the one of the ﬁlteredprocess is called recoloring. Several simulation results justify prewhitening at the beginning and recoloring at the end of
HAR estimation by indicating better ﬁnite-sample properties.
Our focus is on VAR-based prewhitening by Andrews and Monahan (1992), which is most popular in practice. Below the
HAR estimation with prewhitening and recoloring is described in a step-by-step manner.
1. Prewhitening: Fit a low-order VAR model to the original process gt , and obtain OLS estimates of autoregressive coˆ gt−1 + ˆt , t = 2, . . . , T , A
ˆ is the OLS
eﬃcient matrices and the VAR residuals. In the example of a VAR(1) model gt = A
estimate of the coeﬃcient matrix and the ˆt are the VAR residuals.
2. HAR estimation: Compute the LRV estimate of the residual ˆt by HAR estimation in combination of a prespeciﬁed kernelkJ (· ) with one of its corresponding bandwidth choice methods discussed in Section 4.2.2. The resulting LRV estimateˆˆ.
using is denoted as J
3. Recoloring: Obtain the LRV estimate of gt as
−1 ˆ ˆ J −PW = I − A
ˆ
ˆ
ˆ


I − AJ−1.

There are a few remarks in order. First, the sole purpose of prewhitening is to ﬂatten the spectral density of the ﬁltered process in the neighborhood of frequency zero. Therefore, the time series model used for prewhitening does not need
to be correctly speciﬁed, unlike the VARHAC estimator in Section 3.2. Second, Lee and Phillips (1994) propose yet anotherprewhitening method based on an autoregressive moving-average (ARMA) model. Orders and coeﬃcient estimates of the
ARMA model are determined by the Hannan-Rissanen recursion (Hannan and Rissanen, 1982). This procedure is rarely applied in practice, because it is designed for a scalar process. Extending it to a vector process and ﬁtting a VARMA model
would incur additional computational cost in the MA part, as mentioned in Section 3.1.
4.3. Estimator without Smoothing
Robinson (1998) is concerned with using a single bandwidth for all elements of a kernel-smoothed HAR estimator (13).
So far it has been implicitly assumed that all elements of a HAR estimator are smoothed by a single bandwidth. This practice
looks unavoidable, because varying bandwidth numbers across elements destroys the property of generating PSD estimateseven when the Bartlett, Parzen or QS kernel is employed. Actually, the practice may be suboptimal in that the same, ﬁxed
amount of smoothing is imposed on all elements in the spectral matrix that are likely to have a variety of shapes. The
concern motivates Robinson (1998) to propose yet another LRV estimator that is based on sample autocovariances but freeof kernel smoothing.
Robinson’s (1998) LRV estimator is built on the process corresponding to (5) with Zt = Rt  Ir . More speciﬁcally, letthe  × 1 process be gt = ut  Rt , where ut is an r × 1 vector of regression errors with E (ut ) = 0r×1 , Rt is an h × 1 vec
ˆ u = (1/T ) (T + j )∧T ut uˆR =tor of instruments, and  = hr. Let the sample autocovariances of ut and Rt be 
and 
t− j
j
j
t=1∨ (1+ j )(T + j )∧T 
T(1/T ) t=1∨(1+ j ) Rt − R̄ Rt− j − R̄ , respectively, where R̄ =(1/T ) t=1 Rt . The LRV estimator is then deﬁned as

ˆR =


T −1


ˆu  
ˆ R.
jjj=−(T −1 )

ˆ R are that it is T 1/2 -consistent and that it necessarily generates PSD estimates.
Further advantages of 
A disadvantage of this estimator is that its consistency relies on the conditionE (ut  Rt ) ut− j  Rt− j









= E ut ut− j  E Rt Rt− j

for all j, which does not allow for heteroskedasticity in ut . Furthermore, all elements of ut and Rt must be stochastic, andthus the natural instrument “1” cannot enter as an element of Rt .
4.4. A Real Data Example
It is of practical importance and interest to see how different HAR estimators produce different LRV estimates. The test foruncovered interest parity (UIP) by Hansen and Hodrick (1980) is chosen as an example. The test is based on the regression45

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

Table 2
Results of Predictive Regressions and UIP Tests
Regressand
GBP:

βˆ0

βˆ1

βˆ2

0.00330.03000.3002

TR-HH
QS-A
BT-NW
PZ-H

(0.0058 )(0.0055 )(0.0051 )(0.0052 )(0.0920 )(0.0724 )(0.0840 )(0.0859 )(0.2278 )(0.2385 )(0.2055 )(0.2098 )

TR-HH
QS-A
BT-NW
PZ-H

(0.0051 )(0.0047 )(0.0041 )(0.0045 )(0.0971 )(0.0916 )(0.0827 )(0.0889 )(0.1797 )(0.1721 )(0.1549 )(0.1665 )0.01190.2990

TR-HH
QS-A
BT-NW
PZ-H

(0.0054 )(0.0052 )(0.0036 )(0.0051 )(0.1046 )(0.1080 )(0.0805 )(0.1046 )(0.2286 )(0.2105 )(0.1679 )(0.2228 )

TR-HH
QS-A
BT-NW
PZ-H

(0.0063 )(0.0052 )(0.0051 )(0.0053 )(0.1154 )(0.1069 )(0.1062 )(0.1079 )(0.3073 )(0.3302 )(0.3070 )(0.3132 )FFR:

0.0017−0.0025

DM:

JPY:

0.0074−0.0269−0.1214

βˆ3

βˆ5−0.1473(0.1714 )(0.1697 )(0.1567 )(0.1599 )−0.0344(0.1416 )(0.1270 )(0.1278 )(0.1341 )−0.1524(0.1695 )(0.1664 )(0.1320 )(0.1757 )−0.3210(0.2280 )(0.2182 )(0.2086 )(0.2118 )0.27600.2015Wald[ p-value]BW(0.0582 )(0.0545 )(0.0517 )(0.0531 )3.95344.45484.93564.7120[0.5561][0.4860][0.4238][0.4520]

2.00
9.61
3.94
5.47

(0.0596 )(0.0536 )(0.0500 )(0.0541 )5.86456.74288.5216
7.1868[0.3196][0.2405][0.1297][0.2071]

2.00
8.18
2.28
4.61

3.07114.23876.13523.6028[0.6890][0.5156][0.2933][0.6079]

2.00
9.47
1.30
11.90

48.260827.705435.720334.8245[0.0 0 0 0]
[0.0 0 0 0]
[0.0 0 0 0]
[0.0 0 0 0]

2.00
11.11
9.55
12.92

−0.0859−0.0283−0.0059(0.0774 )(0.0543 )(0.0546 )(0.0600 )0.5616(0.0950 )(0.1222 )(0.1071 )(0.1085 )Note: Numbers in parentheses, “Wald”, “p-value”, and “BW/LO” for each LRV estimator are standard errors, values of Wald statistics, their associated pvalues from χ 2 (5 ), and automatic bandwidths SˆT† (for the QS, Bartlett and Parzen estimators) or lag orders n (for the truncated estimator), respectively.

of the forecast error of the exchange rate of a currency c on lagged forecast errors of own and other currenciesccstc − ft−3
,3 = β0 +

4


 j

j
β cj st−3− ft−6+ tc ,,3

j=1

where stc andft,c 3 are logarithms of spot and 3-month-forward exchange rates of a currency c at time t, respectively. We

work with 352 monthly observations from January 1973 to December 1999 in the data set “fxeff_m.prn” available on
the webpage of Mark (2001). Four exchange rates (expressed in the US dollar), namely, pound sterling (GBP; c = 1), Frenchfranc (FFR; c = 2), Deutsche mark (DM; c = 3), and Japanese yen (JPY; c = 4), are examined.
The regression can be estimated by OLS, and then the UIP test for the null hypothesis H0c : β cj = 0 for j = 0, . . . , 4 isconducted. To compute standard errors of OLS estimates and Wald statistics for the testing of H0c , we consider the following
four HAR estimators: (i) the truncated estimator with 2 lags employed by Hansen and Hodrick (1980), reﬂecting that theerror term tc behaves like a MA(2) process [TR-HH]; (ii) the QS HAR estimator using Andrews’ (1991) automatic bandwidth[QS-A]; (iii) the Bartlett HAR estimator using Newey andWest’s (1994)
 automatic bandwidth, where the bandwidth for

Tthe truncated normalized curvature estimator (16) is n = 4 100

2 / 9[BT-NW]; and (iv) the Parzen HAR estimator using

Hirukawa’s (2010) SEPI bandwidth [PZ-H].
Table 2 presents OLS estimates, standard errors, Wald statistics with associated p-values from χ 2 (5 ), and automatic
bandwidths or lag orders. While it is not possible to judge which HAR estimator is most reliable, it is safe to say that theresults are largely inﬂuenced by the choice of the estimator.
4.5. Discussion
Practitioners’ preference has long been given to kernel HAR estimators, and the Bartlett and QS estimators are mostpopular among all such estimators. This reﬂects that the convergence rate of a LRV estimator is not well transmitted toits ﬁnite-sample properties. Although the truncated and VARHAC estimators have faster convergence rates, simulation results indicate that the former often and the latter sometimes perform more poorly than the slower QS estimator, which isoccasionally outperformed by the even slower (and thus theoretically inferior) Bartlett estimator.
However, test statistics using the Bartlett and QS estimators often exhibit unsatisfactory ﬁnite-sample performances. In
particular, substantial size distortions of test statistics in the presence of strong positive dependence are reported and extensively discussed in recent work (e.g., Kiefer and Vogelsang, 2005; Sul et al., 2005).
While the bandwidth determines ﬁnite-sample performance of a HAR estimator, there is uncertainty in the bandwidth
choice itself. Even if the bandwidth is chosen in a data-driven manner, practitioners must make some kind of arbitrary
decision (e.g., which rule given in Section 4.2.2 to choose and how to estimate the normalized curvature). Moreover, there
is no guarantee that the AMSE-optimal bandwidth is equally optimal for parameter estimation or statistical inference. From
this viewpoint, two object-oriented (i.e., test- and estimation-optimal) bandwidth choices, as well as several other topics,
will be discussed in the next section.
46

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

5. Recent Developments in HAR Estimation and Inference
While much material in the previous sections is included in past surveys, research on HAR estimation and inference
continues. This section is devoted to developments over the past two decades that are not covered before.
5.1. Kernel HAR Estimators with Faster Convergence Rates
A faster convergence rate in kernel HAR estimation is a theoretically appealing property. A HAR estimator having an
accelerated convergence rate is more valuable if it generates PSD estimates. Below we present two examples of such HAR
estimators.
5.1.1. Nonparametrically Prewhitened Estimator
The ﬁrst estimator, proposed by Xiao and Linton (20 02) and Hirukawa (20 06), attains the rate improvement via a multiplicative bias correction (MBC) of the spectral matrix estimator in the frequency domain. The MBC
 method, originallyproposed by Jones et al. (1995) for probability density estimation, is based on the identity f (· ) ≡ f˜(· ) f (· )/ f˜(· ) = f˜(· )α (· )for f˜(· ) > 0, where f (· ) and f˜(· ) are the true density and its initial estimator, respectively, and α (· ) can be recognized as thebias correction term. It is anticipated that replacing α (· ) with its consistent estimate can yield a less-biased estimator thanf˜(· ), and non-negativity of the bias-corrected estimator is guaranteed by construction. The MBC procedure is reminiscentto prewhitening when it is applied to spectral estimation. Then, Xiao and Linton (2002) call the MBC-HAR estimator the
nonparametrically prewhitened (NPW) estimator.
Because the MBC is originally designed for scalar quantities, transplanting it to spectral matrix estimation should be donewith care. Hirukawa (2006) proposes the following procedure of computing the NPW estimator. Throughout the kernel k (· )is assumed to have the characteristic exponent q = 2 and the spectral window generator K (· ) ≥ 0. For such K (· ) and thebandwidth ST , the amplitude window (Parzen, 1963) is given by KST (λ ) = ST ∞j=−∞ K {ST (λ + 2π j )} ≈ ST K (ST λ ). Then, thespectral matrix of gt evaluated at frequency ω ∈ (−π , π ) can be estimated as

2π
f˜gg (ω ) =
T





  

KST λ j − ω Igg λ j ,λ j ∈B ( ω )    ∗ 
Twhere Igg λ j = ζg λ j ζg λ j is the periodogram with ζg λ j = (2π T )−1/2 t=1
gt e−it λ j being the ﬁnite Fourier transform
of gt evaluated at fundamental frequencies λ j = 2π j/T for j = 0, ±1, . . . , ±T /2, and B(ω ) = λ j ω − π < λ j < ω + π .
1/2
Observe that f˜gg (ω ) is Hermitian and positive deﬁnite (PD). Hence, f˜gg (ω ) has eigenvalues λ1 , . . . , λ > 0, and thus f˜gg
(ω )1/2
is well-deﬁned by the unitary decomposition f˜gg(ω ) = U1/2 U∗ , where 1/2 = diag λ11/2 , . . . , λ1 /2

is the diagonal matrix∗

containing the square roots of the eigenvalues, and U is the unitary matrix, i.e., UU = I . In the end, the NPW estimator canbe deﬁned as

ˆ NPW = ˜ 1/2 α˜ (0 )
˜ 1/2 ,
˜ 1/2 =
where α˜ (0 ) =√1/2
2π f˜gg(0 ) is the square root of the initial HAR estimator, and  2πKST λ jλ j ∈B ( 0 )T   

 

−1/2
−1/2
f˜gg
λ j Igg λ j f˜gg
λj

ˆ NPW is PSD by construction.
serves as the bias correction term. ˆ NPW , let the bandwidth satisfy 1/ST + S4 /T → 0. Also deﬁne the MSE of 
ˆ NPW
For consistency of T
ˆ NPW − E vec 

ˆ NPW − W vec 

suﬃcient smoothness of the spectral density at frequency zero, the AMSE is given by

ABias2 + AV ar = (2π )2 k42 vec {ϒ(0 )} W vec {ϒ(0 )}ST−8 + tr {W(  )(I2 + K )}
where

as

for some 2 × 2 symmetric PSD weighting matrix W, as in Andrews (1991). Under∞−∞

tk2 (x )dx

S

T

T

,

 2

d
1/2
1/2
ϒ(ω ) = fgg
ω
fgg
(
(ω )
)
(ω ),
dω2
 2

d
−1/2
−1/2
f
ω
fgg
(ω ) = fgg
(ω )
(
)
(ω ),
gg
dω2

∞
tk (x ) = −∞ TK (λ )eixλ dλ, and TK (λ ) = 2K (λ ) − K ◦ K (λ ) is the fourth-order spectral window obtained by twicing (Stuetzle∞and Mittal, 1979) with K ◦ K (λ ) = −∞ K (t )K (λ − t )dt being the convolution of K (λ ). Observe that the NPW estimator ac47

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

 −2 

celerates the bias convergence from O ST
optimal bandwidth is


ST† =

 −4 

to O ST

while maintaining the variance convergence of O(ST /T ). The AMSE-

8k42 vec {ϒ(0 )} W vec {ϒ(0 )}∞
tr {W( fgg (0 ) fgg (0 ) )(I2 + K )} −∞ tk2 (x )dx

1 / 9T 1/9 .

ˆ NPW can attain the rate of convergence T 4/9 when best implemented. It is worth remarking that this result
Therefore, does not contradict Andrews’ (1991) optimality argument of the QS kernel, because the NPW estimator falls outside hisframework.
Simulation results in Hirukawa (2006) indicate that the NPW estimator often yields more accurate LRV estimates thanconventional HAR estimators. A disadvantage of this estimator is that its computational burden increases linearly with T , as

  T / 2 

1/2
f˜gg
λj

j=−T /2are required for the bias-correction term α˜ (0 ).

5.1.2. Flat-Top Kernel Estimator
Politis (2011) investigates yet another HAR estimator with a faster convergence rate. The entire estimation procedure
contains three novel features. These are: (i) using the ﬂat-top kernel; (ii) allowing for different bandwidths in different
elements of the estimate; and (iii) proposing a new procedure of estimating the optimal bandwidth.
We take a closer look at each feature. For (i), the ﬂat-top kernel is deﬁned askF T ( x ) =1g( x )

for |x| ≤ c,otherwise

where c > 0 is a parameter, the function g : R → [−1, 1] is symmetric, continuous at all but a ﬁnite number of points and∞
satisfying g(c ) = 1 and −∞ g2 (x )dx < ∞. A special case of kF T (x ) is the trapezoidal kernel

kT Z ( x ) =for |x| ≤ c
for c < |x| ≤ 1 .
otherwise

1(|x| − 1 )/(c − 1 )0

The property of this kernel becomes closer to those of the truncated and Bartlett kernels if c ≈ 1 and c ≈ 0, respectively.
Politis (2011, pp.715-716) presents other special cases of kF T (x ) including analogs to the Parzen and QS kernels.
˜ F T is given by
For (ii) and (iii), the (a, b) element of the ﬂat-top HAR estimator ˜ F T,a,b =


T −1

j=−(T −1 )



kF T

j
ST,a,bˆ j,a,b ,ˆ j,a,b is the (a, b) element of 
ˆ j . The bandwidth ST,a,b is chosen elementwise by inspecting correlograms or crosswhere correlograms, and the convergence rate of the HAR estimator is also derived in an element-by-element manner. This may
˜ F T can attain thebe thought of as a resolution for Robinson’s (1998) concern in Section 4.3. Politis (2011) establishes that 
1/
2
rate of convergence near T .
˜ F T does not necessarily
The ﬂat-top kernel has the characteristic exponent q = ∞, like the truncated kernel. Hence, generate PSD estimates under a single bandwidth. Assigning different bandwidths for different elements makes the nonPSD problem even more complicated. To generate PSD estimates while maintaining the near-parametric convergence rate,
Politis (2011) advocates replacing
the diagonal
matrix of eigenvalues  = diag {λ1 , . . . , λ } in the unitary decomposition−α˜ F T = UU∗ with + = diag λ+ , . . . , λ+ , where U is the unitary matrix, λ+
a = max {λa , T } for a = 1, . . . , , and T = T

1
+ ∗
ˆ
for some α ∈ [1, 2]. Politis (2011) ﬁnally demonstrates that the clipped ﬂat-top estimator F T = U U still maintains theconvergence rate near T 1/2 .
The ﬂat-top estimator has several theoretically attractive properties. In addition, extensive simulations in Politis (2011) indicate that the estimator often exhibits better ﬁnite-sample performances than those using conventional kernels. Nonetheless, a diﬃculty in using it in practice is to choose quite a few tuning parameters. These include: c and g(x ) for the shape of
the kernel; C0 and KT in Politis (2011, p.718) for the bandwidth selection; and the exponent α in T for the PSD correction.
5.2. Non-Autocovariance-Based Estimators
Taking a (weighted) sum of autocovariances is not the only way of handling heteroskedasticity and temporal dependenceof unknown form. Each of the nonparametric LRV estimators below is free of computing sample autocovariances.
5.2.1. Smoothed-Moments Estimator
The ﬁrst class of LRV estimators that does not depend on sample autocovariances is investigated by Smith (2005) in
the framework of eﬃcient GMM estimation. It still relies on a kernel, but smoothing is made for the error process (or the48

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

moment function) gt = gt (θ ) itself, not for its autocovariances. For a kernel k(· ) satisfying Assumption 1 and a bandwidth
ST , the smoothed moment function is deﬁned asḡt = ḡt (θ ) =

t−1

j=t−T



j
k
ST



gt− j (θ ), t = 1, . . . , T .

The idea of smoothing gt (θ ) can be also found in the context of the generalized empirical likelihood estimation using timeseries data. Because the likelihood function by construction ignores temporal dependence in gt (θ ), we run the nonparametric likelihood based on its smoothed counterpart ḡt (θ ) for valid inference; see Kitamura (1997), Kitamura and Stutzer(1997), and Smith (1997, 2011) for more details.
In the context of eﬃcient GMM estimation, Smith (2005) deﬁnes the LRV estimator as the normalized outer-product of
the smoothed process   T
  (1/T ) t=1
ḡt θ˜ ḡt θ˜
ˆS =ˆ S θ˜ =

,
T −1
2
j=1−T k ( j/ST )where θ˜ is some initial consistent GMM estimate. By construction this estimator necessarily generates PSD estimates. Anˆ −1 as the optimal weighting
interesting feature of this LRV estimator is that the eﬃcient GMM estimator (6) based on Smatrix is ﬁrst-order asymptotically equivalent to the one based on the inverse of the HAR estimator using the same initial∞∞
estimate θ˜ , the induced kernel k (x ) = −∞ k(u − x )k(u )du/ −∞ k2 (u )du and the same bandwidth ST . Smith (2011) shows thatinduced kernels from the truncated and Bartlett kernels are the Bartlett and Parzen kernels, respectively, via the relation
between their spectral window generators K  (· ) = 2π |K (· )|2 .
Smith (2005) assumes that a single kernel and bandwidth are used for all elements of gt (θ ). In fact, even if we employ
different kernels and/or bandwidths across elements like the VARHAC estimator, the resulting estimator is still PSD.
A potential disadvantage of this estimator is that Smith (2005) does not consider its implementation method. The bandwidth ST for the estimator using a kernel k(· ) may be chosen in an analogous manner to its ﬁrst-order asymptoticallyequivalent HAR estimator using the induced kernel k (· ). However, no simulation results of the LRV estimator are available
in Smith (2005), and thus how such a bandwidth choice method is transmitted to ﬁnite-sample properties of the estimator
is unknown.
5.2.2. Orthonormal Series Estimator
The second class of LRV estimators with no autocovariances is based on an orthonormal series (OS) approximation to theprocess gt , and it can be computed as the outer-product of its explained part. This class of estimators is ﬁrst proposed in
the literature by Phillips (2005) and later investigated in the context of HAR inference by Müller (2007) and Sun (2011).
∞
The OS-HAR estimator can be constructed as follows. Let {φk (r )}k=1 be an orthonormal sequence. For a collection √  TK
ˆ k = 1/ T
{φk (t/T )}k=1 with K ∈ N formed by the ﬁrst K members of the orthonormal sequence, deﬁne 
t=1 φk (t/T )gt
√
ˆk =
ˆ k
ˆ . 
ˆ k is approximately the regression coeﬃcient obtained by regressing gt on the regressor φk (t/T )/ T , and
and 
k

ˆ k can be interpreted as the part of the total sum of squares T gt g explained by this regression. The OS-HAR estimator
t=1t
is ﬁnally given by

1 ˆ
k .
K
K

ˆ OS =


(20)

k=1

As the sequence φ√k (t/T ), Phillips (2005), Müller (2007) and Sun (2011) take√the eigenvectors of the covariance kernel ofBrownian √motion 2 sin{π (k − 1/2 )t/T }, the type II discrete cosine transform 2 cos{π k(t − 1/2 )/T }, and a subset of cosine
ˆ OS is PSD by construction.
functions 2 cos(2π kt/T ), respectively. ˆˆ OS asForconsistency of OS , let the number of basis functions K = KT satisfy 1/KT + KT /T → 0. Also deﬁne the MSE of 22ˆˆE {vec OS − ) W vec(OS − ) for some  ×  symmetric PSD weighting matrix W, as in Andrews (1991). The AMSE
implied by this deﬁnition is

K

ABias2 + AV ar = vec (B ) W vec (B )

T

4

T

+ tr {W(  )(I2 + K )}KT−1 ,ˆ OS ) − }] is the dominant bias term of ˆ OS . The AMSE-optimal number of basis functions
where B = limT →∞ [(T /KT )2 {E (
is



KT† =

tr {W(  )(I2 + K )}4 vec (B ) W vec (B )

1 / 5

T 4/5 .

ˆ OS attains the best possible convergence rate of T 2/5 , which is the same as the one for
Under this optimal number, †
the HAR estimator using the Parzen or QS kernel. Implementation of KT can be conducted analogously to the methods
in Section 4.2.2.
49

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

Asymptotic equivalence of kernel and OS-HAR√ estimators is argued in Phillips (2005) and Sun (2011). In particular,ˆ OS using φk (t/T ) = 2 sin {π (k − 1/2 )t/T } is asymptotically equivalent to the HAR estimator
Phillips (2005) shows that using the Daniell kernelkDA (x ) =sin (π x ).
πxTherefore, approximations to the bias and variance of this OS-HAR estimator are the same as those of the Daniell HAR
estimator.
5.3. HAR Inference
5.3.1. Fixed-b Asymptotics
Unsatisfactory ﬁnite-sample performance of test statistics using HAR estimators has long been a concern. Specializing in
hypothesis testing of coeﬃcients in linear regression models, Kiefer et al. (20 0 0) and Kiefer and Vogelsang (2002) propose
an alternative method of constructing test statistics from the Bartlett HAR estimator by setting the bandwidth ST = T . Kiefer
and Vogelsang (2005) extend the analysis to HAR estimators using a general class of kernels and the bandwidth ST = bT forsome constant b ∈ (0, 1]. Such HAR estimators are no longer consistent, and thus have non-degenerate limiting distributions.
dˆ →To be more precise, LQ(b)L for an  ×  matrix L satisfying  = LL , where, for instance,

Q (b ) =2  1 1−bdr − 1b 0
0 B r B r
1 1
− b12 0 0 kQS r−s
B r B sbb( ) ( )



B(r + b)B(r ) + B(r )B(r + b)( ) ( ) drds

drfor Bartlett
for QS,and B(r ) = W(r ) − rW(1 ) and W(r ) for r ∈ [0, 1] are  × 1 vectors of Brownian bridges and standard Wiener processes,
respectively. Nonetheless, the test statistics using inconsistent HAR estimators are well-deﬁned and shown to obey nonstandard limiting distributions that explicitly depend on the kernel and b. This alternative approximation is referred to as
the ﬁxed-b asymptotics. The conventional approximation corresponds to the case with b → 0 ⇔ ST = o(T ), and it is called
the small-b asymptotics. The limiting distribution in this scenario is normal or chi-squared, regardless of the kernel andbandwidth.
Kiefer and Vogelsang (2005) in their simulations show that the ﬁxed-b limiting distributions of test statistics can approximate their ﬁnite-sample distributions more accurately than the standard small-b limiting distributions. Jansson (2004) andSun et al. (2008) conﬁrm this ﬁnding theoretically in a Gaussian location model (i.e., yt = β + ut with a linear Gaussian process ut ), demonstrating that the error in rejection probability (ERP) of the t-test using critical values from the nonstandard
ﬁxed-b limiting distribution is O(log T /T ) and O(1/T ), respectively. Moreover, the ﬁxed-b theory is extended beyond locationand linear regression models. Examples include Sun (2011, 2014a, 2014b) for trend regression, ﬁrst-step and eﬃcient GMM
frameworks, respectively, as well as the applications to be discussed in Section 5.4.
Kiefer and Vogelsang (2005) also ﬁnd the size-power trade-off between small and large values of b. More speciﬁcally,a small b results in higher power at the cost of larger size distortions, whereas a large b leads to greater accuracy in sizein exchange for lower power. A rationale is that size distortions (power losses) come from the bias (variance) of the LRV
estimate due to a small (large) b. The size-power trade-off leads to the proposal of yet another HAR inference below.
5.3.2. Fixed-ρ Asymptotics for Exponentiated Kernels
To improve power of test statistics while retaining accuracy in size in the framework of no truncation in the bandwidth,
Phillips et al. (20 06, 20 07) propose to exponentiate a kernel. For a kernel k(· ) (called the mother kernel) and an exponentρ ∈ N, the exponentiated kernel is given by kρ (· ) = kρ (· ). Widely used kernels including the Bartlett, Parzen and QS kernels
can be employed as the mother kernel.
The bandwidth for HAR inference using exponentiated kernels is set equal to T , as in Kiefer et al. (20 0 0) and Kiefer
and Vogelsang (2002). The long bandwidth contributes to attractive size properties. The exponent ρ controls the amount of
downweighting; as ρ increases, kρ (· ) becomes more concentrated in the vicinity of the origin than the mother kernel k(· ).
Therefore, power properties improve with ρ .
Phillips et al. (20 06, 20 07) investigate two approximations analogous to small- and ﬁxed-b asymptotics. First, lettingρ → ∞ at a suitable rate (called the large-ρ asymptotics) leads to consistency of the HAR estimators and size distortions oftest statistics, as in small-b cases. Second, a ﬁxed ρ (called the ﬁxed-ρ asymptotics) yields inconsistent HAR estimators. Like
ﬁxed-b cases, test statistics have greater accuracy in size than the conventional and large-ρ -based ones while preserving
power. An intuition is that the asymptotic theory that reﬂects that ρ is always ﬁnite in practice may be able to provide a
better approximation to the ﬁnite-sample distribution. Sun et al. (2011) establish that the ERP of the t-test for a Gaussian
location model using critical values from the nonstandard ﬁxed-ρ limiting distribution is again O (1/T ).
5.3.3. Test-Optimal Bandwidth Selection
Practitioners prefer a method that is simple and straightforward to implement, and HAR inference is not an exception.
Although the ﬁxed-b (and ρ ) asymptotics theoretically reﬁne HAR inference, practitioners may wonder how to choose b andwhat numbers should be reported as standard errors for parameter estimates. As regards the second question, under the50

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

literally ﬁxed-b setup, the HAR estimator is inconsistent and thus the resulting standard errors cannot be interpreted as
estimates for standard deviations of the distributions of parameter estimates.
Sun et al. (2008) propose a method of choosing b for two-sided tests based on HAR estimators using a general class ofkernels in a Gaussian location model. Assuming that b = bT satisﬁes bT + 1/(bT T ) → 0, Sun et al. (2008) adopt a weighted
sum of Type I and Type II errors as the loss function. For the
exponent of the kernel q, the loss function is characteristic−q†of O bT + (bT T ), and thus the test-optimal bT is bT = O T −q/(q+1 ) . The corresponding test-optimal bandwidth is givenby ST = bT T = O T 1/(q+1 ) . Observe that this is larger by an order of magnitude than the AMSE-optimal bandwidth for the
†

†





same HAR estimator ST = O T 1/(2q+1 ) .
While the test-optimal bandwidth by Sun et al. (2008) may be viewed as a purely theoretical result, Lazarus et al.
(2018) take their idea one step further. A weighted sum of the squared
 size distortions and the squared size-adjusted power−q
is chosen as their loss function, which is again of O bT + (bT T ). Then, the test-optimal formula of bT is derived undersome weighting scheme and dependent structure of the error term in the Gaussian location model. From a theoretical pointof view and extensive simulation evidence, Lazarus et al. (2018) recommend√ two practices. First, the combination of the OSHAR estimator (20) using the type II discrete cosine transform φk (t/T ) = 2 cos {π k(t − 1/2 )/T }, the test-optimal number†of basis functions KˆT = 0.4T 2/3 and ﬁxed-b critical values should be used for inference about a single parameter. Second,†the combination of the Bartlett HAR estimator, the test-optimal bandwidth Sˆ = 1.3T 1/2 and ﬁxed-b critical values should†T

be preferred for high-dimensional inference with multiple restrictions. Indeed, simulation results in Lazarus et al. (2018) indicate that the recommended procedures substantially reduce size distortions observed in HAR inference under automatic
bandwidths and normal approximations. It is also worth emphasizing that the HAR estimators based on the procedures by
Lazarus et al. (2018) are consistent, and thus we may safely report the standard errors computed from the corresponding
LRV estimates.
As a by-product, Lazarus et al. (2018, p.547) ﬁnd that the QS kernel attains the minimum value of their loss function
among all nonnegative kernels. Sun and Yang (2020) demonstrate that the test-optimality of the QS kernel continues tohold under the loss function by Sun et al. (2008).
5.3.4. HAR Inference Based on a Fixed Number of Basis Functions
The test-optimal method of choosing the bandwidth or the number of basis functions has been established in the asymptotic framework. Whether asymptotic results can be well transmitted in ﬁnite samples is a concern. Two inference procedures below are designed in the context of ﬁnite samples; these are built on some OS-HAR estimator using a ﬁxed numberof basis functions K ∈ N.
The ﬁrst procedure is Ibragimov and Müller’s (2010) split-sample t-test for a single parameter. Their inference
procedure starts from splitting T observations into K non-overlapping sub-samples of (approximately) equal length.
Let β be the parameter of interest and βˆ j be the estimator from the sub-sample j ∈ {1, . . . , K }. Then, the usual t√ statistic tβ = K βˆ − β0 /sβˆ is constructed to test for the null hypothesis H0 : β = β0 , where βˆ = K −1 Kj=1 βˆ j and s2ˆ =(K − 1 )−1 Kj=1β2βˆ j − βˆ . A critical value of t (K − 1 ) can be used for the t-test as it attains the upper bound for rejection

probability under the null; to put it another way, this test becomes conservative under different variances across K subsamples and has exact size under equal variances (as in Student’s original approach). Superior ﬁnite-sample properties ofthis procedure in comparison with ﬁxed-b approaches are conﬁrmed in Monte Carlo simulations by Ibragimov and Müller
(2010) and Müller (2014). It is worth remarking that Ibragimov and Müller’s (2010) test can be interpreted as a special caseof HAR inference. At a ﬁrst glance, the relation of this test with HAR inference is not obvious because of no explicit use
ˆ SS = (T /K )s2 can be expressed as anof a LRV estimate. However, Lazarus et al. (2021, Proposition S1) demonstrate that ˆβequal-weighted OS-HAR estimator (20) using K basis functions of some orthonormal sequence.
Second, Dou (2020) studies optimal HAR inference in ﬁnite samples for a Gaussian location model. The test statistic isagain the t-statistic, where the standard error is computed by the OS-HAR estimator (20) using K type II discrete cosine
transforms. The critical value for a given level of signiﬁcance can be obtained as an adjustment factor times the criticalvalue from t (K ), where the adjustment factor is greater than 1 and thus the usual Student-t critical value is necessarily
inﬂated for implementation. An advantage of these two tests over ﬁxed-b tests is no need to run simulations for criticalvalues, because Ibragimov and Müller’s (2010) test is Student-t based and the adjustment factor for Dou’s (2020) test canbe easily tabulated.
5.4. LRV-Based Inference for Nonstationary Data
Apart from HAR inference discussed above, inference using LRV estimates for (possibly) nonstationary data has evolved
uniquely. This section specializes in unit-root and stationarity tests and structural break tests using LRV estimates.
Ng and Perron (2001) and Vogelsang and Wagner (2013) are two improvements of the Phillips-Perron unit-root test
(Phillips and Perron, 1988) in different directions. For an AR(1) model yt = α yt−1 + ut with a serially correlated error ut , thePhillips-Perron test statistic requires a consistent estimate of the LRV of ut . Perron and Ng (1996) ﬁnd that the test statistic
using an AR spectral density estimator (i.e., a univariate version of the VARHAC estimator (12)) has less size distortions than51

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

the one using a kernel HAR estimator. To improve size and power properties of the Phillips-Perron test for the regression

yt = zt β + α yt−1 + ut

(21)

with deterministic trends zt = () for m ∈ N and a serially correlated error ut , Ng and Perron (2001) propose tocombine the detrending procedure by Elliott et al. (1996) and a modiﬁed information criterion for the lag order selection of(12) that replaces the existing criteria such as AIC (10) and BIC (11).
In contrast, Vogelsang and Wagner (2013) extend the ﬁxed-b asymptotics to the Phillips-Perron test for (21) using kernel
HAR estimators. It is demonstrated that one- and two-step detrending procedures for (21) yield different non-pivotal ﬁxedb limits. Based on this, Vogelsang and Wagner (2013) propose modiﬁed Phillips-Perron tests that allow for asymptotically
pivotal ﬁxed-b inference.
Next, to reduce ﬁnite-sample bias in the KPSS test statistic for stationarity (Kwiatkowski et al., 1992) with prewhitened
HAR estimators employed, Sul et al. (2005) propose to combine two practices. First, Sul et al. (2005) introduce recursive
demeaning procedures to correct ﬁnite-sample bias in prewhitening coeﬃcient estimates. Second, Andrews and Monahan(1992) recommend the so-called ‘0.97’ rule that replaces roots of the ﬁtted characteristic
equation by 0.97 whenever they√exceed this number. In place of 0.97, Sul et al. (2005) advocate the ‘1 − 1/ T ’ boundary condition rule. These practices
jointly improve size and power properties of the KPSS test while maintaining its consistency.
Finally, the test statistic for the null of stability in the mean-shift model
1, t, . . . , t m

yt = μ + δ 1(t ≥ t ∗ ) + ut

(22)

with an unknown break date t ∗ and a serially correlated error ut

again requires a consistent estimate of the LRV of ut . It
has long been known that empirical rejection rates of the test statistic using a kernel HAR estimate computed under thenull do not increase with the magnitude of the break |δ|. Employing a HAR estimate computed under the alternative of abreak may be a remedy for the problem of non-monotonic power, whereas this comes at the cost of severe  size distortions.
Taking these aspects into account, Kejriwal (2009) proposes to estimate the instantaneous variance E ut2 and bandwidth

ST using residuals under the alternative while estimating autocovariances E ut ut− j using those under the null. Kejriwal(2009) reports that the test statistic using the hybrid HAR estimate can bypass the non-monotone power problem while
maintaining adequate size properties.
Furthermore, Wenger and Leschinski (2021) extend the ﬁxed-b asymptotics to the CUSUM test for a break in the model(22) with a long-memory error ut . Comparing two alternative kernel HAR estimates using residuals under the null andalternative, Wenger and Leschinski (2021) ﬁnd that those under the alternative result in good size and power of their test.
5.5. Estimation-Optimal Bandwidth Selection
The problem of bandwidth selection in kernel HAR estimation is not limited to statistical inference. Wilhelm (2015) proposes yet another objective-oriented bandwidth choice method. Assume that  > p, i.e., the parameter of interest θ0 is overidentiﬁed. Based on a higher-order approximation to the MSE of the two-step GMM estimator, the method selects the op†
timal bandwidth ST for the parameter estimator of interest, not the one for the kernel HAR estimator. Speciﬁcally, Wilhelm(2015, Proposition 1) establishes that the eﬃcient GMM estimator (6) admits the stochastic expansion√ T θˆ − θ0 = ZT + BT ST−q + VT

!

"

ST
+ o p ST−q +
T

!

#

ST
T

,

where q is the characteristic exponent of the kernel,



ZT = − D −1 D

BT = D −1 D

VT = D 

−1

D

−1

−1
−1

√
D −1 T G(θ0 ),ˆ −D −1 STq E 

!

D 

−1



T ˆ
ˆ
−E 
ST

√

−1
−1 − −1 D D −1 D D −1 T G(θ0 ), and√−1
−1 − −1 D D −1 D D −1 T G(θ0 ).

Three O p (1 ) quantities ZT , BT and VT correspond to an asymptotically normal term, and terms involving the bias and vari$ˆ , respectively. It follows from BT S−q + VT ST /T = o p (1 ) that estimation errors of ˆ do notance of the HAR estimator Taffect the ﬁrst-order asymptotic result of the eﬃcient GMM estimator θˆ, as is well known.
Following Andrews (1991), Wilhelm (2015) deﬁnes the (scaled) MSE of θˆ as E {T (θˆ − θ0 ) W(θˆ − θ0 )} for some p × p
symmetric PSD weighting matrix W. Then, it can be found that the optimal bandwidth that balances the squared ‘bias’ and†−2q‘variance’ terms of order ST
and ST /T in the MSE is ST = C † (q )T 1/(2q+1 ) for some constant C † (q ). Therefore, the optimal
bandwidth in eﬃcient GMM estimation shares the same expansion rate with the optimal bandwidth in HAR estimation.
Nevertheless, the coeﬃcient C † (q ) substantially differs, for instance, from what Andrews (1991) derives; see Section 2.1 ofWilhelm (2015) for a comparison of two coeﬃcients.
Because C † (q ) is highly complicated, Wilhelm (2015) simpliﬁes its estimation procedure by making a few additional
assumptions on the dependent structure of gt . After assuming that gt obeys a linear Gaussian process, as in Andrews (1991),
52

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

Wilhelm (2015) puts W = I p and ﬁts an AR(1) model to each element of gt to obtain proxies of several unknown quantities
in C † (q ). Wilhelm (2015) also reports nice ﬁnite-sample properties of the eﬃcient GMM estimator θˆ using the estimationoptimal bandwidth. However, the simpliﬁed implementation method could be a disadvantage if gt deviates from Gaussianityor the ﬁtted AR(1) model turns out to be misspeciﬁed, as is the case with Andrews’ (1991) AMSE-optimal bandwidth.
5.6. LRV Estimation beyond Covariance Stationarity
Consistency of LRV estimators is typically built on covariance stationarity with ﬁnite fourth-order moments of gt in theearly literature. Attempts to relax these conditions have been made since then. Examples include consistency proofs byHansen (1992b) for non-covariance stationary, heterogenous processes and by De Jong and Davidson (20 0 0) and Davidson
and De Jong (2002) for near-epoch dependent functions on mixing processes. LRV estimators in the recent literature assume
yet other classes of dependent structure. Two examples are presented below.
5.6.1. Locally Stationary Processes
Kawka (2020) studies spectral density estimation when gt obeys a locally stationary process. The process can describe
slowly-changing dependent structure over time and is characterized by time-varying second-order moments and behavior
like a stationary process in a small neighborhood of each time point t ∈ {1, . . . , T }; see Dahlhaus (2012) for a comprehensive
review of this class of processes. To examine statistical properties of the HAR estimator smoothed by a general class ofˆ J given in (13), deﬁne u = t/T so that the time domain t ∈ {1, . . . , T } is rescaled to the unit interval u ∈ [0, 1]. Also
kernels 
let fgg (u, ω ) be the  ×  time-varying spectral density matrix of gt evaluated at rescaled time u and frequency ω ∈ (−π , π ).
1
p
ˆJ →
Kawka (2020) establishes that under some regularity conditions, 2π
fgg (u, 0 )du holds, i.e., the probability limit of0ˆ J is 2π times the time-averaged spectral density matrix evaluated at frequency zero.


5.6.2. Long-Memory Processes
Robinson (2005) proposes a smoothed nonparametric LRV estimator that is consistent against long-memory or antipersistence. Suppose that the spectral density matrix of gt satisﬁes fgg (ω ) ∼ h(ω )Sh(ω ) as ω → 0+, where S is an  ×  ﬁnite
PD matrix, h(ω ) = diag eid1 π /2 ω−d1 , . . . , eid π /2 ω−d , d1 , . . . , d ∈ (−1/2, 1/2 ) are memory parameters, and ∼ signiﬁes thatthe ratio of real parts, and of imaginary parts, of corresponding elements of matrices on the left- and right-hand sides converges to unity. The ath element of gt is a short-memory, long-memory and antipersistence processes if da = 0, da ∈ (0, 1/2 )and da ∈ (−1/2, 0 ), respectively. Then,



 −d

V ar diag T

1,...,T



T
 1 
−d

gt
√
T t=1



→ V = 2π S ◦ Q,

where the (a, b) element of the  ×  matrix Q = Q(d1 , . . . , d ) issin (π da ) + sin (π db ),(da + db + 2 ) sin {π (da + db )}(· ) is the gamma function, and ◦ signiﬁes the Hadamard product.
Robinson’s (2005) memory and autocorrelation con
ˆ = 2π Sˆ ◦ Q
ˆ = Q dˆ1 , . . . , dˆ , dˆ1 , . . . , dˆ are consistent estimates of memory
ˆ , where Q
sistent (MAC) estimator is deﬁned as V
parameters,Sˆ =12mT − 1

 

m
T −1

   −1−1hˆ λ jIgg λ j hˆ λ j

j=−(mT −1 )

 

 

for the bandwidth mT is a natural estimator of S, Igg λ j is the periodogram evaluated at frequency λ j = 2π j/T , and hˆ λ j =
−dˆ

−dˆ

ˆ is PSD by construction. Moreover, it can be demonstrated that if
diag eid1 π /2 λ j 1 , . . . , eid π /2 λ j  . The MAC estimator V
ˆ

ˆ



p

2

ˆ → V.
mT → ∞ and mT = o T / log T , then V
Abadir et al. (2009) compare the asymptotic expansion of the HAR estimator using the Bartlett kernel with that of the†
MAC estimator for a scalar process having the memory parameter d ∈ (−1/2, 1/2 ). The AMSE-optimal bandwidth ST for the
HAR estimator is shown to depend on the memory d so thatST† =

  1/(3+4d ) O T 1/2−d O T

for d ∈ (−1/2, 1/4 ).
for d ∈ [1/4, 1/2 )† 4/5 Observe that the optimal bandwidth reduces to the familiar one ST = O T 1/3 for the short-memory case (d = 0). In contrast,††

the AMSE-optimal bandwidth mT for the MAC estimator is mT = O T
regardless of d. Abadir et al. (2009) conclude that
the MAC estimator is more robust to the bandwidth choice than the HAR estimator as the expansion rate of the optimal
bandwidth for the latter is inﬂuenced by d.
53

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

5.7. Implementation in Statistical Packages
This section concludes by noting the commands for implementing HAR estimation and inference in standard statisticalpackages. In R, the package cointReg contains two functions getLongRunVar and getBandwidth for conventional
kernel HAR estimation. The former computes kernel HAR estimates, and the latter yields automatic bandwidths by Andrews(1991) and Newey and West (1994). In Stata, the command lrcov by Wang and Wu (2012) corresponds to these twofunctions.
The commands for ﬁxed-b inference are also available in Stata. Ye and Sun (2018) develop several commands based on
test-optimal bandwidths and ﬁxed-smoothing critical values. More speciﬁcally, F and t tests for linear regression models (foreﬃcient linear GMM estimation) can be implemented with their commands har (gmmhar) and hart (gmmhart).
6. Monte Carlo Simulations
While test-optimal bandwidth selection methods have been explored for a few decades, estimation-optimal ones are still
under development. Focusing on the latter, we conduct two independent simulation studies. One is based on eﬃcient GMM
estimation and the other on cointegrating regressions. In these applications, choices of LRV estimators do not matter in theﬁrst-order asymptotic sense but may affect ﬁnite-sample performances of estimates of model parameters. Therefore, each of
two Monte Carlo studies is concerned with how a variety of LRV estimators and/or bandwidth formulae can inﬂuence thequality of parameter estimates in ﬁnite samples.
6.1. Design #1: GMM
The focus of the ﬁrst Monte Carlo study is on eﬃcient GMM estimation. As already seen in the previous section, Wilhelm(2015) proposes an estimation-optimal bandwidth choice method for HAR estimators using a general class of kernels in the
GMM framework. However, its implementation is based on several stringent assumptions. In particular, practitioners may
wonder whether gt is not a linear Gaussian process or whether the ﬁtted AR(1) model turns out to be misspeciﬁed.
In light of these concerns, the Monte Carlo design incorporates the cases in which gt deviates from Gaussianity or AR(1),whereas it basically follows Wilhelm (2015). Consider a slope-only linear regression model yt = θ xt + vt with the true valueθ = 1. The regression error vt falls into one of the two cases



vt =

u1t
1
x u
4 | t | 1t

[HOM]
,
[HET]

where “HOM” and “HET” stand for “homoskedastic” and “heteroskedastic” errors, respectively. The regressor xt has the reduced form xt = π zt + u2t with a ﬁve-dimensional vector of unity π . The bivariate process ut = (u1t , u2t ) obeys either aVAR(1) modelut = 0.8I2 ut−1 +t [VAR]or a VMA(1) modelut =t + 0.8I2 t−1 ,where the innovation  

iidt ∼ N2[VMA]

t = ( 1t , 2t )

0
1
,
0
0.8



0.8
1

is generated by.

Due to the correlation between 1t and 2t , xt is correlated with vt and thus endogenous. Finally, the vector of instruments zt = (z1t , . . . , z5t ) is also generated by a VAR(1) model zt = z zt−1 + wt , where z = diag {0, 0.8, . . . , 0.8}, wt =(w1t , . . . , w5t ) , and w1t , . . . , w5t are iid copies of N (0, 1 ), χ 2 (2 ) − 2 or t (5 ) and independent of s for all t, s. For eachcombination of vt , ut and wt , 50 0 0 data sets of sample sizes T = 128 are simulated.
Observe that gt = gt (θ ) = zt (yt − θ xt ) satisﬁes the orthogonality condition E (gt ) = 05×1 . Taking the two-stage leastsquares (2SLS) estimator θˆ2SLS as the initial estimator, we have the eﬃcient GMM estimator of θ in a closed form"θˆGMM =

T

t=1

#

ˆ −1
xt zt 

"

T

t=1

#−1 "

zt xt

T


#

"

ˆ −1
xt zt 

t=1



T


#

zt yt

t=1



ˆ computed from gˆ t = gt θˆ2SLS = zt yt − θˆ2SLS xt .
for some LRV estimator ˆ , in addition to QS-A, BT-NW and PZ-H used in Section 4.4, we consider three versions of Smith’s (2005) smoothedAs 
moments estimator using the truncated kernel [S-BT-A, S-BT-NW, S-BT-H] and the QS, Bartlett and Parzen HAR estimators
54

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

using Wilhelm’s (2015) estimation-optimal bandwidths [QS-W, BT-W, PZ-W]. Here are some details on Smith’s (2005) estimator. As in Example 2.1 of Smith (2011), the smoothed moment function using the truncated kernel can be constructedbyḡt = ḡt (θ ) =

(t−1
)∧mT


2
2mT + 1

gt− j (θ ), t = 1, . . . , T .

j=(t−T )∨(−mT )†The optimal bandwidth mT is implied by the induced kernel of the truncated kernel, which is the Bartlett kernel. Because no†further information on implementing the smoothed-moments estimator is available, we estimate mT indirectly with the aidof automatic bandwidths for the Bartlett HAR estimator by Andrews (1991), Newey and West (1994) and Hirukawa (2010),††

denoted by BT-A, BT-NW and BT-H, respectively. More speciﬁcally, given the relation ST = 2mT + 1 /2 and an automatic
†

†

†

ˆT =
bandwidth SˆT , mT can be estimated as m



†



2SˆT − 1 /2 .

Performance of each estimator is measured by the root MSE (RMSE), bias (BIAS) and standard deviation (SD). Monte
Carlo averages of bandwidths (BW) are also presented for convenience. The results of θˆ2SLS are provided as a benchmark.
Table 3 reports simulation results. Except the three cases of VMA-HOM, RMSEs of GMM estimates are smaller than thoseof the corresponding 2SLS estimates, as the theory predicts. Although there is no dominant estimator, PZ-W, QS-W and BT-W
tend to attain the smallest RMSE for VAR-HOM, VAR-HET and VMA-HET, respectively. Therefore, Wilhelm’s (2015) bandwidthappears to have some degree of robustness against the deviations from Gaussianity and AR(1).
A few interesting phenomena can be also found. First, Wilhelm’s (2015) bandwidths are in general smaller than thecorresponding automatic bandwidths. This ﬁnding agrees with simulation results in Wilhelm (2015), and it is a sharp contrast to test-optimal bandwidths, which tend to be longer than the automatic bandwidths. Second, ﬁnite-sample properties
of Smith’s (2005) smoothed-moments estimators do not surpass those of the ﬁrst-order equivalent Bartlett HAR estimator. Third but not least importantly, a careful inspection reveals that regardless of the data generation processes, average
bandwidths of S-BT-H and PZ-H are close to those of BT-W and PZ-W, respectively. As a consequence, RMSEs of the former
are comparable to those of the latter. Invoke that both S-BT-H and PZ-H are computed by the SEPI algorithm by Hirukawa(2010). This may be a hint for the direction of improving Wilhelm’s (2015) estimation-optimal approach.
6.2. Design #2: Cointegrating Regressions
The second Monte Carlo study deals with a linear cointegrating regression model with an endogenous I (1 ) regressor.
Because of the endogeneity and the serially correlated regression error, the OLS estimator of the cointegrating vector is
T -consistent but ineﬃcient due to the so-called second-order bias. To gain full eﬃciency, several authors have proposedbias-correction estimation methods of the cointegrating vector. This study concentrates on two bias-correction procedures
using LRV estimates, namely, fully modiﬁed least squares (FMLS) by Phillips and Hansen (1990) and canonical cointegrating regression (CCR) estimation by Park (1992). Finite-sample properties of FMLS and CCR are compared with those of
integrated modiﬁed OLS (IM-OLS) estimation, a LRV estimation-free bias-corrected procedure proposed by Vogelsang and
Wagner (2014).
T
The Monte Carlo design is similar to Hirukawa (2011). The bivariate process {(yt , xt )}t=1∈ R2 admits the triangular representationyt = xt θ + u1t ,xt = u2t ,(23)

where xt = (1, xt ) , θ = (θ0 , θ1 ) with true values θ0 = θ1 = 1, and θ1 is the parameter of interest. The error term ut =(u1t , u2t ) obeys either a VAR(1) model



ut =



ρ

0
u
+
0 t−1

0

or a VMA(1) model



ut =

t +

0.3

ρ

t [VAR]



−0.4
0.6

t−1 ,

[VMA]

where ρ ∈ {0.4, 0.8}. The innovation  

iidt ∼ N2

0
,
0

1

σ21

σ21
1



t = ( 1t , 2t )is generated by,where σ21 ∈ {0.4, 0.8} for VAR and σ21 ∈ {−0.4, −0.8} for VMA. For each combination of (ρ , σ21 ), 50 0 0 data sets of sample
sizes T = 128 are simulated.
55

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

Table 3
Simulation Results on GMM
HOM

HET

Model

Distribution of w

Estimator

RMSE

BIAS

SD

BW

RMSE

BIAS

SD

BW

VAR

N ( 0, 1 )

2SLS
GMM:

0.08710.08330.08390.08370.08420.08360.08340.08320.08480.08300.04390.0426
0.04300.04280.04280.04270.04280.04270.04340.0426
0.0671
0.0662
0.06620.06590.06630.06630.06620.06580.0666
0.06570.0431
0.04490.04440.04430.04490.04470.04450.04440.0441
0.04440.02180.02290.02260.02260.02290.02270.02270.02270.02250.02270.03410.03560.03530.03530.03570.03550.03530.03530.03510.03530.01470.01450.01500.01480.01430.01450.01460.01480.01510.01480.00370.00400.00390.00400.00410.00390.00390.00390.00380.00400.00710.00700.00750.00720.00700.00720.00720.00730.00730.00730.00420.00420.00420.00420.00400.00440.00440.00440.00430.00430.00110.00120.00120.00120.00130.00120.00120.00120.00120.00120.00170.00180.00180.00180.00170.00180.00180.00180.00180.00180.08590.08210.08250.08240.08290.08240.08210.08190.08340.08170.04370.04240.04280.04260.04260.04260.04260.04250.04330.04240.06670.06590.06580.06550.06600.06590.06580.06540.06620.06530.04290.04470.04420.0441
0.04470.04450.04420.04420.04390.04420.02180.02290.02260.02260.02280.02270.02270.02270.02250.0226
0.03410.03560.03520.03530.03560.03550.03530.03530.03510.0353

9.69
6.27
11.43
10.40
6.27
4.89
5.53
4.12
11.14
9.78
6.24
11.36
10.45
6.24
4.86
5.60
4.14
11.27
9.74
6.27
11.38
10.43
6.27
4.87
5.55
4.12
11.18
5.10
5.10
6.43
5.80
5.10
2.37
2.98
2.30
5.99
5.17
5.16
6.42
5.86
5.16
2.36
3.03
2.34
6.10
5.15
5.24
6.41
5.84
5.24
2.36
3.01
2.33
6.07

0.17280.15800.15720.15740.16090.15870.1576
0.15690.15800.15710.13910.12860.12820.12800.13030.12900.12840.12820.12900.12810.15550.14350.14290.14280.14570.14380.14320.14270.14370.1427
0.09650.09130.09070.09040.09180.09180.09020.09000.08940.09030.07450.07070.07010.07000.07130.07080.0701
0.06990.06950.07000.08580.08150.08080.08080.08200.08180.08080.08080.08030.08090.12570.11080.10930.10980.11390.11150.10990.10900.10890.10950.06690.0596
0.05870.0592
0.06140.05990.05930.05850.05840.05890.09850.08700.08580.08610.08940.08740.08640.08560.08540.08600.07490.06790.06780.06740.06890.06900.06700.0667
0.06620.06710.03890.03650.03610.03580.03700.03670.03580.03570.03530.03590.0581
0.0531
0.05290.05280.05390.05370.05230.05230.05190.05260.11850.1126
0.11300.11280.1136
0.11300.11290.11290.11440.11260.12190.11400.11390.11350.11490.1142
0.11390.11400.1151
0.11380.12030.11420.1142
0.11390.11510.11420.11420.11420.11550.11390.06080.06090.06030.06020.06060.06050.06030.06040.06010.06040.06350.06060.06000.06010.06090.06060.06030.06010.05990.0601
0.0632
0.06180.0612
0.0612
0.06180.0617
0.0616
0.0616
0.0612
0.0615

7.18
5.81
9.74
7.92
5.81
4.06
4.12
3.08
8.29
7.16
5.87
9.72
7.85
5.87
4.07
4.15
3.09
8.36
7.24
5.84
9.66
7.94
5.84
4.02
4.16
3.09
8.37
4.33
6.00
5.97
4.87
6.00
1.95
2.58
2.00
5.19
4.43
6.15
6.15
4.94
6.15
2.00
2.67
2.07
5.38
4.41
6.51
6.03
4.92
6.51
1.99
2.65
2.05
5.33

χ 2 (2 ) − 2

t (5 )VMA

N ( 0, 1 )χ 2 (2 ) − 2

t (5 )

2SLS
GMM:

2SLS
GMM:

2SLS
GMM:

2SLS
GMM:

2SLS
GMM:

QS-A
BT-NW
PZ-H
S-BT-A
S-BT-NW
S-BT-H
QS-W
BT-W
PZ-W
QS-A
BT-NW
PZ-H
S-BT-A
S-BT-NW
S-BT-H
QS-W
BT-W
PZ-W
QS-A
BT-NW
PZ-H
S-BT-A
S-BT-NW
S-BT-H
QS-W
BT-W
PZ-W
QS-A
BT-NW
PZ-H
S-BT-A
S-BT-NW
S-BT-H
QS-W
BT-W
PZ-W
QS-A
BT-NW
PZ-H
S-BT-A
S-BT-NW
S-BT-H
QS-W
BT-W
PZ-W
QS-A
BT-NW
PZ-H
S-BT-A
S-BT-NW
S-BT-H
QS-W
BT-W
PZ-W

Note: “RMSE”, “Bias”, “SD”, and “BW” are the RMSE, bias, standard deviation, and Monte Carlo average of bandwidths, respectively. Averages of automaticˆ †T , are reported as those for smoothed-moments estimators.
bandwidths for the Bartlett estimator SˆT† , not those of m56

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

In this setup θˆOLS , the OLS estimator of the cointegrating vector θ in (23), is T -consistent but ineﬃcient due to thesecond-order bias. FMLS and CCR are developed as eﬃcient estimation methods that can eliminate the bias, and these are
ﬁrst-order asymptotically equivalent. Let  and  be the two- and one-sided LRVs of ut , respectively, where the deﬁnitionof  is given in (7) and these LRVs can be partitioned into





12
11
and  =
22
21

11
=
21

  
12
1
=
.
22
2



ˆ and ˆ be their consistent estimators with ut replaced by u
ˆ t = uˆ1t , xt
Also let (23). The FMLS estimator of θ is given by"θˆF MLS =

T


#−1 "

zt zt

t=1

T


ˆ0 =
where as

zt yt+ − T ˆJ+ ,t=1

ˆ2
ˆ −1 
0Tˆtu
ˆ t /T
t=1 u"

θˆCCR =

T

t=1

 



ˆ 21 − ˆ −1 
ˆ
ˆ 22 
0, 
22 21

ˆ t and yt∗ = yt −u

. CCR employs the transformed data
ˆ 2 θˆOLS + 0, 
ˆ 12 
ˆ −1
ˆ −1 

0
22


ˆt,
u



is a consistent estimator of 0 = E (ut ut ). Denoting zt∗ = 1, xt∗#−1 "zt∗ zt∗, where uˆ1t is the OLS residual from#ˆ −1 xt and ˆJ+ =ˆ 12 
where yt+ = yt − 22xt∗ = xt −



T




yields the CCR estimator of θ#zt∗ yt∗.

t=1

A problem is that so far no estimation-optimal bandwidth choice method for a HAR estimator has been proposed in
the context of cointegrating regressions, unlike the GMM case, to the best of our knowledge. As a compromise, we again
ˆtemploy HAR estimators of  using currently available automatic bandwidths, namely, QS-A, BT-NW and PZ-H, where u
may be prewhitened by a VAR(1) ﬁlter before HAR estimation. The one-sided LRV estimator based on VAR(1)-prewhitening
can be found% in &Hansen (1992a, p.47). In addition, the VARHAC estimator of  is employed, where the maximum lag orderˆ t is chosen via BIC (11). The one-sided LRV estimator based on
M = MT = T 1/3 and the lag order in each element of u
VARHAC is implied by Proposition 2 of Park and Ogaki (1991).
Instead of LRV estimation, the second-order bias correction in IM-OLS is made by a partial sum transformation of theregression (23). The IM-OLS estimator of θ is deﬁned as the OLS estimator of the transformed regression 

Sty = Stx

θ + xt γ + Stu1 ,uywhere St = tj=1 yt , Stx and St 1 are deﬁned analogously, and xt is added as an extra regressor to establish asymptotic mixed
normality when (u1t , u2t ) are correlated.
As before, RMSE, BIAS and SD are chosen as performance measures. Monte Carlo averages of bandwidths and lag orders(BW/LO) are also reported. The results of θˆOLS are provided as a benchmark.
Finite-sample performances of OLS, IM-OLS, FMLS, and CCR estimators of θ1 are presented in Table 4. FMLS and CCR
decrease RMSEs from OLS, except those using VARHAC. Performances of FMLS and CCR do not have much difference, andprewhitening often contributes bias reduction. Again results are mixed, and no dominant LRV estimator can be found.
VARHAC does not look an appropriate tool for FMLS or CCR. It yields the highest RMSE for all eight cases. In addition,
RMSEs of VARHAC often exceed those of OLS for four cases of VMA. Disappointing performance of VARHAC is attributed
to both large BIAS (in magnitude) and SD. This reﬂects diﬃculty in estimating AR coeﬃcients precisely, and coeﬃcientestimates with poor quality in turn induce additional bias and/or variability in estimates of cointegrating vectors.
IM-OLS is more effective in bias reduction than FMLS and CCR using non-prewhitened (in all cases) and prewhitened
HAR estimates (frequently). However, its attractive bias properties come at the cost of higher SD and thus RMSE. Theseﬁndings coincide with what are reported in Vogelsang and Wagner (2014).
7. Concluding Remarks
The literature on HAR estimation and inference is still growing. This article has discussed and compared a variety of LRV
estimators that have been proposed in the past few decades, and it covers both parametric and nonparametric estimators.
Particular attention is paid to whether a LRV estimator necessarily generates PSD estimates. Kernel methods have been dominant among all LRV estimation procedures. A section is devoted to recent developments in HAR estimation and inference
including ﬁxed-b asymptotics and test-optimal bandwidth selection. This article also covers estimation-optimal bandwidth
choice methods for HAR estimation, which are still under development. In this connection, we run Monte Carlo simulations
implied by two estimation problems using LRV estimators, namely, GMM and cointegrating regressions.
57

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

Table 4
Simulation Results on Cointegrating Regressions

ρ = 0.4

ρ = 0.8

Model

σ21

Estimator

RMSE

BIAS

SD

BW/LO

RMSE

BIAS

SD

BW/LO

VAR

0.4

OLS
IM-OLS
FMLS:

0.05030.05700.04150.04080.04060.04080.04170.04080.04600.04140.04070.04050.04070.04150.04070.04230.07270.04770.03420.03540.03410.02900.03660.0291
0.09820.03530.03600.03520.03210.03750.03210.0396
0.03740.03570.02600.02690.02610.02750.02680.02750.05190.02610.02690.02610.02620.02620.02630.04380.05100.03640.02530.02520.02430.02280.02300.02230.09930.02560.02540.02470.02210.02300.02220.0496

0.02640.00280.00670.00820.00780.00350.00550.00350.00570.00650.00800.00760.00360.00550.00360.0061
0.05320.00520.01330.01630.01540.00560.01120.0056−0.00850.01510.01730.01670.01070.01400.01070.0084−0.0225−0.0021−0.0044−0.0072−0.0058
0.00980.00320.0094−0.0115−0.0048−0.0073−0.0060
0.00780.00220.0076−0.0137−0.0357−0.0036−0.0073−0.0089−0.0076
0.00620.00080.0036−0.0068−0.0080−0.0093−0.0082
0.0020−0.00190.0002−0.0208

0.04280.05700.04100.04000.03990.04070.04130.04070.04560.04090.03990.03980.04050.04120.04050.04190.04950.04740.03150.03140.03040.02850.03490.0286
0.09780.03190.0316
0.03100.03030.03480.03030.03870.02990.03560.02560.02590.02540.02570.02660.02590.05060.02560.02590.02540.02500.02620.0251
0.0416
0.03650.03620.02420.02350.02300.02190.02290.02200.09900.02430.02370.02330.02200.02290.02220.0450

4.48
4.49
4.96
0.96
5.18
0.14
1.01,1.01
4.48
4.49
4.96
0.96
5.18
0.14
1.01,1.01
4.67
4.41
5.05
0.94
5.19
0.10
1.02,1.01
4.67
4.41
5.05
0.94
5.19
0.10
1.02,1.01
4.27
4.40
5.83
1.86
4.75
2.18
1.55,1.49
4.27
4.40
5.83
1.86
4.75
2.18
1.55,1.49
4.91
4.56
5.69
1.96
4.49
4.62
1.61,1.12
4.91
4.56
5.69
1.96
4.49
4.62
1.61,1.12

0.12980.16610.12570.12160.12240.12580.12700.12600.13080.12510.12130.12200.12340.12490.12360.12810.18230.14570.13210.12730.1245
0.09620.10690.09630.21800.13600.13120.12980.11820.12280.11760.17480.04280.03810.02890.03100.02960.02970.02880.02910.04950.02900.03110.02970.02880.02850.02840.04600.07760.05020.03990.04170.04100.03490.03510.03480.11170.04050.04220.04170.03360.03530.03530.08850.06950.0216
0.03640.0406
0.03820.01840.02130.01870.06580.03770.04090.03890.02320.02550.02340.06690.13930.04180.07440.07980.07460.02530.03490.02470.10700.08210.08440.0816
0.06380.06580.06310.1269−0.0269−0.0026−0.0056−0.0103−0.00820.00930.00350.0071−0.0168−0.0060−0.0105−0.0084
0.00790.00260.0058−0.0182−0.0565−0.0049−0.0129−0.0181−0.0169
0.0073−0.0014−0.0025−0.0486−0.0143−0.0189−0.0179
0.0020−0.0049−0.0057−0.0529

0.10960.16460.12030.11460.11630.12450.12520.1246
0.11310.11930.11420.11560.12120.12220.12130.1092
0.11760.1396
0.10920.09930.0997
0.09280.10100.09310.18990.10840.10050.10100.09960.1036
0.09930.12030.03330.03800.02840.02920.02840.02820.02860.02820.04650.02840.02930.02850.02770.02830.02780.04230.05320.05000.03780.03750.03740.03410.03510.03470.1006
0.03790.03780.03770.0336
0.03490.03480.0710

13.77
7.57
13.69
1.20
4.30
0.43
1.01,1.01
13.77
7.57
13.69
1.20
4.30
0.43
1.01,1.01
15.04
7.68
14.08
1.10
4.40
0.26
1.02,1.01
15.04
7.68
14.08
1.10
4.40
0.26
1.02,1.01
3.81
4.38
5.84
2.00
4.33
4.37
1.39,2.47
3.81
4.38
5.84
2.00
4.33
4.37
1.39,2.47
5.12
4.30
4.36
1.66
4.69
7.12
1.35,1.89
5.12
4.30
4.36
1.66
4.69
7.12
1.35,1.89

CCR:

0.8

OLS
IM-OLS
FMLS:

CCR:

VMA

−0.4

OLS
IM-OLS
FMLS:

CCR:

−0.8

OLS
IM-OLS
FMLS:

CCR:

QS-A
BT-NW
PZ-H
QS-A-PW
BT-NW-PW
PZ-H-PW
VARHAC
QS-A
BT-NW
PZ-H
QS-A-PW
BT-NW-PW
PZ-H-PW
VARHAC

QS-A
BT-NW
PZ-H
QS-A-PW
BT-NW-PW
PZ-H-PW
VARHAC
QS-A
BT-NW
PZ-H
QS-A-PW
BT-NW-PW
PZ-H-PW
VARHAC

QS-A
BT-NW
PZ-H
QS-A-PW
BT-NW-PW
PZ-H-PW
VARHAC
QS-A
BT-NW
PZ-H
QS-A-PW
BT-NW-PW
PZ-H-PW
VARHAC

QS-A
BT-NW
PZ-H
QS-A-PW
BT-NW-PW
PZ-H-PW
VARHAC
QS-A
BT-NW
PZ-H
QS-A-PW
BT-NW-PW
PZ-H-PW
VARHAC

Note: HAR estimators with “PW” indicate prewhitened ones. “RMSE”, “Bias”, “SD”, and “BW/LO” are the RMSE, bias, standard deviation, and Monte Carloˆ t are presented.
average of bandwidths or lag orders, respectively. For VARHAC, average lag orders of the ﬁrst and second elements of u58

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

This article concludes with possible directions of future research suggested by the Monte Carlo simulations. The ﬁrst
simulation study conﬁrms some degree of robustness in Wilhelm’s (2015) estimation-optimal bandwidth choice method for
GMM. While it is implemented like Andrews’ (1991) automatic bandwidth, it in general works well even when the underlying process is not Gaussian or AR(1). On the other hand, Wilhelm’s (2015) and Hirukawa’s (2010) bandwidths are very
close on average, and the latter appears to have the potential to match the former. In addition, Hirukawa’s (2010) SEPI bandwidth can be implemented under less stringent conditions than Wilhelm’s (2015). Therefore, natural extensions of Wilhelm’s(2015) approach would be: (i) relaxing the conditions for its implementation; (ii) investigating applicability of prewhitening (which is not supported theoretically); and (iii) improving its computation in combination with the SEPI algorithm, if
applicable.
The second simulation study indicates that while the parametric VARHAC estimator is easy to implement, estimation errors of AR coeﬃcients induce additional bias and/or variability in the estimates of cointegrating vectors. Hence, preference isstill given to nonparametric HAR estimators, whereas no estimation-optimal bandwidth choice method has been proposed in
the framework of semiparametric bias-corrected estimators of cointegrating regressions. It is worth tailoring the estimationoptimal bandwidth choice method to FMLS and CCR. As in Wilhelm (2015), this task may start from stochastic expansionsof the FMLS and CCR estimators. The contribution of prewhitening to bias reduction should be also justiﬁed theoretically.
Acknowledgment
The author would like to thank the Editor, an Associate Editor, two anonymous referees, Liudas Giraitis, Eiji Kurozumi,
Artem Prokhorov, Ken West, and the participants of ISF2021 for their constructive comments and suggestions, which helped
improve this article substantially. The author also gratefully acknowledges ﬁnancial support from Japan Society for the Promotion of Science (grant number 19K01595).
References
Abadir, K.M., Distaso, W., Giraitis, L., 2009. Two estimators of the long-run variance: Beyond short memory. Journal of Econometrics 150, 56–70.
Ahn, S.C., Gadarowski, C., Perez, M.F., 2012. Robust two-pass cross-sectional regressions: A minimum distance approach. Journal of Financial Econometrics
10, 669–701.
Akaike, H., 1973. Information theory and an extension of the maximum likelihood principle. In: Petrov, B.N., Csaki, F. (Eds.), Proceedings of the Second
International Symposium on Information Theory. Akademiai Kiado, Budapest, pp. 267–281.
Anderson, T.W., 1971. The Statistical Analysis of Time Series. John Wiley & Sons, New York.
Andrews, D.W.K., 1991. Heteroskedasticity and autocorrelation consistent covariance matrix estimation. Econometrica 59, 817–858.
Andrews, D.W.K., 1993. Tests for parameter instability and structural change with unknown change point. Econometrica 61, 821–856.
Andrews, D.W.K., Monahan, J.C., 1992. An improved heteroskedasticity and autocorrelation consistent covariance matrix estimator. Econometrica 60,
953–966.
Basu, S., Michailidis, G., 2015. Regularized estimation in sparse high-dimensional time series models. Annals of Statistics 43, 1535–1567.
Berk, K.N., 1974. Consistent autoregressive spectral estimates. Annals of Statistics 2, 489–502.
Blackman, R.B., Tukey, J.W., 1958. The measurement of power spectra from the point of view of communications engineering - Part I. Bell System Technical
Journal 37, 185–282.
Brillinger, D.R., 1975. Time Series: Data Analysis and Theory. Holt, Rinehart and Winston, New York.
Creel, M., Farell, M., 1996. SUR estimation of multiple time-series models with heteroscedasticity and serial correlation of unknown form. Economics Letters53, 239–245.
Cumby, R.E., Huizinga, J., Obstfeld, M., 1983. Two-step two-stage least squares estimation in models with rational expectations. Journal of Econometrics 21,
333–355.
Cushing, M.J., McGarvey, M.G., 1999. Covariance matrix estimation. In: Mátyás, L. (Ed.), Generalized Method of Moments Estimation. Cambridge University
Press, New York, pp. 63–95.
Dahlhaus, R., 2012. Locally stationary processes. In: Rao, T.S., Rao, S.S., Rao, C.R. (Eds.), Handbook of Statistics Volume 30: Time Series Analysis: Methods
and Applications. Elsevier, Amsterdam, pp. 351–413.
Davidson, J., De Jong, R.M., 2002. Consistency of kernel variance estimators for sums of semiparametric linear processes. Econometrics Journal 5, 160–175.
De Jong, R.M., Davidson, J., 20 0 0. Consistency of kernel estimators of heteroskedastic and autocorrelated covariance matrices. Econometrica 68, 407–423.
Den Haan, W.J., Levin, A.T., 1996. Inferences from parametric and non-parametric covariance matrix estimation procedures. Technical Working Paper 195.
National Bureau of Economic Research.
Den Haan, W.J., Levin, A.T., 1997. A practitioner’s guide to robust covariance matrix estimation. In: Maddala, G.S., Rao, C.R. (Eds.), Handbook of Statistics
Volume 15: Robust Inference. Elsevier, Amsterdam, pp. 299–342.
Den Haan, W.J., Levin, A.T., 20 0 0. Robust covariance matrix estimation with data-dependent VAR prewhitening order. Technical Working Paper 255. National
Bureau of Economic Research.
Diebold, F.X., Mariano, R.S., 1995. Comparing predictive accuracy. Journal of Business & Economic Statistics 13, 253–263.
Dou, L., 2020. Optimal HAR Inference. Unpublished Manuscript, Chinese University of Hong Kong.
Eichenbaum, M.S., Hansen, L.P., Singleton, K.J., 1988. A time series analysis of representative agent models of consumption and leisure choice under uncertainty. Quarterly Journal of Economics 103, 51–78.
Elliott, G., Rothenberg, T.J., Stock, J.H., 1996. Eﬃcient tests for an autoregressive unit root. Econometrica 64, 813–836.
Gallant, A.R., 1987. Nonlinear Statistical Models. John Wiley & Sons, New York.
Hall, A.R., 2005. Generalized Method of Moments. Oxford University Press, New York.
Hannan, E.J., 1957. The variance of the mean of a stationary process. Journal of the Royal Statistical Society, Series B 19, 282–285.
Hannan, E.J., 1970. Multiple Time Series. John Wiley & Sons, New York.
Hannan, E.J., Rissanen, J., 1982. Recursive estimation of mixed autoregressive-moving average order. Biometrika 69, 81–94.
Hansen, B.E., 1992a. Tests for parameter instability in regressions with I (1 ) processes. Journal of Business & Economic Statistics 10, 45–59.
Hansen, B.E., 1992b. Consistent covariance matrix estimation for dependent heterogenous processes. Econometrica 60, 967–972.
Hansen, L.P., 1982. Large sample properties of generalized method of moments estimators. Econometrica 50, 1029–1054.
Hansen, L.P., Hodrick, R.J., 1980. Forward exchange rates as optimal predictors of future spot rates: An econometric analysis. Journal of Political Economy
88, 829–853.
Hansen, L.P., Singleton, K.J., 1982. Generalized instrumental variables estimation of nonlinear rational expectations models. Econometrica 50, 1269–1286.
59

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

Hirukawa, M., 2006. A modiﬁed nonparametric prewhitened covariance estimator. Journal of Time Series Analysis 27, 441–476.
Hirukawa, M., 2010. A two-stage plug-in bandwidth selection and its implementation for covariance estimation. Econometric Theory 26, 710–743.
Hirukawa, M., 2011. How useful is yet another data-driven bandwidth in long-run variance estimation?: A simulation study on cointegrating regressions.
Economics Letters 111, 170–172.
Hodrick, R.J., 1992. Dividend yields and expected stock returns: Alternative procedures for inference and measurement. Review of Financial Studies 5,357–386.
Ibragimov, R., Müller, U.K., 2010. t-statistic based correlation and heterogeneity robust inference. Journal of Business & Economic Statistics 28, 453–468.
Jansson, M., 2002. Consistent covariance matrix estimation for linear processes. Econometric Theory 18, 1449–1459.
Jansson, M., 2004. The error in rejection probability of simple autocorrelation robust tests. Econometrica 72, 937–946.
Jones, M.C., Linton, O., Nielsen, J.P., 1995. A simple bias reduction method for density estimation. Biometrika 82, 327–338.
Jowett, G.H., 1955. The comparison of means of sets of observations from sections of independent stochastic series. Journal of the Royal Statistical Society,
Series B 17, 208–227.
Kawka, R., 2020. Convergence of spectral density estimators in the locally stationary framework. Econometrics and Statistics. forthcoming
Kejriwal, M., 2009. Tests for a mean shift with good size and monotonic power. Economics Letters 102, 78–82.
Kiefer, N.M., Vogelsang, T.J., 2002. Heteroskedasticity-autocorrelation robust standard errors using the Bartlett kernel without truncation. Econometrica 70,
2093–2095.
Kiefer, N.M., Vogelsang, T.J., 2005. A new asymptotic theory for heteroskedasticity-autocorrelation robust tests. Econometric Theory 21, 1130–1164.
Kiefer, N.M., Vogelsang, T.J., Bunzel, H., 20 0 0. Simple robust testing of regression hypotheses. Econometrica 68, 695–714.
Kitamura, Y., 1997. Empirical likelihood methods with weakly dependent processes. Annals of Statistics 25, 2084–2102.
Kitamura, Y., Stutzer, M., 1997. An information-theoretic alternative to generalized method of moments estimation. Econometrica 65, 861–874.
Kwiatkowski, D., Phillips, P.C.B., Schmidt, P., Shin, Y., 1992. Testing the null hypothesis of stationarity against the alternative of a unit root: How sure are
we that economic time series have a unit root? Journal of Econometrics 54, 159–178.
Lazarus, E., Lewis, D.J., Stock, J.H., 2021. Supplement to The Size-Power Tradeoff in HAR Inference (Econometrica, Vol.89, No.5, September 2021, 2497-2516).
Econometrica Supplementary Material.
Lazarus, E., Lewis, D.J., Stock, J.H., Watson, M.W., 2018. HAR inference: Recommendations for practice. Journal of Business & Economic Statistics 36, 541–559.
Lee, C.C., Phillips, P.C.B., 1994. An ARMA prewhitened long-run variance estimator. Unpublished Manuscript, Yale University.
Mark, N.C., 2001. International Macroeconomics and Finance: Theory and Econometric Methods. Blackwell, Malden, MA.
McCracken, M.W., 20 0 0. Robust out-of-sample inference. Journal of Econometrics 99, 195–223.
Müller, U.K., 2007. A theory of robust long-run variance estimation. Journal of Econometrics 141, 1331–1352.
Müller, U.K., 2014. HAC corrections for strongly autocorrelated time series. Journal of Business & Economic Statistics 32, 311–322.
Newey, W.K., West, K.D., 1987. A simple, positive semi-deﬁnite, heteroskedasticity and autocorrelation consistent covariance matrix. Econometrica 55,
703–708.
Newey, W.K., West, K.D., 1994. Automatic lag selection in covariance matrix estimation. Review of Economic Studies 61, 631–653.
Ng, S., Perron, P., 2001. Lag length selection and the construction of unit root tests with good size and power. Econometrica 69, 1519–1554.
Park, J.Y., 1992. Canonical cointegrating regressions. Econometrica 60, 119–143.
Park, J.Y., Ogaki, M., 1991. Inference in cointegrated models using VAR prewhitening to estimate shortrun dynamics. Rochester Center for Economic ResearchWorking Paper No. 281. University of Rochester.
Parzen, E., 1957. On consistent estimates of the spectrum of a stationary time series. Annals of Mathematical Statistics 28, 329–348.
Parzen, E., 1961. Mathematical considerations in the estimation of spectra. Technometrics 3, 167–190.
Parzen, E., 1963. Notes on Fourier analysis and spectrum windows. Technical Report No. 48. Department of Statistics, Stanford University. (Published inParzen, E. (1967): Time Series Analysis Papers. Holden-Day, San Francisco).
Perron, P., Ng, S., 1996. Useful modiﬁcations to some unit root tests with dependent errors and their local asymptotic properties. Review of Economic
Studies 63, 435–463.
Phillips, P.C.B., 1987. Time series regression with a unit root. Econometrica 55, 277–301.
Phillips, P.C.B., 1991. Spectral regression for cointegrated time series. In: Barnett, W.A., Powell, J., Tauchen, G.E. (Eds.), Nonparametric and Semiparametric
Methods in Econometrics and Statistics: Proceedings of the Fifth International Symposium in Economic Theory and Econometrics. Cambridge University
Press, Cambridge, UK, pp. 413–435.
Phillips, P.C.B., 2005. HAC estimation by automated regression. Econometric Theory 21, 116–142.
Phillips, P.C.B., Hansen, B.E., 1990. Statistical inference in instrumental variables regression with I (1 ) processes. Review of Economic Studies 57, 99–125.
Phillips, P.C.B., Ouliaris, S., 1990. Asymptotic properties of residual based tests for cointegration. Econometrica 58, 165–193.
Phillips, P.C.B., Perron, P., 1988. Testing for a unit root in time series regression. Biometrika 75, 335–346.
Phillips, P.C.B., Sun, Y., Jin, S., 2006. Spectral density estimation and robust hypothesis testing using steep origin kernels without truncation. International
Economic Review 47, 837–894.
Phillips, P.C.B., Sun, Y., Jin, S., 2007. Long run variance estimation and robust regression testing using sharp origin kernels with no truncation. Journal of
Statistical Planning and Inference 137, 985–1023.
Politis, D.N., 2011. Higher-order accurate, positive semideﬁnite estimation of large-sample covariance and spectral density matrices. Econometric Theory 27,
703–744.
Press, H., Tukey, J.W., 1956. Power spectral methods of analysis and their application to problems in airplane dynamics. In: Flight Test Manual, NATO
Advisory Group for Aeronautical Research and Development, IV-C, pp. 1-41. (reprinted as Bell System Monograph No. 2606.).
Priestley, M.B., 1981. Spectral Analysis and Time Series. Academic Press, San Diego.
Robinson, P.M., 1998. Inference-without-smoothing in the presence of nonparametric autocorrelation. Econometrica 66, 1163–1182.
Robinson, P.M., 2005. Robust covariance matrix estimation: HAC estimates with long memory/antipersistence correction. Econometric Theory 21, 171–180.
Rossi, B., Inoue, A., 2012. Out-of-sample forecast tests robust to the choice of window size. Journal of Business & Economic Statistics 30, 432–453.
Schwarz, G., 1978. Estimating the dimension of a model. Annals of Statistics 6, 461–464.
Sheather, S.J., Jones, M.C., 1991. A reliable data-based bandwidth selection method for kernel density estimation. Journal of the Royal Statistical Society,
Series B 53, 683–690.
Silverman, B.W., 1986. Density Estimation for Statistics and Data Analysis. Chapman & Hall, London.
Smith, R.J., 1997. Alternative semi-parametric likelihood approaches to generalised method of moments estimation. Econometrics Journal 107, 503–519.
Smith, R.J., 2005. Automatic positive semideﬁnite HAC covariance matrix and GMM estimation. Econometric Theory 21, 158–170.
Smith, R.J., 2011. GEL criteria for moment condition models. Econometric Theory 27, 1192–1235.
Stuetzle, W., Mittal, Y., 1979. Some comments on the asymptotic behavior of robust smoothers. In: Gasser, T., Rosenblatt, M. (Eds.), Smoothing Techniques
for Curve Estimation: Proceedings of a Workshop Held in Heidelberg, April 2-4, 1979. Springer-Verlag, Berlin, pp. 191–195.
Sul, D., Phillips, P.C.B., Choi, C.Y., 2005. Prewhitening bias in HAC estimation. Oxford Bulletin of Economics and Statistics 67, 517–546.
Sun, Y., 2011. Robust trend inference with series variance estimator and testing-optimal smoothing parameter. Journal of Econometrics 164, 345–366.
Sun, Y., 2014a. Let’s ﬁx it: Fixed-b asymptotics versus small-b asymptotics in heteroskedasticity and autocorrelation robust inference. Journal of Econometrics
178, 659–677.
Sun, Y., 2014b. Fixed-smoothing asymptotics in a two-step generalized method of moments framework. Econometrica 82, 2327–2370.
Sun, Y., Phillips, P.C.B., Jin, S., 2008. Optimal bandwidth selection in heteroskedasticity-autocorrelation robust testing. Econometrica 76, 175–194.

60

M. Hirukawa

Econometrics and Statistics 27 (2023) 36–61

Sun, Y., Phillips, P.C.B., Jin, S., 2011. Power maximization and size control in heteroskedasticity and autocorrelation robust tests with exponentiated kernels.
Econometric Theory 27, 1320–1368.
Sun, Y., Yang, J., 2020. Testing-optimal kernel choice in HAR inference. Journal of Econometrics 219, 123–136.
Vogelsang, T.J., Wagner, M., 2013. A ﬁxed-b perspective on the Phillips-Perron unit root tests. Econometric Theory 29, 609–628.
Vogelsang, T.J., Wagner, M., 2014. Integrated modiﬁed OLS estimation and ﬁxed-b inference for cointegrating regressions. Journal of Econometrics 178,
741–760.
Wang, Q., Wu, N., 2012. Long-run covariance and its applications in cointegrating regression. Stata Journal 12, 515–542.
Wenger, K., Leschinski, C., 2021. Fixed-bandwidth CUSUM tests under long memory. Econometrics and Statistics 20, 46–61.
West, K.D., 1996. Asymptotic inference about predictive ability. Econometrica 64, 1067–1084.
West, K.D., 1997. Another heteroskedasticity- and autocorrelation-consistent covariance matrix estimator. Journal of Econometrics 76, 171–191.
West, K.D., 2008. Heteroskedasticity and autocorrelation corrections. In: Durlauf, S.N., Brume, L.E. (Eds.), The New Palgrave Dictionary of Economics, Second
Edition, Volume 4. Palgrave Macmillan, Hampshire, UK, pp. 6–12.
White, H., Domowitz, I., 1984. Nonlinear regression with dependent observations. Econometrica 52, 143–162.
Wilhelm, D., 2015. Optimal bandwidth selection for robust generalized method of moments estimation. Econometric Theory 31, 1054–1077.
Xiao, Z., Linton, O., 2002. A nonparametric prewhitened covariance estimator. Journal of Time Series Analysis 23, 215–250.
Ye, X., Sun, Y., 2018. Heteroscedasticity and autocorrelation robust F and t tests in Stata. Stata Journal 18, 951–980.

61