  fe_ar1_pcse <- function(object, f=NULL, pV=NULL, tV=NULL, data=NULL, includeM4=FALSE, ... ) {
  # message("This is used to calculate panel-corrected standard errors for
  #         two-way fixed-effects panel data utilizing the prais-winsten
  #         AR1 correction.   When using both panel and time dummies,
  #         the original panelAR() summary function gives an error
  #         when it attempts to calculate the Wald Statistic, and therefore
  #         doesn't return the standard errors and other statistics.")
    source("myFunctions/summary_panelAR.R")
    `%||%` <- function(a, b) if (!is.null(a)) a else b
    
    data <- (data %||% df)
    pV <- (pV %||% "ID")
    tV <- (tV %||% "YR")  
    pS <- summary_panelAR(object)
    pSc <- pS$call
    f <- (f %||% eval(pSc$formula))
    
    
  results <- list()
  require(lmtest)
  require(panelAR)
  
  #METHOD #1.  Calculating manually using the vcov from the panelAR object
  ##This is taken directly from the summary.panelAR() function
 # if(class(object)== "panelAR") {
    rdf <- object$df.residual
    se <- sqrt(diag(object$vcov))
    coef <- object$coefficients
    t.stat <- (coef)/se
    p.val <- 2 * pt(abs(t.stat), rdf, lower.tail = FALSE)
    
    dfP <- as.data.frame(cbind(coef, se, t.stat, p.val))
    names(dfP) <- c("coefficients", "pcses", "tvals", "pvals")
    results$method1 <- dfP
    #} else {results$method1 <- NA}
  
  ##METHOD #2.  Use the modified summary.panelAR() function with Wald stat removed.
  ##Object needs to be of class
  
  #if(class(object)== "panelAR") {
    p <- summary_panelAR(object)$coefficients
    coefs <- p[,1]
    ses <- p[,2]
    tvals <- p[,3]
    pvals <- p[,4]
    dfP <- as.data.frame(cbind(coefs, ses, tvals, pvals))
    names(dfP) <- c("coefficients", "pcses", "tvals", "pvals")
    results$method2 <- dfP
   # } else {results$method2 <- NA}
  
  #METHOD #3.  Use function coeftest() from lmtest package.
  vc <- object$vcov
  coefM <- coeftest(object, vcov.=vc)
  coefs <- coefM[,1]
  ses <- coefM[,2]  ##this works!
  tvals <- coefM[,3]
  pvals <- coefM[,4]
  
  dfP <- as.data.frame(cbind(coefs, ses, tvals, pvals))
  names(dfP) <- c("coefficients", "pcses", "tvals", "pvals")
  results$method3 <- dfP
  
  ##Method #4.  Use modified panelAR() function, which returns an object of class "lm"
  ##  which then allows you to use pcse::vcovPC() function.
if(includeM4){
  source("myFunctions/pcse_revised.R")

  #m <- (data %||% pSc$data)
  m <- data
  m[,pV] <- as.factor(m[,pV])
  p <- pcse_r(object, groupN=m[,pV], groupT=m[,tV], pairwise = TRUE)
  coefs <- p$b
  ses <- p$pcse
  tvals <- p$tstats
  pvals <- p$pval
  dfP <- as.data.frame(cbind(coefs, ses, tvals, pvals))
  names(dfP) <- c("coefficients", "pcses", "tvals", "pvals")
  results$method4 <- dfP
  } else {results$method4 <- NA}
  
  
  return(results)
}