
panelAR_r <- function (formula, data, panelVar, timeVar, 
                       autoCorr = c("ar1", "none", "psar1"), 
                      panelCorrMethod = c("none", "phet", "pcse",                                                                                               "pwls", "parks"), rhotype = "breg", bound.rho = FALSE, rho.na.rm = FALSE, 
          panel.weight = c("t-1", "t"), dof.correction = FALSE, complete.case = FALSE, 
          seq.times = FALSE, singular.ok = TRUE) 
{
  
  file.sources = list.files(pattern="*.R", path="myFunctions/panelAR/")
  fsources <- paste("myFunctions/panelAR/", file.sources, sep="")
  sapply(fsources,source)

  env.base <- environment()
  
  call <- match.call()
  ind <- match(c("formula", "data", "timeVar"), names(call), 
               nomatch = 0)
  if (ind[1] == 0) {
    stop("A formula argument is required.", call. = FALSE)
  }
  if (ind[2] == 0) {
    stop("A data argument is required.", call. = FALSE)
  }
  if (ind[3] == 0) {
    stop("You must specify time ID using timeVar.", call. = FALSE)
  }
  autoCorr.method <- match.arg(autoCorr)
  panelCorr.method <- match.arg(panelCorrMethod)
  pMethod <- switch(panelCorr.method, none = "OLS", phet = "OLS", 
                    pcse = "OLS", pwls = "GLS", parks = "GLS")
  rhotype <- match.arg(rhotype, c("breg", "freg", "dw", "theil-nagar", 
                                  "scorr", "theil"))
  panel.weight <- match.arg(panel.weight)
  if (!is.logical(rho.na.rm)) {
    stop("rho.na.rm should be logical argument.")
  }
  if (!(timeVar %in% colnames(data))) {
    stop("Please make sure that timeVar is in the data object.", 
         call. = FALSE)
  }
  time.vec <- data[, timeVar]
  if (any(is.na(time.vec))) {
    stop("You cannot have NA values for the time variable.", 
         call. = FALSE)
  }
  if (!is.integer(time.vec) & !is.numeric(time.vec)) {
    stop("The time variable must be defined as an integer.", 
         call. = FALSE)
  }
  if (!is.null(panelVar)) {
    if (!(panelVar %in% colnames(data))) {
      stop("Please make sure that panelVar is in the data object.", 
           call. = FALSE)
    }
    else {
      panel.vec <- as.character(as.vector(data[, panelVar]))
      if (any(is.na(panel.vec))) {
        stop("You cannot have NA values for the panel ID variable.", 
             call. = FALSE)
      }
    }
  }
  else {
    panel.vec <- rep("1", nrow(data))
    if (panelCorr.method != "none") 
      warning("Without specifying a panel ID variable, data is assumed to come from a single panel and variance is assumed to be homoskedastic within that panel. Panel heteroskedasticity and/or correlation is ignored.", 
              call. = FALSE)
  }
  if (!is.null(timeVar)) {
    if (any(by(time.vec, panel.vec, function(x) any(table(x) > 
                                                    1)))) {
      stop("You must specify unique times for each observation in a given panel.", 
           call. = FALSE)
    }
  }
  order.index <- order(panel.vec, time.vec)
  data <- data[order.index, ]
  panel.vec <- panel.vec[order.index]
  time.vec <- time.vec[order.index]
  lm.out <- lm(formula = formula, data = data, singular.ok = singular.ok)
  mterms <- lm.out$terms
  aliased <- is.na(coef(lm.out))
  X <- model.matrix(lm.out)[, !aliased]
  mf <- model.frame(lm.out)
  y <- model.response(mf)
  yX <- cbind(y, X)
  original.e <- residuals(lm.out)
  var.names <- colnames(X)
  N <- length(y)
  rank <- lm.out$rank
  rdf <- N - rank
  obs.dropped <- lm.out$na.action
  if (!is.null(obs.dropped)) {
    data <- data[-obs.dropped, ]
    panel.vec <- panel.vec[-obs.dropped]
    time.vec <- time.vec[-obs.dropped]
  }
  if (seq.times) {
    time.vec <- as.vector(unlist(by(data, panel.vec, function(x) 1:nrow(x))))
    data[, timeVar] <- time.vec
  }
  units <- sort(unique(panel.vec))
  times <- sort(unique(time.vec))
  N.units <- length(units)
  N.times <- length(times)
  N.avgperpanel <- N/N.units
  if (N.times < 2) {
    stop("More than two time periods required.", call. = FALSE)
  }
  if (panelCorr.method == "parks" & (N.units > N.times)) {
    stop("Cannot estimate Parks-Kmenta method because of singularity.")
  }
  NT <- N.units * N.times
  balanced <- ifelse(N == NT, TRUE, FALSE)
  obs.mat <- reshape(cbind(data[, c(panelVar, timeVar)], TRUE), 
                     timevar = timeVar, idvar = panelVar, direction = "wide", 
                     new.row.names = units)[, -1]
  col.order <- order(as.integer(gsub("TRUE.", "", colnames(obs.mat))))
  obs.mat <- obs.mat[, col.order]
  colnames(obs.mat) <- times
  obs.mat[is.na(obs.mat)] <- FALSE
  obs.mat <- as.matrix(obs.mat)
  N.runs.panel <- apply(obs.mat, MARGIN = 1, function(x) sum(rle(x)$values == 
                                                               TRUE))
  if (any(N.runs.panel > 1)) {
    message(paste("The following units have non-consecutive observations. Use runs.analysis() on output for additional details: ", 
                  paste(names(N.runs.panel[N.runs.panel > 1]), collapse = ", "), 
                  ".", sep = ""))
  }
  e.mat <- matrix(NA, nrow = ncol(obs.mat), ncol = nrow(obs.mat))
  e.mat[t(obs.mat) == TRUE] <- original.e
  pw.output <- prais.correct(method = autoCorr.method, env.base = env.base)
  transformed.resids <- pw.output$pw.lm$residuals
  model.mat.pw <- model.matrix(pw.output$pw.lm)
  obs.mat.pw <- t(obs.mat)
  if (!is.null(pw.output$pw.lm$na.action)) {
    panel.vec.pw <- panel.vec[-pw.output$pw.lm$na.action]
    time.vec.pw <- time.vec[-pw.output$pw.lm$na.action]
    obs.mat.pw[obs.mat.pw == TRUE][pw.output$pw.lm$na.action] <- FALSE
  }
  else {
    panel.vec.pw <- panel.vec
    time.vec.pw <- time.vec
  }
  if (panelCorr.method == "none") {
    sigma <- mean(transformed.resids^2)
    Sigma <- diag(sigma, N.units)
    Omega <- diag(sigma, nrow(model.mat.pw))
    N.cov <- 1
    res <- switch(pMethod, OLS = ols(env.base), GLS = gls(env.base))
  }
  else if (panelCorr.method == "phet" | panelCorr.method == 
           "pwls") {
    sigma.vec <- as.vector(by(transformed.resids, panel.vec.pw, 
                              function(x) mean(x^2)))
    if (length(sigma.vec) > 1) {
      Sigma <- diag(sigma.vec)
    }
    else {
      Sigma <- sigma.vec
    }
    Omega <- diag(rep(sigma.vec, times = as.integer(table(panel.vec.pw))))
    N.cov <- length(unique(panel.vec.pw))
    res <- switch(pMethod, OLS = ols(env.base), GLS = gls(env.base))
  }
  else {
    if (balanced) {
      E <- matrix(transformed.resids, nrow = N.times, ncol = N.units, 
                  byrow = FALSE)
      E.E <- crossprod(E)
      weight.mat <- crossprod(obs.mat.pw)
    }
    else {
      E <- obs.mat.pw
      E[E == TRUE] <- transformed.resids
      E[E == FALSE] <- 0
      if (complete.case) {
        I.com.case <- apply(obs.mat.pw, MARGIN = 1, function(x) prod(x) == 
                              1)
        if (!any(I.com.case)) {
          stop("Unable to compute correlated SEs / PCSEs because there are no time periods in common across all units. Instead, consider setting complete.case=FALSE.", 
               call. = FALSE)
        }
        else {
          if (sum(I.com.case) < (0.5 * N.avgperpanel)) {
            warning(paste("The number of time periods used for the calculation of correlated SEs / PCSEs (", 
                          as.character(sum(I.com.case)), ") is less than half the average number of time periods per panel (", 
                          as.character(round(N.avgperpanel, digits = 2)), 
                          "). Consider setting complete.case=FALSE.", 
                          sep = ""), call. = FALSE)
          }
          E[!I.com.case, ] <- 0
          E.E <- crossprod(E)
          weight.mat <- matrix(data = sum(I.com.case), 
                               nrow = nrow(E.E), ncol = nrow(E.E))
        }
      }
      else {
        E.E <- crossprod(E)
        weight.mat <- crossprod(obs.mat.pw)
      }
    }
    Sigma <- E.E/weight.mat
    N.cov <- length(Sigma[lower.tri(Sigma, diag = TRUE)])
    Sigma <- replace(Sigma, is.na(Sigma), 0)
    Omega <- kronecker(Sigma, diag(1, N.times))
    if (!balanced) {
      Omega <- Omega[as.vector(obs.mat.pw), as.vector(obs.mat.pw)]
    }
    res <- switch(pMethod, OLS = ols(env.base), GLS = gls(env.base))
  }
  coef <- as.vector(res$coef)
  vcov <- res$vcov
  colnames(vcov) <- rownames(vcov) <- names(coef) <- var.names
  yhat <- as.vector(X %*% coef)
  names(yhat) <- row.names(X)
  resids <- y - yhat
  if (dof.correction) {
    vcov <- vcov * (N/(N - rank))
  }
  panelStructure <- list(obs.mat = obs.mat, rho = pw.output$pw.rho, 
                         Sigma = Sigma, N.cov = N.cov)
  if (autoCorr.method == "psar1") {
    names(panelStructure$rho) <- rownames(obs.mat)
  }
  if (pMethod == "OLS") {
    transform.y.vec <- model.response(model.frame(pw.output$pw.lm))
    r2 <- 1 - sum(transformed.resids^2)/sum((transform.y.vec - 
                                               mean(transform.y.vec))^2)
  }
  else {
    r2 <- NULL
  }
  mf[, panelVar] <- panel.vec
  mf[, timeVar] <- time.vec
  fit <- list(coefficients = coef, residuals = resids, fitted.values = yhat, 
              rank = rank, df.residual = rdf, call = call, terms = mterms, 
              model = mf, aliased = aliased, na.action = obs.dropped, 
              vcov = vcov, r2 = r2, panelStructure = panelStructure)
  #class(fit) <- "panelAR"
  class(fit) <- "lm"
  return(fit)
}