1## vcovHAC() is the general workhorse for HAC estimation 2## although the essential part is now moved to meatHAC() 3## interfacing the sandwich() function 4 5vcovHAC <- function(x, ...) { 6 UseMethod("vcovHAC") 7} 8 9vcovHAC.default <- function(x, order.by = NULL, prewhite = FALSE, 10 weights = weightsAndrews, adjust = TRUE, diagnostics = FALSE, 11 sandwich = TRUE, ar.method = "ols", data = list(), ...) 12{ 13 rval <- meatHAC(x, order.by = order.by, prewhite = prewhite, 14 weights = weights, adjust = adjust, diagnostics = diagnostics, 15 ar.method = ar.method, data = data) 16 17 if(sandwich) { 18 diagn <- attr(rval, "diagnostics") 19 rval <- sandwich(x, meat. = rval, ...) 20 attr(rval, "diagnostics") <- diagn 21 } 22 23 return(rval) 24} 25 26meatHAC <- function(x, order.by = NULL, prewhite = FALSE, 27 weights = weightsAndrews, adjust = TRUE, diagnostics = FALSE, 28 ar.method = "ols", data = list(), ...) 29{ 30 ## ensure that NAs are omitted 31 if(is.list(x) && !is.null(x$na.action)) class(x$na.action) <- "omit" 32 33 prewhite <- as.integer(prewhite) 34 35 umat <- estfun(x, ...)[, , drop = FALSE] 36 if(is.zoo(umat)) umat <- as.matrix(coredata(umat)) 37 n.orig <- n <- nrow(umat) 38 k <- ncol(umat) 39 40 if(!is.null(order.by)) 41 { 42 if(inherits(order.by, "formula")) { 43 z <- model.matrix(order.by, data = data) 44 z <- as.vector(z[,ncol(z)]) 45 } else { 46 z <- order.by 47 } 48 index <- order(z) 49 } else { 50 index <- 1:n 51 } 52 umat <- umat[index, , drop = FALSE] 53 54 if(prewhite > 0) { 55 var.fit <- try(ar(umat, order.max = prewhite, demean = FALSE, aic = FALSE, method = ar.method)) 56 if(inherits(var.fit, "try-error")) stop(sprintf("VAR(%i) prewhitening of estimating functions failed", prewhite)) 57 if(k > 1) D <- solve(diag(ncol(umat)) - apply(var.fit$ar, 2:3, sum)) 58 else D <- as.matrix(1/(1 - sum(var.fit$ar))) 59 umat <- as.matrix(na.omit(var.fit$resid)) 60 n <- n - prewhite 61 } 62 63 if(is.function(weights)) 64 weights <- weights(x, order.by = order.by, prewhite = prewhite, ar.method = ar.method, data = data) 65 66 if(length(weights) > n) { 67 warning("more weights than observations, only first n used") 68 weights <- weights[1:n] 69 } 70 71 utu <- 0.5 * crossprod(umat) * weights[1] 72 wsum<-n*weights[1]/2 73 w2sum<-n*weights[1]^2/2 74 75 if(length(weights) > 1) { 76 for (ii in 2:length(weights)) { 77 utu <- utu + weights[ii] * crossprod(umat[1:(n-ii+1),,drop=FALSE], umat[ii:n,,drop=FALSE]) 78 wsum <- wsum + (n-ii+1) * weights[ii] 79 w2sum <- w2sum + (n-ii+1) * weights[ii]^2 80 } 81 } 82 83 utu <- utu + t(utu) 84 85 ## Andrews: multiply with df n/(n-k) 86 if(adjust) utu <- n.orig/(n.orig-k) * utu 87 88 if(prewhite > 0) { 89 utu <- crossprod(t(D), utu) %*% t(D) 90 } 91 92 wsum <- 2*wsum 93 w2sum <- 2*w2sum 94 bc <- n^2/(n^2 - wsum) 95 df <- n^2/w2sum 96 97 rval <- utu/n.orig 98 99 if(diagnostics) attr(rval, "diagnostics") <- list(bias.correction = bc, df = df) 100 return(rval) 101} 102 103 104## weightsAndrews() and bwAndrews() implement the HAC estimation 105## procedure described in Andrews (1991) and Andrews & Monahan (1992) 106## kernHAC() is the convenience interface. 107## (Note, that bwNeweyWest() can also be used with weightsAndrews()) 108 109weightsAndrews <- function(x, order.by = NULL, bw = bwAndrews, 110 kernel = c("Quadratic Spectral", "Truncated", "Bartlett", "Parzen", "Tukey-Hanning"), 111 prewhite = 1, ar.method = "ols", tol = 1e-7, data = list(), verbose = FALSE, ...) 112{ 113 ## ensure that NAs are omitted 114 if(is.list(x) && !is.null(x$na.action)) class(x$na.action) <- "omit" 115 116 kernel <- match.arg(kernel) 117 if(is.function(bw)) 118 bw <- bw(x, order.by = order.by, kernel = kernel, 119 prewhite = prewhite, data = data, ar.method = ar.method, ...) 120 if(verbose) cat(paste("\nBandwidth chosen:", format(bw), "\n")) 121 122 n <- NROW(estfun(x)) - as.integer(prewhite) 123 124 weights <- kweights(0:(n-1)/bw, kernel = kernel) 125 weights <- weights[1:max(which(abs(weights) > tol))] 126 return(weights) 127} 128 129bwAndrews <- function(x, order.by = NULL, kernel = c("Quadratic Spectral", "Truncated", 130 "Bartlett", "Parzen", "Tukey-Hanning"), approx = c("AR(1)", "ARMA(1,1)"), 131 weights = NULL, prewhite = 1, ar.method = "ols", data = list(), ...) 132{ 133 ## ensure that NAs are omitted 134 if(is.list(x) && !is.null(x$na.action)) class(x$na.action) <- "omit" 135 136 kernel <- match.arg(kernel) 137 approx <- match.arg(approx) 138 prewhite <- as.integer(prewhite) 139 140 umat <- estfun(x)[,, drop = FALSE] 141 if(is.zoo(umat)) umat <- as.matrix(coredata(umat)) 142 n <- nrow(umat) 143 k <- ncol(umat) 144 145 if(!is.null(order.by)) 146 { 147 if(inherits(order.by, "formula")) { 148 z <- model.matrix(order.by, data = data) 149 z <- as.vector(z[,ncol(z)]) 150 } else { 151 z <- order.by 152 } 153 index <- order(z) 154 } else { 155 index <- 1:n 156 } 157 158 umat <- umat[index, , drop = FALSE] 159 160 ## compute weights (try to set the intercept weight to 0) 161 #### could be ignored by using: weights = 1 162 163 if(is.null(weights)) { 164 weights <- rep(1, k) 165 unames <- colnames(umat) 166 if(!is.null(unames) && "(Intercept)" %in% unames) 167 weights[which(unames == "(Intercept)")] <- 0 168 else { 169 res <- try(as.vector(rowMeans(estfun(x)/model.matrix(x), na.rm = TRUE)), silent = TRUE) 170 if(inherits(res, "try-error")) res <- try(residuals(x), silent = TRUE) 171 if(!inherits(res, "try-error")) weights[which(colSums((umat - res)^2) < 1e-16)] <- 0 172 } 173 if(isTRUE(all.equal(weights, rep(0, k)))) weights <- rep(1, k) 174 } else { 175 weights <- rep(weights, length.out = k) 176 } 177 if(length(weights) < 2) weights <- 1 178 179 if(prewhite > 0) { 180 var.fit <- ar(umat, order.max = prewhite, demean = FALSE, aic = FALSE, method = ar.method) 181 if(inherits(var.fit, "try-error")) stop(sprintf("VAR(%i) prewhitening of estimating functions failed", prewhite)) 182 umat <- as.matrix(na.omit(var.fit$resid)) 183 n <- n - prewhite ##?? 184 } 185 186 if(approx == "AR(1)") { 187 fitAR1 <- function(x) { 188 rval <- ar(x, order.max = 1, aic = FALSE, method = "ols") 189 rval <- c(rval$ar, sqrt(rval$var.pred)) 190 names(rval) <- c("rho", "sigma") 191 return(rval) 192 } 193 194 ar.coef <- apply(umat, 2, fitAR1) 195 196 denum <- sum(weights * (ar.coef["sigma",]/(1-ar.coef["rho",]))^4) 197 alpha2 <- sum(weights * 4 * ar.coef["rho",]^2 * ar.coef["sigma",]^4/(1-ar.coef["rho",])^8) / denum 198 alpha1 <- sum(weights * 4 * ar.coef["rho",]^2 * ar.coef["sigma",]^4/((1-ar.coef["rho",])^6 * (1+ar.coef["rho",])^2)) / denum 199 200 } else { 201 202 fitARMA11 <- function(x) { 203 rval <- arima(x, order = c(1, 0, 1), include.mean = FALSE) 204 rval <- c(rval$coef, sqrt(rval$sigma2)) 205 names(rval) <- c("rho", "psi", "sigma") 206 return(rval) 207 } 208 209 arma.coef <- apply(umat, 2, fitARMA11) 210 211 denum <- sum(weights * ((1 + arma.coef["psi",]) * arma.coef["sigma",]/(1-arma.coef["rho",]))^4) 212 alpha2 <- sum(weights * 4 * ((1 + arma.coef["rho",] * arma.coef["psi",]) * ( 213 arma.coef["rho",] + arma.coef["psi",]))^2 * arma.coef["sigma",]^4/ 214 (1-arma.coef["rho",])^8) / denum 215 alpha1 <- sum(weights * 4 * ((1 + arma.coef["rho",] * arma.coef["psi",]) * ( 216 arma.coef["rho",] + arma.coef["psi",]))^2 * arma.coef["sigma",]^4/ 217 ((1-arma.coef["rho",])^6 * (1+arma.coef["rho",])^2)) / denum 218 } 219 220 rval <- switch(kernel, 221 "Truncated" = {0.6611 * (n * alpha2)^(1/5)}, 222 "Bartlett" = {1.1447 * (n * alpha1)^(1/3)}, 223 "Parzen" = {2.6614 * (n * alpha2)^(1/5)}, 224 "Tukey-Hanning" = {1.7462 * (n * alpha2)^(1/5)}, 225 "Quadratic Spectral" = {1.3221 * (n * alpha2)^(1/5)}) 226 227 return(rval) 228} 229 230kernHAC <- function(x, order.by = NULL, prewhite = 1, bw = bwAndrews, 231 kernel = c("Quadratic Spectral", "Truncated", "Bartlett", "Parzen", "Tukey-Hanning"), 232 approx = c("AR(1)", "ARMA(1,1)"), adjust = TRUE, diagnostics = FALSE, sandwich = TRUE, 233 ar.method = "ols", tol = 1e-7, data = list(), verbose = FALSE, ...) 234{ 235 myweights <- function(x, order.by = NULL, prewhite = FALSE, ar.method = "ols", data = list()) 236 weightsAndrews(x, order.by = order.by, prewhite = prewhite, bw = bw, 237 kernel = kernel, approx = approx, ar.method = ar.method, tol = tol, 238 data = data, verbose = verbose, ...) 239 vcovHAC(x, order.by = order.by, prewhite = prewhite, 240 weights = myweights, adjust = adjust, diagnostics = diagnostics, 241 sandwich = sandwich, ar.method = ar.method, data = data) 242} 243 244 245 246## weightsLumley() implements the WEAVE estimators from 247## Lumley & Heagerty (1999) 248## weave() is a convenience interface 249 250weightsLumley <- function(x, order.by = NULL, C = NULL, 251 method = c("truncate", "smooth"), acf = isoacf, tol = 1e-7, data = list(), ...) 252{ 253 ## ensure that NAs are omitted 254 if(is.list(x) && !is.null(x$na.action)) class(x$na.action) <- "omit" 255 256 method <- match.arg(method) 257 res <- residuals(x, "response") #FIXME# available for which models? 258 n <- length(res) 259 260 if(!is.null(order.by)) 261 { 262 if(inherits(order.by, "formula")) { 263 z <- model.matrix(order.by, data = data) 264 z <- as.vector(z[,ncol(z)]) 265 } else { 266 z <- order.by 267 } 268 index <- order(z) 269 } else { 270 index <- 1:n 271 } 272 res <- res[index] 273 274 rhohat <- acf(res) 275 276 switch(method, 277 "truncate" = { 278 if(is.null(C)) C <- 4 279 lag <- max((1:length(rhohat))[rhohat^2*n > C]) 280 weights <- rep(1, lag) 281 }, 282 "smooth" = { 283 if(is.null(C)) C <- 1 284 weights <- C * n * rhohat^2 285 weights <- ifelse(weights > 1, 1, weights) 286 weights <- weights[1:max(which(abs(weights) > tol))] 287 }) 288 289 return(weights) 290} 291 292 293 294weave <- function(x, order.by = NULL, prewhite = FALSE, C = NULL, 295 method = c("truncate", "smooth"), acf = isoacf, adjust = FALSE, 296 diagnostics = FALSE, sandwich = TRUE, tol = 1e-7, data = list(), ...) 297{ 298 myweights <- function(x, order.by = NULL, prewhite = FALSE, data = list(), ...) 299 weightsLumley(x, order.by = order.by, prewhite = prewhite, C = C, 300 method = method, acf = acf, tol = tol, data = data) 301 vcovHAC(x, order.by = order.by, prewhite = prewhite, 302 weights = myweights, adjust = adjust, diagnostics = diagnostics, 303 sandwich = sandwich, data = data) 304} 305 306 307## bwNeweyWest() implements the procedure from Newey & West (1994) 308## It works for Bartlett/Parzen/QS kernels and can thus be passed 309## to weightsAndrews() and kernHAC() respectively. 310## A convenience interface NeweyWest() to only the Bartlett kernel 311## is also available. 312 313bwNeweyWest <- function(x, order.by = NULL, kernel = c("Bartlett", "Parzen", 314 "Quadratic Spectral", "Truncated", "Tukey-Hanning"), weights = NULL, prewhite = 1, 315 ar.method = "ols", data = list(), ...) 316{ 317 ## ensure that NAs are omitted 318 if(is.list(x) && !is.null(x$na.action)) class(x$na.action) <- "omit" 319 320 kernel <- match.arg(kernel) 321 if(kernel %in% c("Truncated", "Tukey-Hanning")) 322 stop(paste("Automatic bandwidth selection only available for ", 323 dQuote("Bartlett"), ", ", dQuote("Parzen"), " and ", dQuote("Quadratic Spectral"), 324 " kernel. Use ", sQuote("bwAndrews"), " instead.", sep = "")) 325 prewhite <- as.integer(prewhite) 326 327 umat <- estfun(x)[,, drop = FALSE] 328 if(is.zoo(umat)) umat <- as.matrix(coredata(umat)) 329 n <- nrow(umat) 330 k <- ncol(umat) 331 332 if(!is.null(order.by)) 333 { 334 if(inherits(order.by, "formula")) { 335 z <- model.matrix(order.by, data = data) 336 z <- as.vector(z[,ncol(z)]) 337 } else { 338 z <- order.by 339 } 340 index <- order(z) 341 } else { 342 index <- 1:n 343 } 344 345 umat <- umat[index, , drop = FALSE] 346 347 ## compute weights (try to set the intercept weight to 0) 348 #### could be ignored by using: weights = 1 349 350 if(is.null(weights)) { 351 weights <- rep(1, k) 352 unames <- colnames(umat) 353 if(!is.null(unames) && "(Intercept)" %in% unames) 354 weights[which(unames == "(Intercept)")] <- 0 355 else { 356 res <- try(as.vector(rowMeans(estfun(x)/model.matrix(x), na.rm = TRUE)), silent = TRUE) 357 if(inherits(res, "try-error")) res <- try(residuals(x), silent = TRUE) 358 if(!inherits(res, "try-error")) weights[which(colSums((umat - res)^2) < 1e-16)] <- 0 359 } 360 if(all(weights <= 0)) weights <- rep(1, length.out = k) 361 } else { 362 weights <- rep(weights, length.out = k) 363 } 364 if(length(weights) < 2) weights <- 1 365 366 ## select lag truncation according to Table II C. from Newey & West (1994) 367 mrate <- switch(kernel, 368 "Bartlett" = 2/9, 369 "Parzen" = 4/25, 370 "Quadratic Spectral" = 2/25) 371 m <- floor(ifelse(prewhite > 0, 3, 4) * (n/100)^mrate) 372 373 if(prewhite > 0) { 374 var.fit <- ar(umat, order.max = prewhite, demean = FALSE, aic = FALSE, method = ar.method) 375 if(inherits(var.fit, "try-error")) stop(sprintf("VAR(%i) prewhitening of estimating functions failed", prewhite)) 376 umat <- as.matrix(na.omit(var.fit$resid)) 377 n <- n - prewhite 378 } 379 380 ## compute weighted variances 381 hw <- umat %*% weights 382 sigmaj <- function(j) sum(hw[1:(n-j)] * hw[(j+1):n])/n 383 sigma <- sapply(0:m, sigmaj) 384 s0 <- sigma[1] + 2*sum(sigma[-1]) 385 s1 <- 2 * sum(1:m * sigma[-1]) 386 s2 <- 2 * sum((1:m)^2 * sigma[-1]) 387 388 ## use parameters as in Table I B. 389 ## choose 1/(2*q + 1) 390 qrate <- 1/(2 * ifelse(kernel == "Bartlett", 1, 2) + 1) 391 ## compute gamma 392 rval <- switch(kernel, 393 "Bartlett" = { 1.1447 * ((s1/s0)^2)^qrate }, 394 "Parzen" = { 2.6614 * ((s2/s0)^2)^qrate }, 395 "Quadratic Spectral" = { 1.3221 * ((s2/s0)^2)^qrate }) 396 ## compute bandwidth 397 rval <- rval * (n + prewhite)^qrate 398 399 ## rval is not truncated. This is done in NeweyWest(), 400 ## but bwNeweyWest() can also be used without truncation. 401 402 return(rval) 403} 404 405NeweyWest <- function(x, lag = NULL, 406 order.by = NULL, prewhite = TRUE, adjust = FALSE, 407 diagnostics = FALSE, sandwich = TRUE, ar.method = "ols", data = list(), 408 verbose = FALSE) 409{ 410 if(is.null(lag)) lag <- floor(bwNeweyWest(x, 411 order.by = order.by, prewhite = prewhite, 412 ar.method = ar.method, data = data)) 413 if(verbose) cat(paste("\nLag truncation parameter chosen:", lag, "\n")) 414 415 myweights <- seq(1, 0, by = -(1/(lag + 1))) 416 vcovHAC(x, order.by = order.by, prewhite = prewhite, 417 weights = myweights, adjust = adjust, diagnostics = diagnostics, 418 sandwich = sandwich, ar.method = ar.method, data = data) 419} 420 421