Econometric Computing with HC and HAC Covariance Matrix Estimators

This introduction to the R package sandwich is a (slightly) modiﬁed version of Zeileis (2004a), published in the Journal of Statistical Software . Data described by econometric models typically contains autocorrelation and/or het-eroskedasticity of unknown form and for inference in such models it is essential to use covariance matrix estimators that can consistently estimate the covariance of the model parameters. Hence, suitable heteroskedasticity-consistent (HC) and heteroskedasticity and autocorrelation consistent (HAC) estimators have been receiving attention in the econometric literature over the last 20 years. To apply these estimators in practice, an implementation is needed that preferably translates the conceptual properties of the underlying theoretical frameworks into computational tools. In this paper, such an implementation in the package sandwich in the R system for statistical computing is described and it is shown how the suggested functions provide reusable components that build on readily existing functionality and how they can be integrated easily into new inferential procedures or applications. The toolbox contained in sandwich is extremely ﬂexible and comprehensive, including speciﬁc functions for the most important HC and HAC estimators from the econometric literature. Several real-world data sets are used to illustrate how the functionality can be integrated into applications.


Introduction
This paper combines two topics that play an important role in applied econometrics: computational tools and robust covariance estimation.
Without the aid of statistical and econometric software modern data analysis would not be possible: hence, both practitioners and (applied) researchers rely on computational tools that should preferably implement state-of-the-art methodology and be numerically reliable, easy to use, flexible and extensible.
In many situations, economic data arises from time-series or cross-sectional studies which typically exhibit some form of autocorrelation and/or heteroskedasticity.If the covariance structure were known, it could be taken into account in a (parametric) model, but more often than not the form of autocorrelation and heteroskedasticity is unknown.In such cases, model parameters can typically still be estimated consistently using the usual estimating functions, but for valid inference in such models a consistent covariance matrix estimate is essential.Over the last 20 years several procedures for heteroskedasticity consistent (HC) and for heteroskedasticity and autocorrelation consistent (HAC) covariance estimation have been suggested in the econometrics literature (White 1980;MacKinnon and White 1985;Newey andWest 1987, 1994;Andrews 1991, among others) and are now routinely used in econometric analyses.
Many statistical and econometric software packages implement various HC and HAC estimators for certain inference procedures, so why is there a need for a paper about econometric computing with HC and HAC estimators?Typically, only certain special cases of such estimators-and not the general framework they are taken from-are implemented in statistical and econometric software packages and sometimes they are only available as options to certain inference functions.It is desirable to improve on this for two reasons: First, the literature suggested conceptual frameworks for HC and HAC estimation and it would only be natural to translate these conceptual properties into computational tools that reflect the flexibility of the general framework.Second, it is important, particularly for applied research, to have covariance matrices not only as options to certain tests but as stand-alone functions which can be used as modular building blocks and plugged into various inference procedures.This is becoming more and more relevant, because today, as Cribari-Neto and Zarkos (2003) point out, applied researchers typically cannot wait until a certain procedure becomes available in the software package of their choice but are often forced to program new techniques themselves.Thus, just as suitable covariance estimators are routinely plugged into formulas in theoretical work, programmers should be enabled to plug in implementations of such estimators in computational work.Hence, the aim of this paper is to present an econometric computing approach to HC and HAC estimation that provides reusable components that can be used as modular building blocks in implementing new inferential techniques and in applications.All functions described are available in the package sandwich implemented in the R system for statistical computing (R Development Core Team 2004) which is currently not the most popular environment for econometric computing but which is finding increasing attention among econometricians (Cribari-Neto and Zarkos 1999;Racine and Hyndman 2002).Both R itself and the sandwich package (as well as all other packages used in this paper) are available at no cost under the terms of the general public licence (GPL) from the comprehensive R archive network (CRAN, http://CRAN.R-project.org/).R has no built-in support for HC and HAC estimation and at the time we started writing sandwich there was only one package that implements HC (but not HAC) estimators (the car package Fox 2002) but which does not allow for as much flexibility as the tools presented here.sandwich provides the functions vcovHC and vcovHAC implementing general classes of HC and HAC estimators.The names of the functions are chosen to correspond to vcov, R's generic function for extracting covariance matrices from fitted model objects.Below, we focus on the general linear regression model estimated by ordinary least squares (OLS), which is typically fitted in R using the function lm from which the standard covariance matrix (assuming spherical errors) can be extracted by vcov.Using the tools from sandwich, HC and HAC covariances matrices can now be extracted from the same fitted models using vcovHC and vcovHAC.Due to the object orientation of R, these functions are not only limited to the linear regression model but can be easily extended to other models.The HAC estimators are already available for generalized linear models (fitted by glm) and robust regression (fitted by rlm in package MASS).Another important feature of R that is used repeatedly below is that functions are first-level objects-i.e., functions can take functions as arguments and return functionswhich is particularly useful for defining certain procedures for data-driven computations such as the definition of the structure of covariance matrices in HC estimation and weighting schemes for HAC estimation.The remainder of this paper is structured as follows: To fix notations, Section 2 describes the linear regression model used and motivates the following sections.Section 3 gives brief literature reviews and describes the conceptual frameworks for HC and HAC estimation respectively and then shows how the conceptual properties are turned into computational tools in sandwich.Section 4 provides some illustrations and applications of these tools before a summary is given in Section 5.More details about the R code used are provided in an appendix.

The linear regression model
To fix notations, we consider the linear regression model with dependent variable y i , k-dimensional regressor x i with coefficient vector β and error term u i .
In the usual matrix notation comprising all n observations this can be formulated as y = Xβ + u.
In the general linear model, it is typically assumed that the errors have zero mean and variance VAR[u] = Ω.Under suitable regularity conditions (see e.g., Greene 1993;White 2000), the coefficients β can be consistently estimated by OLS giving the well-known OLS estimator β with corresponding OLS residuals ûi : where I n is the n-dimensional identity matrix and H is usually called hat matrix.The estimates β are unbiased and asymptotically normal (White 2000).Their covariance matrix Ψ is usually denoted in one of the two following ways: where Φ = n −1 X ΩX is essentially the covariance matrix of the scores or estimating functions The estimating functions evaluated at the parameter estimates Vi = V i ( β) have then sum zero.
For inference in the linear regression model, it is essential to have a consistent estimator for Ψ.What kind of estimator should be used for Ψ depends on the assumptions about Ω: In the classical linear model independent and homoskedastic errors with variance σ 2 are assumed yielding Ω = σ 2 I n and Ψ = σ 2 (X X) −1 which can be consistently estimated by plugging in the usual OLS estimator σ2 = (n − k) −1 n i=1 û2 i .But if the independence and/or homoskedasticity assumption is violated, inference based on this estimator Ψconst = σ(X X) −1 will be biased.HC and HAC estimators tackle this problem by plugging an estimate Ω or Φ into (4) or (5) respectively which are consistent in the presence of heteroskedasticity and autocorrelation respectively.Such estimators and their implementation are described in the following section.

Dealing with heteroskedasticity
If it is assumed that the errors u i are independent but potentially heteroskedastic-a situation which typically arises with cross-sectional data-their covariance matrix Ω is diagonal but has nonconstant diagonal elements.Therefore, various HC estimators ΨHC have been suggested which are constructed by plugging an estimate of type Ω = diag(ω 1 , . . ., ω n ) into Equation (4).These estimators differ in their choice of the ω i , an overview of the most important cases is given in the following: where h i = H ii are the diagonal elements of the hat matrix, h is their mean and δ i = min{4, h i / h}.
The first equation above yields the standard estimator Ψconst for homoskedastic errors.All others produce different kinds of HC estimators.The estimator HC0 was suggested in the econometrics literature by White (1980) and is justified by asymptotic arguments.The estimators HC1, HC2 and HC3 were suggested by MacKinnon and White (1985) to improve the performance in small samples.
A more extensive study of small sample behaviour was carried out by Long and Ervin (2000) which arrive at the conclusion that HC3 provides the best performance in small samples as it gives less weight to influential observations.Recently, Cribari-Neto ( 2004) suggested the estimator HC4 to further improve small sample performance, especially in the presence of influential observations.
All of these HC estimators ΨHC have in common that they are determined by ω = (ω 1 , . . ., ω n ) which in turn can be computed based on the residuals û, the diagonal of the hat matrix h and the degrees of freedom n−k.To translate these conceptual properties of this class of HC estimators into a computational tool, a function is required which takes a fitted regression model and the diagonal elements ω as inputs and returns the corresponding ΨHC .In sandwich, this is implemented in the function vcovHC which takes the following arguments: vcovHC(lmobj, omega = NULL, type = "HC3", ...) The first argument lmobj is an object as returned by lm, R's standard function for fitting linear regression models.The argument omega can either be the vector ω or a function for data-driven computation of ω based on the residuals û, the diagonal of the hat matrix h and the residual degrees of freedom n−k.Thus, it has to be of the form omega(residuals, diaghat, df): e.g., for computing HC3 omega is set to function(residuals, diaghat, df) residuals^2/(1 -diaghat)^2.
As a convenience option, a type argument can be set to "const", "HC0" (or equivalently "HC"), "HC1", "HC2", "HC3" (the default) or "HC4" and then vcovHC uses the corresponding omega function.As soon as omega is specified by the user, type is ignored.
In summary, by specfying ω-either as a vector or as a function-vcovHC can compute arbitrary HC covariance matrix estimates from the class of estimators outlined above.In Section 4, it will be illustrated how this function can be used as a building block when doing inference in linear regression models.

Dealing with autocorrelation
If the error terms u i are not independent, Ω is not diagonal and without further specification of a parametic model for the type of dependence it is typically burdensome to estimate Ω directly.However, if the form of heteroskedasticity and autocorrelation is unknown, a solution to this problem is to estimate Φ instead which is essentially the covariance matrix of the estimating functions1 .This is what HAC estimators do: ΨHAC is computed by plugging an estimate Φ into Equation ( 5) with where w = (w 0 , . . ., w n−1 ) is a vector of weights.An additional finite sample adjustment can be applied by multiplication with n/(n − k).For many data structures, it is a reasonable assumption that the autocorrelations should decrease with increasing lag = |i − j|-otherwise β can typically not be estimated consistently by OLS-so that it is rather intuitive that the weights w should also decrease.Starting from White and Domowitz (1984) and Newey and West (1987), different choices for the vector of weights w have been suggested in the econometrics literature which have been placed by Andrews (1991) in a more general framework of choosing the weights by kernel functions with automatic bandwidth selection.Andrews and Monahan (1992) show that the bias of the estimators can be reduced by prewhitening the estimating functions Vi using a vector autoregression (VAR) of order p and applying the estimator in Equation ( 6) to the VAR(p) residuals subsequently.Lumley and Heagerty (1999) suggest an adaptive weighting scheme where the weights are chosen based on the estimated autocorrelations of the residuals û.
All the estimators mentioned above are of the form (6), i.e., a weighted sum of lagged products of the estimating functions corresponding to a fitted regression model.Therefore, a natural implementation for this class of HAC estimators is the following: vcovHAC(lmobj, weights, prewhite = FALSE, adjust = TRUE, sandwich = TRUE, order.by,ar.method, data) The most important arguments are again the fitted linear model2 lmobj-from which the estimating functions Vi can easily be extracted using the generic function estfun(lmobj)-and the argument weights which specifys w.The latter can be either the vector w directly or a function to compute it from lmobj. 3 The argument prewhite specifies wether prewhitening should be used or not4 and adjust determines wether a finite sample correction by multiplication with n/(n − k) should be made or not.By setting sandwich it can be controlled wether the full sandwich estimator ΨHAC or only the "meat" Φ/n of the sandwich should be returned.The remaining arguments are a bit more technical: order.byspecifies by which variable the data should be ordered (the default is that they are already ordered, as is natural with time series data), which ar.method should be used for fitting the VAR(p) model (the default is OLS) and data provides a data frame from which order.by can be taken (the default is the environment from which vcovHAC is called).5 As already pointed out above, all that is required for specifying an estimator ΨHAC is the appropriate vector of weights (or a function for data-driven computation of the weights).For the most important estimators from the literature mentioned above there are functions for computing the corresponding weights readily available in sandwich.They are all of the form weights(lmobj, order.by,prewhite, ar.method, data), i.e., functions that compute the weights depending on the fitted model object lmobj and the arguments order.by,prewhite, data which are only needed for ordering and prewhitening.The function weightsAndrews implements the class of weights of Andrews (1991) and weightsLumley implements the class of weights of Lumley and Heagerty (1999).Both functions have convenience interfaces: kernHAC calls vcovHAC with weightsAndrews (and different defaults for some parameters) and weave calls vcovHAC with weightsLumley.Finally, a third convenience interface to vcovHAC is available for computing the estimator(s) of Newey andWest (1987, 1994).
• Newey and West (1987) suggested to use linearly decaying weights where L is the maximum lag, all other weights are zero.This is implemented in the function NeweyWest(lmobj, lag = NULL, ...) where lag specifies L and ... are (here, and in the following) further arguments passed to other functions, detailed information is always available in the reference manual.If lag is set to NULL (the default) the non-parametric bandwidth selection procedure of Newey and West (1994) is used.This is also available in a stand-alone function bwNeweyWest, see also below.Newey and West (1987) in Equation ( 7) when the bandwidth B is set to L + 1.The kernel recommended by Andrews (1991) and probably most used in the literature is the quadratic spectral kernel which leads to the following weights: where z = 6π/5 • /B.The definitions for the remaining kernels can be found in Andrews (1991).All kernel weights mentioned above are available in weightsAndrews(lmobj, kernel, bw, ...) where kernel specifies one of the kernels via a character string ("Truncated", "Bartlett", "Parzen", "Tukey-Hanning" or "Quadratic Spectral") and bw the bandwidth either as a scalar or as a function.The automatic bandwidth selection described in Andrews (1991) via AR(1) or ARMA(1,1) approximations is implemented in a function bwAndrews which is set as the default in weightsAndrews.For the Bartlett, Parzen and quadratic spectral kernels, Newey and West (1994) suggested a different nonparametric bandwidth selection procedure, which is implemented in bwNeweyWest and which can also be passed to weightsAndrews.As the flexibility of this conceptual framework of estimators leads to a lot of knobs and switches in the computational tools, a convenience function kernHAC for kernel-based HAC estimation has been added to sandwich that calls vcovHAC based on weightsAndrews and bwAndrews with defaults as motivated by Andrews (1991) and Andrews and Monahan (1992): by default, it computes a quadratic spectral kernel HAC estimator with VAR(1) prewhitening and automatic bandwidth selection based on an AR(1) approximation.But of course, all the options described above can also be changed by the user when calling kernHAC.
• Lumley and Heagerty (1999) suggested a different approach for specifying the weights in (6) based on some estimate ˆ of the autocorrelation of the residuals ûi at lag 0 = 1, . . ., n − 1.They suggest either to use truncated weights w = I{n ˆ 2 > C} (where I(•) is the indicator function) or smoothed weights w = min{1, C n ˆ 2 }, where for both a suitable constant C has to be specified.Lumley and Heagerty (1999) suggest using a default of C = 4 and C = 1 for the truncated and smoothed weights respectively.Note, that the truncated weights are equivalent to the truncated kernel from the framework of Andrews (1991) but using a different method for computing the truncation lag.To ensure that the weights |w | are decreasing, the autocorrelations have to be decreasing for increasing lag which can be achieved by using isotonic regression methods.In sandwich, these two weighting schemes are implemented in a function weightsLumley with a convenience interface weave (which stands for weighted empirical adaptive variance estimators) which again sets up the weights and then calls vcovHAC.Its most important arguments are weave(lmobj, method, C, ...) where method can be either "truncate" or "smooth" and C is by default 4 or 1 respectively.
To sum up, vcovHAC provides a simple yet flexible interface for general HAC estimation as defined in Equation ( 6).Arbitrary weights can be supplied either as vectors or functions for data-driven computation of the weights.As the latter might easily become rather complex, in particular due to the automatic choice of bandwidth or lag truncation parameters, three strategies suggested in the literature are readily available in sandwich: First, the Bartlett kernel weights suggested by Newey andWest (1987, 1994) are used in NeweyWest which by default uses the bandwidth selection function bwNeweyWest.Second, the weighting scheme introduced by Andrews (1991) for kernelbased HAC estimation with automatic bandwidth selection is implemented in weightsAndrews and bwAndrews with corresponding convenience interface kernHAC.Third, the weighted empirical adaptive variance estimation scheme suggested by Lumley and Heagerty (1999) is available in weightsLumley with convenience interface weave.
It is illustrated in the following section how these functions can be easily used in applications.

Applications and illustrations
In econometric analyses, the practitioner is only seldom interested in the covariance matrix Ψ (or Ω or Φ) per se, but mainly wants to compute them to use them for inferential procedures.Therefore, it is important that the functions vcovHC and vcovHAC described in the previous section can be easily supplied to other procedures such that the user does not necessarily have to compute the variances in advance.A typical field of application for HC and HAC covariances are partial t or z tests for assessing whether a parameter β j is significantly different from zero.Exploiting the (asymptotic) normality of the estimates, these tests are based on the t ratio βj / Ψjj and either use the asymptotic normal distribution or the t distribution with n − k degrees of freedom for computing p values (White 2000).This procedure is available in the R package lmtest (Zeileis and Hothorn 2002) in the generic function coeftest which has a default method applicable to fitted "lm" objects.coeftest(lmobj, vcov = NULL, df = NULL, ...) where vcov specifies the covariances either as a matrix (corresponding to the covariance matrix estimate) or as a function computing it from lmobj (corresponding to the covariance matrix estimator).By default, it uses the vcov method which computes Ψconst assuming spherical errors.The df argument determines the degrees of freedom: if df is finite and positive, a t distribution with df degrees of freedom is used, otherwise a normal approximation is employed.The default is to set df to n − k.Inference based on HC and HAC estimators is illustrated in the following using three real-world data sets: testing coefficients in two models from Greene (1993) and a structural change problem from Bai and Perron (2003).To make the results exactly reproducible for the reader, the commands for the inferential procedures is given along with their output within the text.A full list of commands, including those which produce the figures in the text, are provided (without output) in the appendix along with the versions of R and the packages used.Before we start with the examples, the sandwich and lmtest package have to be loaded: > library(sandwich) > library(lmtest)

Testing coefficients in cross-sectional data
A quadratic regression model for per capita expenditures on public schools explained by per capita income in the United States in 1979 has been analyzed by Greene (1993) and re-analyzed in Cribari-Neto (2004).The corresponding cross-sectional data for the 51 US states is given in Table 14.1 in Greene (1993) and available in sandwich in the data frame PublicSchools which can be loaded by: > data(PublicSchools) > ps <-na.omit(PublicSchools)> ps$Income <-ps$Income * 1e-04 where the second line omits a missing value (NA) in Wisconsin and assigns the result to a new data frame ps and the third line transforms the income to be in USD 10, 000.The quadratic regression can now easily be fit using the function lm which fits linear regression models specified by a symbolic formula via OLS.
> fm.ps <-lm(Expenditure ~Income + I(Income^2), data = ps) The fitted "lm" object fm.ps now contains the regression of the variable Expenditure on the variable Income and its sqared value, both variables are taken from the data frame ps.The question in this data set is whether the quadratic term is really needed, i.e., whether the coefficient of I(Income^2) is significantly different from zero.The partial quasi-t tests (or z tests) for all coefficients can be computed using the function coeftest.Greene (1993) assesses the significance using the HC0 estimator of White (1980).
> coeftest(fm.ps,df = Inf, vcov = vcovHC(fm.ps,type = "HC0")) The vcov argument specifies the covariance matrix as a matrix (as opposed to a function) which is returned by vcovHC(fm.ps,type = "HC0").As df is set to infinity (Inf) a normal approximation is used for computing the p values which seem to suggest that the quadratic term might be weakly significant.In his analysis, Cribari-Neto (2004) uses his HC4 estimator (among others) giving the following result: > coeftest(fm.ps,df = Inf, vcov = vcovHC(fm.ps,type = "HC4")) The quadratic term is clearly non-significant.The reason for this result is depicted in Figure 2 which shows the data along with the fitted linear and quadratic model-the latter being obviously heavily influenced by a single outlier: Alaska.Thus, the improved performance of the HC4 as compared to the HC0 estimator is due to the correction for high leverage points.q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.6 0.7 0.8 0.9 1.0 1.1 300 400 500 600 700 800 per capita income per capita spending on public schools Alaska

Testing coefficients in time-series data
Greene (1993) also anayzes a time-series regression model based on robust covariance matrix estimates: his Table 15.1 provides data on the nominal gross national product (GNP), nominal gross private domestic investment, a price index and an interest rate which is used to formulate a model that explains real investment by real GNP and real interest.The corresponding transformed variables RealInv, RealGNP and RealInt are stored in the data frame Investment in sandwich which can be loaded by: > data(Investment) Subsequently, the fitted linear regression model is computed by: > fm.inv <-lm(RealInv ~RealGNP + RealInt, data = Investment) and the significance of the coefficients can again be assessed by partial z tests using coeftest.Greene (1993) uses the estimator of Newey and West (1987) without prewhitening and with lag L = 4 for this purpose which is here passed as a matrix (as opposed to a function) to coeftest.
> coeftest (fm.inv,df = Inf,vcov = NeweyWest(fm.inv,lag = 4 For illustration purposes, we show how a new function implementing a particular HAC estimator can be easily set up using the tools provided by sandwich.This is particularly helpful if the same estimator is to be applied several times in the course of an analysis.Suppose, we want to use a Parzen kernel with VAR(2) prewhitening, no finite sample adjustment and automatic bandwidth selection according to Newey and West (1994).First, we set up the function parzenHAC and then pass this function to coeftest.The three estimators leads to slightly different standard errors, but all tests agree that real GNP has a highly significant influence while the real interest rate has not.The data along with the fitted regression are depicted in Figure 3.

Testing and dating structural changes in the presence of heteroskedasticity and autocorrelation
To illustrate that the functionality provided by the covariance estimators implemented in sandwich cannot only be used in simple settings, such as partial quasi-t tests, but also for more complicated tasks, we employ the real interest time series analyzed by Bai and Perron (2003).This series contains changes in the mean (see Figure 4, right panel) which Bai and Perron (2003) detect using several structural change tests based on F statistics and date using a dynamic programming algorithm.As the visualization suggests, this series exhibits both heteroskedasticity and autocorrelation, hence Bai and Perron (2003) use a quadratic spectral kernel HAC estimator in their analysis.Here, we use the same dating procedure but assess the significance using an OLS-based CUSUM test (Ploberger and Krämer 1992)  RealGNP RealInt RealInv q q q q q q q q q q q q q q q q q q q If autocorrelation and heteroskedasticity are present in the data, a robust variance estimator should be used: if x i is equal to unity, this can simply be achieved by replacing σ2 with Φ or n Ψ respectively.Here, we use the quadratic spectral kernel HAC estimator of Andrews (1991) with VAR(1) prewhitening and automatic bandwidth selection based on an AR(1) approximation as implemented in the function kernHAC.The p values for the OLS-based CUSUM test can be computed from the distribution of the supremum of a Brownian bridge (see e.g., Ploberger and Krämer 1992).This and other methods for testing, dating and monitoring structural changes are implemented in the R package strucchange (Zeileis, Leisch, Hornik, and Kleiber 2002) which contains the function gefp for fitting and assessing fluctuation processes including OLS-based CUSUM processes (see Zeileis 2004b, for more details).
After loading the package and the data, > library(strucchange) > data(RealInt) the command > ocus <-gefp(RealInt ~1, fit = lm, vcov = kernHAC) fits the OLS-based CUSUM process for a regression on the mean (RealInt ~1), using the function lm and estimating the variance using the function kernHAC.The fitted OLS-based CUSUM process can then be visualized together with its 5% critical value (horizontal lines) by plot(scus) which leads to a similar plot as in the left panel of Figure 4 (see the appendix for more details).As the process crosses its boundary, there is a significant change in the mean, while the clear peak in the process conveys that there is at least one strong break in the early 1980s.A formal significance test can also be carried out by sctest(ocus) which leads to a highly significant p value of 0.0082.Similarly, the same quadratic spectral kernel HAC estimator could also be used for computing and visualizing the supF test of Andrews (1993), the code is provided in the appendix.
Finally, the breakpoints in this model along with their confidence intervals can be computed by: > bp <-breakpoints(RealInt ~1) > confint(bp, vcov = kernHAC) Confidence intervals for breakpoints of optimal 3-segment partition: Call: confint.breakpointsfull(object= bp, vcov.= kernHAC) Breakpoints at observation number: 2.5 % breakpoints 97.The dating algorithm breakpoints implements the procedure described in Bai and Perron (2003) and estimates the timing of the structural changes by OLS.Therefore, in this step no covariance matrix estimate is required, but for computing the confidence intervals using a consistent covariance matrix estimator is again essential.The confint method for computing confidence intervals takes again a vcov argument which has to be a function (and not a matrix) because it has to be applied to several segments of the data.By default, it computes the breakpoints for the minimum BIC partition which gives in this case two breaks. 6The fitted three-segment model along with the breakpoints and their confidence intervals is depicted in the right panel of Figure 4.

Summary
This paper briefly reviews a class of heteroskedasticity-consistent (HC) and a class of heteroskedasticity and autocorrelation consistent (HAC) covariance matrix estimators suggested in the econometric literature over the last 20 years and introduces unified computational tools that reflect the flexibility and the conceptual ideas of the underlying theoretical frameworks.Based on these general tools, a number of special cases of HC and HAC estimators is provided including the most popular in applied econometric research.All the functions suggested are implemented in the package sandwich in the R system for statistical computing and designed in such a way that they build on readily available model fitting functions and provide building blocks that can be easily integrated into other programs or applications.To achieve this flexibility, the object orientation mechanism of R and the fact that functions are first-level objects are of prime importance.

Figure 1 :
Figure 1: Kernel functions for kernel-based HAC estimation

Figure 2 :
Figure 2: Expenditure on public schools and income with fitted models

Figure 3 :
Figure 3: Investment equation data with fitted model

Figure 4 :
Figure 4: OLS-based CUSUM test (left) and fitted model (right) for real interest data test of coefficients of "lm" object 'fm.ps': test of coefficients of "lm" object 'fm.ps': Newey and West (1994)automatic bandwidth selection procedure ofNewey and West (1994)with prewhitening should be used, this can be passed as a function to coeftest.