R/test.hetero.R

In sheeputech/econometrics: Functions for Econometrics

#' Breusch-Pagun Test
#' @author Kazumasa Hirata
#' @description This is function for Breusch-Pagun Test for heteroskedasticity.
#' @param model lm
#' @return results list
#' @importFrom dplyr select mutate %>%
#' @export
#'
test.hetero.BreuschPagun <- function(model) {
    # prepare data
    regressors <- attr(terms(model), "term.labels")
    data <- model.frame(model) %>% dplyr::select(regressors)

    # auxiliary regression
    aux_data <- dplyr::mutate(data, residsq = residuals(model)^2)
    aux_model <- lm(residsq ~ ., data = aux_data)
    aux_sum <- summary(aux_model)

    # prepare result list
    obs <- nrow(aux_data)
    rsq <- aux_sum$r.squared
    LM <- obs * rsq
    df <- length(regressors)
    p <- pchisq(LM, df, lower.tail = F)

    result <- matrix(c(obs, round(rsq, 6), round(LM, 6), df, round(p, 6)), ncol = 1)
    rownames(result) <- c("aux.Observations", "aux.R-squared", "result.LM", "result.Chi.df", "result.Chi.p")
    colnames(result) <- c("results")
    return(as.data.frame(result))
}


#' White Test (Original)
#' @author Kazumasa Hirata
#' @description This is function for White Test for heteroskedasticity (original version).
#' @param model lm
#' @return results list
#' @import stats utils
#' @importFrom dplyr select %>%
#' @export
#'
test.hetero.White <- function(model) {
    # prepare for auxiliary regression
    regressors <- attr(terms(model), "term.labels")
    aux_data <- dplyr::select(model.frame(model), regressors)
    aux_vars <- combn(x = regressors, m = 2)
    for (var in regressors) {
        # Avoid perfect multi-collinearity
        if (any(apply(aux_data, MARGIN = 2, function(x){all(x == aux_data[, var]^2)}) == T)) next()
        aux_vars <- cbind(aux_vars, c(var, var))
    }

    for (i in 1:ncol(aux_vars)) {
        col <- aux_vars[, i]
        colname <- sprintf("%s_%s", col[1], col[2])
        aux_data[colname] <- aux_data[, col[1]] * aux_data[, col[2]]
    }
    aux_data <- dplyr::mutate(aux_data, residsq = residuals(model) ^ 2)

    # perform auxiliary regression
    aux_model <- lm(residsq ~ ., data = aux_data)
    aux_sum <- summary(aux_model)

    # organize result
    obs <- nrow(aux_data)
    rsq <- aux_sum$r.squared
    LM <- obs * rsq
    df <- ncol(aux_data) - 1
    p <- pchisq(LM, df, lower.tail = F)

    result <- matrix(c(obs, round(rsq, 6), round(LM, 6), df, round(p, 6)), ncol = 1)
    rownames(result) <- c("aux.Observations", "aux.R-squared", "result.LM", "result.Chi.df", "result.Chi.p")
    colnames(result) <- c("results")
    return(as.data.frame(result))
}


#' White Test (Special)
#' @author Kazumasa Hirata
#' @description This is function for White Test for heteroskedasticity (special version).
#' @param model lm
#' @return results list
#' @import stats utils
#' @export
#'
test.hetero.White.sp <- function(model) {
    # auxiliary regression
    aux_data <- data.frame(
        residsq = residuals(model)^2,
        pred = predict(model),
        predsq = predict(model)^2
    )
    aux_model <- lm(residsq ~ pred + predsq, data = aux_data)
    aux_sum <- summary(aux_model)

    # organize result
    obs <- nrow(aux_data)
    rsq <- aux_sum$r.squared
    LM <- obs * rsq
    df <- 2
    p <- pchisq(LM, df, lower.tail = F)

    result <- matrix(c(obs, round(rsq, 6), round(LM, 6), df, round(p, 6)), ncol = 1)
    rownames(result) <- c("aux.Observations", "aux.R-squared", "result.LM", "result.Chi.df", "result.Chi.p")
    colnames(result) <- c("results")
    return(as.data.frame(result))
}