diff --git a/DESCRIPTION b/DESCRIPTION
index 32b7fea2..62a6900d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: eventstudyr
 Title: Estimation and Visualization of Linear Panel Event Studies
-Version: 1.1.5
+Version: 1.2.0
 Authors@R: 
     c(person(given = "Simon",
            family = "Freyaldenhoven",
@@ -43,6 +43,7 @@ Imports:
     data.table,
     dplyr,
     estimatr,
+    fixest,
     ggplot2,
     MASS,
     rlang,
@@ -53,7 +54,7 @@ VignetteBuilder:
     knitr
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Suggests: 
     rmarkdown,
     knitr,
diff --git a/NAMESPACE b/NAMESPACE
index 7867b6ac..d3bbc0ab 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -19,12 +19,20 @@ importFrom(data.table,setnames)
 importFrom(data.table,setorder)
 importFrom(data.table,setorderv)
 importFrom(data.table,shift)
+importFrom(dplyr,rename)
+importFrom(dplyr,select)
+importFrom(fixest,coeftable)
+importFrom(fixest,feols)
+importFrom(fixest,fitstat)
+importFrom(fixest,ssc)
 importFrom(pracma,inv)
 importFrom(pracma,pinv)
 importFrom(rlang,.data)
 importFrom(stats,as.formula)
+importFrom(stats,coef)
 importFrom(stats,pnorm)
 importFrom(stats,qchisq)
 importFrom(stats,qnorm)
 importFrom(stats,reformulate)
 importFrom(stats,setNames)
+importFrom(stats,vcov)
diff --git a/R/AddSuptBand.R b/R/AddSuptBand.R
index ccca421b..31f8c7e4 100644
--- a/R/AddSuptBand.R
+++ b/R/AddSuptBand.R
@@ -17,6 +17,8 @@
 #' for each event-study coefficient.
 #' @import estimatr
 #' @importFrom MASS mvrnorm
+#' @importFrom fixest coeftable
+#' @importFrom dplyr rename
 #' @keywords internal
 #' @noRd
 #'
@@ -47,16 +49,18 @@
 #'   eventstudy_coefficients = eventstudy_estimates$arguments$eventstudy_coefficients
 #')
 
-AddSuptBand <- function(estimates, num_sim = 1000, conf_level = .95, eventstudy_coefficients) {
+AddSuptBand <- function(model_estimates, num_sim = 1000, conf_level = .95, eventstudy_coefficients) {
 
-    if (! class(estimates) %in% c("lm_robust", "iv_robust")) {
+    if (! class(model_estimates) %in% c("lm_robust", "iv_robust", "fixest")) {
         stop("estimates is not a data frame with coefficient estimates and standard errors")
     }
     if (! is.numeric(num_sim) | num_sim %% 1 != 0 | num_sim <= 0) {stop("num_sim should be a natural number.")}
     if (! is.numeric(conf_level) | conf_level < 0 | conf_level > 1) {stop("conf_level should be a real number between 0 and 1, inclusive.")}
     if (! is.character(eventstudy_coefficients)) {stop("eventstudy_coefficients should be a character.")}
 
-    vcov_matrix_all <- estimates$vcov
+    fixest = class(model_estimates) == "fixest"
+
+    vcov_matrix_all <- if(fixest){vcov(model_estimates)} else {model_estimates$vcov}
     v_terms_to_keep <- colnames(vcov_matrix_all) %in% eventstudy_coefficients
     vcov_matrix <- vcov_matrix_all[v_terms_to_keep, v_terms_to_keep]
 
@@ -75,12 +79,18 @@ AddSuptBand <- function(estimates, num_sim = 1000, conf_level = .95, eventstudy_
         critical_value = t[floor(conf_level_num_sim) + 1]
     }
 
-    df_estimates_tidy <- estimatr::tidy(estimates)
+    df_estimates_tidy <- if(fixest){
+        coef_table <- model_estimates |> 
+        fixest::coeftable() |> 
+        as.data.frame()
+        coef_table$term <- rownames(coef_table)
+        coef_table |> 
+        dplyr::rename(estimate = Estimate, std.error = `Std. Error`) |> 
+        dplyr::select(term, estimate, std.error)
+    } else {estimatr::tidy(model_estimates)}
 
     df_estimates_tidy["suptband_lower"] <- df_estimates_tidy$estimate - (critical_value * df_estimates_tidy$std.error)
     df_estimates_tidy["suptband_upper"] <- df_estimates_tidy$estimate + (critical_value * df_estimates_tidy$std.error)
 
-
     return(df_estimates_tidy)
-
 }
diff --git a/R/EventStudy.R b/R/EventStudy.R
index 9d26f5b4..c8811b9b 100644
--- a/R/EventStudy.R
+++ b/R/EventStudy.R
@@ -34,6 +34,7 @@
 #' when there are anticipation effects. If set to FALSE, does not make the switch. Defaults to TRUE.
 #' @param allow_duplicate_id If TRUE, the function estimates a regression where duplicated ID-time rows are weighted by their duplication count. If FALSE, the function raises an error if duplicate unit-time keys exist in the input data. Default is FALSE.
 #' @param avoid_internal_copy If TRUE, the function avoids making an internal deep copy of the input data, and instead directly modifies the input data.table. Default is FALSE.
+#' @param kernel Accepts one of "estimatr" or "fixest". If "estimatr" is specified, uses the estimatr package for estimation. If "fixest" is specified, uses the fixest package for estimation. Defaults to "estimatr" (deprecated - will change to "fixest" in a future release).
 #'
 #' @return A list that contains, under "output", the estimation output as an lm_robust object, and under "arguments", the arguments passed to the function.
 #' @import dplyr
@@ -54,7 +55,8 @@
 #'     idvar = "id",
 #'     timevar = "t",
 #'     pre = 0, post = 3,
-#'     normalize = -1
+#'     normalize = -1,
+#'     kernel = "fixest"
 #'   )
 #'
 #' ### Access estimated model
@@ -63,7 +65,8 @@
 #' summary(eventstudy_model$output)
 #'
 #' ### data.frame of estimates
-#' estimatr::tidy(eventstudy_model$output)
+#' fixest::coeftable(eventstudy_model$output) # for kernel='fixest'
+#' # estimatr::tidy(eventstudy_model$output) # for kernel='estimatr'
 #'
 #' ### Access arguments
 #' eventstudy_model$arguments
@@ -83,6 +86,7 @@
 #'     pre  = 2, overidpre  = 4,
 #'     normalize = - 3,
 #'     cluster = TRUE,
+#'     kernel = "fixest",
 #'     anticipation_effects_normalization = TRUE
 #'   )
 #'
@@ -100,7 +104,8 @@
 #'     FE = TRUE, TFE = TRUE,
 #'     post = 0, overidpost = 0,
 #'     pre  = 0, overidpre  = 0,
-#'     cluster = TRUE
+#'     cluster = TRUE,
+#'     kernel = "fixest"
 #'   )
 #'
 #' summary(eventstudy_model_static$output)
@@ -117,7 +122,8 @@
 #'     idvar = "id",
 #'     timevar = "t",
 #'     pre = 0, post = 3,
-#'     normalize = -1
+#'     normalize = -1,
+#'     kernel = "fixest"
 #'   )
 #'
 #' summary(eventstudy_model_unbal$output)
@@ -136,7 +142,8 @@
 #'     post = 2, overidpost = 1,
 #'     pre  = 0, overidpre  = 3,
 #'     normalize = -1,
-#'     cluster = TRUE
+#'     cluster = TRUE,
+#'     kernel = "fixest"
 #'   )
 #'
 #' summary(eventstudy_model_iv$output)
@@ -145,10 +152,13 @@
 EventStudy <- function(estimator, data, outcomevar, policyvar, idvar, timevar, controls = NULL,
                        proxy = NULL, proxyIV = NULL, FE = TRUE, TFE = TRUE, post, overidpost = 1, pre, overidpre = post + pre,
                        normalize = -1 * (pre + 1), cluster = TRUE, anticipation_effects_normalization = TRUE,
-                       allow_duplicate_id = FALSE, avoid_internal_copy = FALSE) {
+                       allow_duplicate_id = FALSE, avoid_internal_copy = FALSE, kernel = "estimatr") {
 
     # Check for errors in arguments
     if (! estimator %in% c("OLS", "FHS")) {stop("estimator should be either 'OLS' or 'FHS'.")}
+    if (! kernel %in% c("estimatr", "fixest")) {stop("kernel should be either 'estimatr' or 'fixest'.")}
+    if (missing(kernel)) {warning("Argument 'kernel' was not specified; using 'estimatr' as default; we strongly recommend explicitly specifying a kernel because the default is scheduled to change.")}
+    if (kernel == "estimatr") {warning("'estimatr' selected as kernel. We no longer maintain it and will depreciate it in a future release. We recommend using 'fixest' instead.")}
     if (! is.data.frame(data))            {stop("data should be a data frame.")}
     for (var in c(idvar, timevar, outcomevar, policyvar)) {
         if ((! is.character(var))) {
@@ -336,14 +346,17 @@ EventStudy <- function(estimator, data, outcomevar, policyvar, idvar, timevar, c
     }
 
     if (estimator == "OLS") {
-        event_study_formula <- PrepareModelFormula(estimator, outcomevar, str_policy_vars,
-                                                   static, controls, proxy, proxyIV)
-
-        output       <- EventStudyOLS(event_study_formula, data, idvar, timevar, FE, TFE, cluster)
+        formula <- PrepareModelFormula(estimator, outcomevar, str_policy_vars,
+                                       static, controls, proxy, proxyIV,
+                                       kernel, idvar, timevar, FE, TFE)
+        
+        if (kernel == "estimatr") {
+            output <- EventStudyOLS(formula, data, idvar, timevar, FE, TFE, cluster)
+        } else if (kernel == "fixest") {
+            output <- EventStudyFEOLS(formula, data, idvar, timevar, FE, TFE, cluster)
+        }
         coefficients <- str_policy_vars
-    }
-    if (estimator == "FHS") {
-
+    } else if (estimator == "FHS") {
         if (is.null(proxyIV)) {
             Fstart <- 0
             str_fd_leads <- str_policy_vars[grepl("^z_fd_lead", str_policy_vars)]
@@ -360,10 +373,15 @@ EventStudy <- function(estimator, data, outcomevar, policyvar, idvar, timevar, c
                            ". To specify a different proxyIV use the proxyIV argument."))
         }
 
-        event_study_formula <- PrepareModelFormula(estimator, outcomevar, str_policy_vars,
-                                                   static, controls, proxy, proxyIV)
-
-        output       <- EventStudyFHS(event_study_formula, data, idvar, timevar, FE, TFE, cluster)
+        formula <- PrepareModelFormula(estimator, outcomevar, str_policy_vars,
+                                       static, controls, proxy, proxyIV,
+                                       kernel, idvar, timevar, FE, TFE)
+        
+        if (kernel == "estimatr") {
+            output <- EventStudyFHS(formula, data, idvar, timevar, FE, TFE, cluster)
+        } else if (kernel == "fixest") {
+            output <- EventStudyFEOLS_FHS(formula, data, idvar, timevar, FE, TFE, cluster)
+        }
         coefficients <- dplyr::setdiff(str_policy_vars, proxyIV)
     }
 
diff --git a/R/EventStudyFHS.R b/R/EventStudyFHS.R
index bccf9d85..36204e62 100644
--- a/R/EventStudyFHS.R
+++ b/R/EventStudyFHS.R
@@ -12,7 +12,8 @@
 #'
 #' @return A data.frame that contains the estimates for the event study coefficients.
 #' @import estimatr
-#' @importFrom stats qnorm pnorm
+#' @importFrom fixest feols
+#' @importFrom stats qnorm pnorm coef
 #' @keywords internal
 #' @noRd
 #'
@@ -127,3 +128,51 @@ EventStudyFHS <- function(prepared_model_formula, prepared_data,
   return(fhs_output)
 }
 
+
+EventStudyFEOLS_FHS <- function(formula, prepared_data,
+                                idvar, timevar, FE, TFE, cluster) {
+
+    if (! inherits(formula, "formula")) {stop("formula should be a formula")}
+    if (! is.data.frame(prepared_data)) {stop("data should be a data frame.")}
+    if (! is.character(idvar)) {stop("idvar should be a character.")}
+    if (! is.character(timevar)) {stop("timevar should be a character.")}
+    if (! is.logical(FE)) {stop("FE should be either TRUE or FALSE.")}
+    if (! is.logical(TFE)) {stop("TFE should be either TRUE or FALSE.")}
+    if (! is.logical(cluster)) {stop("cluster should be either TRUE or FALSE.")}
+    if (FE & !cluster) {stop("cluster=TRUE required when FE=TRUE.")}
+
+    if (cluster) {
+        vcov_arg <- as.formula(paste0("~", idvar))
+    } else {
+        vcov_arg <- "HC1"
+    }
+
+    fhs_output <- fixest::feols(
+        fml = formula,
+        data = prepared_data,
+        vcov = vcov_arg
+    )
+
+    if (FE & cluster) {
+        coefs <- coef(fhs_output)
+        N <- fhs_output$nobs
+        n <- length(unique(prepared_data[[idvar]]))
+        
+        if (TFE) {
+            K <- fhs_output$fixef_sizes[[timevar]] + length(coefs)
+        } else {
+            K <- 1 + length(coefs)
+        }
+
+        adjustment_factor <- (N - K) / (N - n - K + 1)
+        fhs_output$se <- fhs_output$se / sqrt(adjustment_factor)
+        fhs_output$cov.scaled <- fhs_output$cov.scaled / adjustment_factor
+
+        fhs_output$tstat <- coefs / fhs_output$se
+        fhs_output$pvalue <- 2 * stats::pnorm(abs(fhs_output$tstat), lower.tail = FALSE)
+        fhs_output$conf.low <- coefs - stats::qnorm(0.975) * fhs_output$se
+        fhs_output$conf.high <- coefs + stats::qnorm(0.975) * fhs_output$se
+    }
+
+    return(fhs_output)
+}
diff --git a/R/EventStudyOLS.R b/R/EventStudyOLS.R
index a0258bf5..c583f2d8 100644
--- a/R/EventStudyOLS.R
+++ b/R/EventStudyOLS.R
@@ -1,6 +1,6 @@
 #' Runs Ordinary Least Squares (OLS) with optional fixed effects and clustering
 #'
-#' @param prepared_model_formula A formula object created in [PrepareModelFormula()] that is passed to [EventStudy()].
+#' @param prepared_model_formula A formula object created in `PrepareModelFormula()` that is passed to [EventStudy()].
 #' @param prepared_data Data frame containing all of the parameters required for [EventStudy()] plus leads and
 #' lags of the first differenced policy variable and leads and lags of the policy variable.
 #' @param idvar Character indicating column of units.
@@ -12,6 +12,7 @@
 #'
 #' @return A data.frame that contains the estimates for the event study coefficients.
 #' @import estimatr
+#' @importFrom fixest feols ssc
 #' @keywords internal
 #' @noRd
 #'
@@ -123,3 +124,23 @@ EventStudyOLS <- function(prepared_model_formula, prepared_data,
 
     return(ols_output)
 }
+
+EventStudyFEOLS <- function(formula, prepared_data,
+                          idvar, timevar, FE, TFE, cluster) {
+
+    if (cluster) {
+        vcov_fixest <- as.formula(paste0("~", idvar))
+        small_sample_correction <- fixest::ssc(K.fixef = "full")
+    } else {
+        vcov_fixest <- "iid"
+        small_sample_correction <- fixest::ssc()
+    }
+
+    ols_output <- fixest::feols(
+        fml = formula,
+        data = prepared_data,
+        vcov = vcov_fixest,
+        ssc = small_sample_correction
+    )
+    return(ols_output)
+}
diff --git a/R/EventStudyPlot.R b/R/EventStudyPlot.R
index cb6f57e6..c2b5c9f2 100644
--- a/R/EventStudyPlot.R
+++ b/R/EventStudyPlot.R
@@ -29,7 +29,10 @@
 #' @return The Event-Study plot as a ggplot2 object.
 #' @import ggplot2 dplyr
 #' @import estimatr
+#' @importFrom stats vcov
 #' @importFrom rlang .data
+#' @importFrom fixest coeftable
+#' @importFrom dplyr rename select
 #' @importFrom data.table setorder
 #' @export
 #'
@@ -121,10 +124,19 @@ EventStudyPlot <- function(estimates,
 
 # Estimation Elements -----------------------------------------------------
 
-    df_estimates      <- estimates$output
-    df_estimates_tidy <- estimatr::tidy(estimates$output)
-
-    static_model <- nrow(df_estimates_tidy) == 1
+    model_estimates <- estimates$output
+    is_fixest <- class(model_estimates) == "fixest"
+    model_estimates_tidy <- if(is_fixest) {
+        coef_table <- model_estimates |> 
+        fixest::coeftable() |> 
+        as.data.frame()
+        coef_table$term <- rownames(coef_table)
+        coef_table |> 
+        dplyr::rename(estimate = Estimate, std.error = `Std. Error`) |> 
+        dplyr::select(term, estimate, std.error)
+    } else {estimatr::tidy(model_estimates)}
+
+    static_model <- length(coef(model_estimates)) == 1
     if (static_model) {
         stop("EventStudyPlot() does not support static models.")
     }
@@ -146,7 +158,7 @@ EventStudyPlot <- function(estimates,
     plot_supt <- if(!is.null(supt)) TRUE else FALSE
 
     if (plot_supt) {
-        df_estimates_tidy <- AddSuptBand(df_estimates, num_sim = 1000, conf_level = supt,
+        model_estimates_tidy <- AddSuptBand(model_estimates, num_sim = 1000, conf_level = supt,
                                          eventstudy_coefficients = eventstudy_coefficients)
     }
 
@@ -154,7 +166,7 @@ EventStudyPlot <- function(estimates,
 
     if (plot_CI) {
 
-        df_estimates_tidy <- AddCIs(df_estimates_tidy, eventstudy_coefficients, conf_level)
+        model_estimates_tidy <- AddCIs(model_estimates_tidy, eventstudy_coefficients, conf_level)
     }
 
 # Optionally Test For Pretrends/Levelling-Off -----------------------------
@@ -171,21 +183,18 @@ EventStudyPlot <- function(estimates,
 
         if (pre_event_coeffs & post_event_coeffs) {
             text_caption <- paste0(text_pretrends, " -- ", text_levelingoff)
-
         } else if (pre_event_coeffs & !post_event_coeffs) {
             text_caption <- text_pretrends
-
         } else if (!pre_event_coeffs & post_event_coeffs) {
             text_caption <- text_levelingoff
-
         }
     } else {
         text_caption <- NULL
     }
 
-
-    df_plt <- PreparePlottingData(df_estimates_tidy, policyvar,
-                                  post, overidpost, pre, overidpre, normalization_column, proxyIV)
+    df_plt <- PreparePlottingData(
+        model_estimates_tidy, policyvar,
+        post, overidpost, pre, overidpre, normalization_column, proxyIV)
 
 # Construct y breaks ------------------------------------------------------
 
@@ -269,10 +278,13 @@ EventStudyPlot <- function(estimates,
         coefficients <- df_plt$estimate
 
         # Add column and row in matrix of coefficients in index of norm columns
-        covar <- AddZerosCovar(estimates$output$vcov,
-                               eventstudy_coefficients,
-                               df_plt[df_plt$estimate==0, ]$term,
-                               df_plt$term)
+        vcov <- if(is_fixest) {vcov(estimates$output)} else {estimates$output$vcov}
+        covar <- AddZerosCovar(
+          vcov,
+          eventstudy_coefficients,
+          df_plt[df_plt$estimate == 0, ]$term,
+          df_plt$term
+        )
 
         inv_covar <- pracma::pinv(covar)
 
diff --git a/R/PrepareModelFormula.R b/R/PrepareModelFormula.R
index c4962042..65f69d29 100644
--- a/R/PrepareModelFormula.R
+++ b/R/PrepareModelFormula.R
@@ -1,5 +1,5 @@
 
-#' Prepares a formula object for use in [EventStudyOLS()] or [EventStudyFHS()]
+#' Prepares a formula object for use in `EventStudyOLS()` or `EventStudyFHS()`
 
 #' @param estimator Accepts one of "OLS" or "FHS". If "FHS" is specified, implements IV estimator in Freyaldenhoven et al. 2019.
 #' @param outcomevar Character indicating column of outcome variable.
@@ -10,6 +10,11 @@
 #' @param proxyIV Character of column to be used as an instrument. Should be specified if and only if estimator is specified as "FHS".
 #' If NULL, defaults to the strongest lead of the policy variable based on the first stage.
 #' @param static Indicates whether the model to be estimated is static. Defaults to FALSE.
+#' @param kernel Character indicating the estimation kernel. Accepts "estimatr" or "fixest". Defaults to "estimatr".
+#' @param idvar Character indicating the identifier variable for fixed effects. Required when kernel = "fixest".
+#' @param timevar Character indicating the time variable for fixed effects. Required when kernel = "fixest".
+#' @param FE Logical indicating whether to include unit fixed effects. Required when kernel = "fixest". Defaults to FALSE.
+#' @param TFE Logical indicating whether to include time fixed effects. Required when kernel = "fixest". Defaults to FALSE.
 #' @return A formula object to be passed to EventStudy
 #'
 #' @importFrom stats reformulate as.formula
@@ -39,7 +44,9 @@
 
 PrepareModelFormula <- function(estimator, outcomevar,
                                 str_policy_vars, static = FALSE,
-                                controls = NULL, proxy = NULL, proxyIV = NULL) {
+                                controls = NULL, proxy = NULL, proxyIV = NULL,
+                                kernel = "estimatr", idvar = NULL, timevar = NULL,
+                                FE = FALSE, TFE = FALSE) {
 
     if (! estimator %in% c("OLS", "FHS"))      {stop("estimator should be either 'OLS' or 'FHS'.")}
     if (! is.character(outcomevar))            {stop("outcomevar should be a character.")}
@@ -53,26 +60,105 @@ PrepareModelFormula <- function(estimator, outcomevar,
     if (! static & length(str_policy_vars) <= 1) {stop("str_policy_vars must have more than one variable with static = FALSE.")}
     if (  static & !is.null(proxyIV))            {stop("static model is not compatible with FHS estimator.")}
 
-    if (estimator == "OLS") {
-        reg_formula <- stats::reformulate(
-            termlabels = c(str_policy_vars, controls),
-            response = outcomevar,
-            intercept = FALSE
-        )
+    if (! kernel %in% c("estimatr", "fixest")) {stop("kernel should be either 'estimatr' or 'fixest'.")}
+    if (kernel == "fixest") {
+        if (is.null(idvar) | !is.character(idvar)) {stop("idvar must be specified as a character when kernel is 'fixest'.")}
+        if (is.null(timevar) | !is.character(timevar)) {stop("timevar must be specified as a character when kernel is 'fixest'.")}
+        if (! is.logical(FE)) {stop("FE should be a logical.")}
+        if (! is.logical(TFE)) {stop("TFE should be a logical.")}
     }
 
-    if (estimator == "FHS") {
-        exogenous <- c(str_policy_vars, controls)
-        exogenous <- exogenous[exogenous != proxy]
-        exogenous <- exogenous[exogenous != proxyIV]
-
-        reg_formula <- stats::as.formula(
-            paste(outcomevar, "~",
-            paste(c(exogenous, proxy), collapse="+"),
-            "|",
-            paste(c(exogenous, proxyIV), collapse="+"))
-        )
+    if (kernel == "estimatr") {
+        if (estimator == "OLS") {
+            reg_formula <- stats::reformulate(
+                termlabels = c(str_policy_vars, controls),
+                response = outcomevar,
+                intercept = FALSE
+            )
+        }
+
+        if (estimator == "FHS") {
+            exogenous <- c(str_policy_vars, controls)
+            exogenous <- exogenous[exogenous != proxy]
+            exogenous <- exogenous[exogenous != proxyIV]
+
+            reg_formula <- stats::as.formula(
+                paste(outcomevar, "~",
+                paste(c(exogenous, proxy), collapse="+"),
+                "|",
+                paste(c(exogenous, proxyIV), collapse="+"))
+            )
+        }
+    } else if (kernel == "fixest") {
+        regressors <- c(str_policy_vars, controls)
+
+        if (estimator == "OLS") {
+            if (FE | TFE) {
+                fes <- c()
+                if (FE) {
+                    fes <- c(fes, idvar)
+                }
+                if (TFE) {
+                    fes <- c(fes, timevar)
+                }
+
+                formula_str <- paste(
+                    outcomevar,
+                    "~",
+                    paste(regressors, collapse = " + "),
+                    "|",
+                    paste(fes, collapse = " + ")
+                )
+                reg_formula <- stats::as.formula(formula_str)
+            } else {
+                reg_formula <- stats::reformulate(
+                    termlabels = regressors,
+                    response = outcomevar,
+                    intercept = TRUE
+                )
+            }
+        } else if (estimator == "FHS") {
+            exogenous <- c(str_policy_vars, controls)
+            exogenous <- exogenous[exogenous != proxy]
+            exogenous <- exogenous[exogenous != proxyIV]
+
+            if (FE | TFE) {
+                fes <- c()
+                if (FE) {
+                    fes <- c(fes, idvar)
+                }
+                if (TFE) {
+                    fes <- c(fes, timevar)
+                }
+
+                formula_str <- paste(
+                    outcomevar,
+                    "~",
+                    paste(exogenous, collapse = " + "),
+                    "|",
+                    paste(fes, collapse = " + "),
+                    "|",
+                    proxy,
+                    "~",
+                    paste(c(exogenous, proxyIV), collapse = " + ")
+                )
+                reg_formula <- stats::as.formula(formula_str)
+            } else {
+                formula_str <- paste(
+                    outcomevar,
+                    "~",
+                    paste(exogenous, collapse = " + "),
+                    "|",
+                    proxy,
+                    "~",
+                    paste(c(exogenous, proxyIV), collapse = " + ")
+                )
+                reg_formula <- stats::as.formula(formula_str)
+            }
+        }
     }
 
     return(reg_formula)
 }
+
+
diff --git a/R/SmPathHelpers.R b/R/SmPathHelpers.R
index 70a08b33..b9674ccc 100644
--- a/R/SmPathHelpers.R
+++ b/R/SmPathHelpers.R
@@ -1,9 +1,9 @@
 # Add zero where normalized coefficient(s) should be in covar matrix
-AddZerosCovar <- function(vcov_matrix_all, eventstudy_coeffs, norm_column,
+AddZerosCovar <- function(vcov, eventstudy_coeffs, norm_column,
                           coeffs_order) {
 
-    v_terms_to_keep <- colnames(vcov_matrix_all) %in% eventstudy_coeffs
-    covar           <- vcov_matrix_all[v_terms_to_keep, v_terms_to_keep]
+    v_terms_to_keep <- colnames(vcov) %in% eventstudy_coeffs
+    covar           <- vcov[v_terms_to_keep, v_terms_to_keep]
 
     n_coefs      = length(coeffs_order)
     needed_zeros = length(norm_column)
diff --git a/R/TestLinear.R b/R/TestLinear.R
index a5b6e293..72c24102 100644
--- a/R/TestLinear.R
+++ b/R/TestLinear.R
@@ -14,6 +14,7 @@
 #'
 #' @return A data frame containing the F-statistic and p-value for the specified test(s).
 #' @importFrom car linearHypothesis
+#' @importFrom fixest fitstat
 #' @export
 #'
 #' @examples
@@ -36,7 +37,7 @@
 TestLinear <- function(estimates, test = NA, test_name = "User Test", pretrends = TRUE, leveling_off = TRUE){
     if (! is.list(estimates) | length(estimates) != 2){
         stop("estimates should be a list of length two, an output of EventStudy()")}
-    if ((! class(estimates$output) %in% c("lm_robust", "iv_robust")) | ! is.list(estimates$output)) {
+    if ((! class(estimates$output) %in% c("lm_robust", "iv_robust", "fixest")) | ! is.list(estimates$output)) {
         stop("The first element of estimates should be a list of class 'lm_robust' with coefficient estimates and standard errors")
     }
     if (! is.list(estimates$arguments) | ! is.list(estimates$arguments)) {
@@ -46,10 +47,14 @@ TestLinear <- function(estimates, test = NA, test_name = "User Test", pretrends
     if (! is.logical(pretrends)) {stop("pretrends should be a logical. Default value is TRUE")}
     if (! is.logical(leveling_off)) {stop("leveling_off should be a logical. Default value is TRUE")}
 
-    if(estimates$arguments$cluster == TRUE){
+    fixest = class(estimates$output) == "fixest"
 
-        estimates$output$df.residual <- estimates$output$nclusters - 1
+    if(estimates$arguments$cluster == TRUE){
 
+        estimates$output$df.residual <- ifelse(
+            fixest,
+            as.integer(fixest::fitstat(estimates$output, "g")),
+            estimates$output$nclusters) - 1
     }
 
     coefficients <- estimates$arguments$eventstudy_coefficients
diff --git a/README.md b/README.md
index 5d46caa9..1b4e48a1 100644
--- a/README.md
+++ b/README.md
@@ -50,7 +50,8 @@ estimates_ols <- EventStudy(
    idvar = "id",
    timevar = "t",
    controls = "x_r",
-   pre = 0,  post = 4
+   pre = 0,  post = 4,
+   kernel = "fixest"
 )
 
 plt <- EventStudyPlot(estimates = estimates_ols)
diff --git a/cran-comments.md b/cran-comments.md
index c9b60b88..ee2ac6bb 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -15,4 +15,4 @@ We checked 0 reverse dependencies, comparing R CMD check results across CRAN and
 
 ## Package changes
 
- * Implemented patch updates for compatibility with `ggplot2` 3.6.0
+ * Added support for `fixest` as the regression kernel. `fixest` will replace `estimatr` as the default kernel in a future release.
diff --git a/man/EventStudy.Rd b/man/EventStudy.Rd
index 6cb118ca..70a9da5e 100644
--- a/man/EventStudy.Rd
+++ b/man/EventStudy.Rd
@@ -24,7 +24,8 @@ EventStudy(
   cluster = TRUE,
   anticipation_effects_normalization = TRUE,
   allow_duplicate_id = FALSE,
-  avoid_internal_copy = FALSE
+  avoid_internal_copy = FALSE,
+  kernel = "estimatr"
 )
 }
 \arguments{
@@ -79,6 +80,8 @@ when there are anticipation effects. If set to FALSE, does not make the switch.
 \item{allow_duplicate_id}{If TRUE, the function estimates a regression where duplicated ID-time rows are weighted by their duplication count. If FALSE, the function raises an error if duplicate unit-time keys exist in the input data. Default is FALSE.}
 
 \item{avoid_internal_copy}{If TRUE, the function avoids making an internal deep copy of the input data, and instead directly modifies the input data.table. Default is FALSE.}
+
+\item{kernel}{Accepts one of "estimatr" or "fixest". If "estimatr" is specified, uses the estimatr package for estimation. If "fixest" is specified, uses the fixest package for estimation. Defaults to "estimatr" (deprecated - will change to "fixest" in a future release).}
 }
 \value{
 A list that contains, under "output", the estimation output as an lm_robust object, and under "arguments", the arguments passed to the function.
@@ -98,7 +101,8 @@ eventstudy_model <-
     idvar = "id",
     timevar = "t",
     pre = 0, post = 3,
-    normalize = -1
+    normalize = -1,
+    kernel = "fixest"
   )
 
 ### Access estimated model
@@ -107,7 +111,8 @@ eventstudy_model$output
 summary(eventstudy_model$output)
 
 ### data.frame of estimates
-estimatr::tidy(eventstudy_model$output)
+fixest::coeftable(eventstudy_model$output) # for kernel='fixest'
+# estimatr::tidy(eventstudy_model$output) # for kernel='estimatr'
 
 ### Access arguments
 eventstudy_model$arguments
@@ -127,6 +132,7 @@ eventstudy_model_dyn <-
     pre  = 2, overidpre  = 4,
     normalize = - 3,
     cluster = TRUE,
+    kernel = "fixest",
     anticipation_effects_normalization = TRUE
   )
 
@@ -144,7 +150,8 @@ eventstudy_model_static <-
     FE = TRUE, TFE = TRUE,
     post = 0, overidpost = 0,
     pre  = 0, overidpre  = 0,
-    cluster = TRUE
+    cluster = TRUE,
+    kernel = "fixest"
   )
 
 summary(eventstudy_model_static$output)
@@ -161,7 +168,8 @@ eventstudy_model_unbal <-
     idvar = "id",
     timevar = "t",
     pre = 0, post = 3,
-    normalize = -1
+    normalize = -1,
+    kernel = "fixest"
   )
 
 summary(eventstudy_model_unbal$output)
@@ -180,7 +188,8 @@ eventstudy_model_iv <-
     post = 2, overidpost = 1,
     pre  = 0, overidpre  = 3,
     normalize = -1,
-    cluster = TRUE
+    cluster = TRUE,
+    kernel = "fixest"
   )
 
 summary(eventstudy_model_iv$output)
diff --git a/tests/testthat/test-EventStudy.R b/tests/testthat/test-EventStudy.R
index 5c094dd3..57006857 100644
--- a/tests/testthat/test-EventStudy.R
+++ b/tests/testthat/test-EventStudy.R
@@ -172,24 +172,36 @@ test_that("tests that package and STATA output agree when post, overidpost, pre,
     overidpost <- 0
     normalize  <- -1
 
-    outputs <- EventStudy(estimator = "OLS", data = example_data, outcomevar = "y_base",
+    outputs_fixest <- EventStudy(estimator = "OLS", data = example_data, outcomevar = "y_base",
                           policyvar = "z", idvar = "id", timevar = "t",
                           FE = TRUE, TFE = TRUE,
                           post = post, pre = pre, overidpre = overidpre, overidpost = overidpost,
-                          normalize = normalize, cluster = TRUE, anticipation_effects_normalization = TRUE)
+                          normalize = normalize, cluster = TRUE, anticipation_effects_normalization = TRUE,
+                          kernel = "fixest")
+    outputs_estimatr <- EventStudy(estimator = "OLS", data = example_data, outcomevar = "y_base",
+                          policyvar = "z", idvar = "id", timevar = "t",
+                          FE = TRUE, TFE = TRUE,
+                          post = post, pre = pre, overidpre = overidpre, overidpost = overidpost,
+                          normalize = normalize, cluster = TRUE, anticipation_effects_normalization = TRUE,
+                          kernel = "estimatr")
 
-    coef_package <- outputs$output$coefficients[[1]]
-    std_package  <- outputs$output$std.error[[1]]
+    coef_fixest <- coef(outputs_fixest$output)[[1]]
+    std_fixest  <- fixest::se(outputs_fixest$output)[[1]]
+    coef_estimatr <- outputs_estimatr$output$coefficients[[1]]
+    std_estimatr  <- outputs_estimatr$output$std.error[[1]]
 
     STATA_output <- read.csv('./input/df_test_base_STATA_allzero.csv')
     coef_STATA <- STATA_output$coef[[1]]
     std_STATA  <- STATA_output$std_error[[1]]
 
     epsilon <- 10e-7
-    expect_equal(coef_package, coef_STATA, tolerance = epsilon)
-    expect_equal(std_package, std_STATA, tolerance = epsilon)
+    expect_equal(coef_fixest, coef_STATA, tolerance = epsilon)
+    expect_equal(std_fixest, std_STATA, tolerance = epsilon)
+    expect_equal(coef_estimatr, coef_STATA, tolerance = epsilon)
+    expect_equal(std_estimatr, std_STATA, tolerance = epsilon)
 })
 
+
 test_that("does not create shiftvalues of differenced variable when post + overidpost - 1 < 1", {
 
     post       <- 1
diff --git a/tests/testthat/test-EventStudyFHS.R b/tests/testthat/test-EventStudyFHS.R
index 3aac45ef..1bf2cc19 100644
--- a/tests/testthat/test-EventStudyFHS.R
+++ b/tests/testthat/test-EventStudyFHS.R
@@ -315,7 +315,7 @@ test_that("Coefficients and Standard Errors agree with base STATA", {
 
     bools <- c("TTT", "TFT", "FTF", "FTT", "FFF", "FFT")
 
-    for (i in length(bools)) {
+    for (i in 1:length(bools)) {
         bool <- bools[i]
         estimator       <-  "FHS"
         outcomevar      <- "y_base"
@@ -323,9 +323,6 @@ test_that("Coefficients and Standard Errors agree with base STATA", {
         controls        <- "x_r"
         proxy           <- "eta_m"
         proxyIV         <- "z_fd_lead3"
-        event_study_formula <- PrepareModelFormula(estimator, outcomevar, 
-                                               str_policy_vars, FALSE, controls, proxy, proxyIV)
-
         idvar   <- "id"
         timevar <- "t"
 
@@ -333,8 +330,18 @@ test_that("Coefficients and Standard Errors agree with base STATA", {
         TFE <- as.logical(substring(bool, 2, 2))
         cluster <- as.logical(substring(bool, 3, 3))
 
-        reg <- EventStudyFHS(
-            prepared_model_formula = event_study_formula,
+        formula_fixest <- PrepareModelFormula(
+            estimator, outcomevar, str_policy_vars, static = FALSE, controls, 
+            proxy, proxyIV, idvar, timevar, FE, TFE, kernel = "fixest")
+        formula_estimatr <- PrepareModelFormula(
+            estimator, outcomevar, str_policy_vars, static = FALSE, controls,
+            proxy, proxyIV
+        )
+        reg_fixest <- EventStudyFEOLS_FHS(
+            formula_fixest, df_EventStudyFHS_example, idvar, timevar, FE, TFE, cluster
+        )
+        reg_estimatr <- EventStudyFHS(
+            prepared_model_formula = formula_estimatr,
             prepared_data = df_EventStudyFHS_example,
             idvar = idvar,
             timevar = timevar,
@@ -346,23 +353,53 @@ test_that("Coefficients and Standard Errors agree with base STATA", {
         df_test_STATA <- read.csv("./input/df_test_base_STATA_FHS.csv")
         epsilon <- 10e-6
 
-        expect_equal(unname(reg$coefficients["z_fd"]),       df_test_STATA[df_test_STATA[1] ==         "z_fd",][[2*i]],      tolerance = epsilon)
-        expect_equal(unname(reg$coefficients["z_fd_lead2"]), df_test_STATA[df_test_STATA[1] ==   "z_fd_lead2",][[2*i]],      tolerance = epsilon)
-        expect_equal(unname(reg$coefficients["eta_m"]),      df_test_STATA[df_test_STATA[1] ==        "eta_m",][[2*i]],      tolerance = epsilon)
-        expect_equal(unname(reg$coefficients["z_fd_lag1"]),  df_test_STATA[df_test_STATA[1] ==    "z_fd_lag1",][[2*i]],      tolerance = epsilon)
-        expect_equal(unname(reg$coefficients["z_fd_lag2"]),  df_test_STATA[df_test_STATA[1] ==    "z_fd_lag2",][[2*i]],      tolerance = epsilon)
-        expect_equal(unname(reg$coefficients["z_lead3"]),    df_test_STATA[df_test_STATA[1] ==      "z_lead3",][[2*i]]*(-1), tolerance = epsilon)
-        expect_equal(unname(reg$coefficients["z_lag3"]),     df_test_STATA[df_test_STATA[1] ==       "z_lag3",][[2*i]],      tolerance = epsilon)
-        expect_equal(unname(reg$coefficients["x_r"]),        df_test_STATA[df_test_STATA[1] ==          "x_r",][[2*i]],      tolerance = epsilon)
-
-        expect_equal(unname(reg$std.error["z_fd"]),          df_test_STATA[df_test_STATA[1] ==       "z_fd",][[2*i+1]],      tolerance = epsilon)
-        expect_equal(unname(reg$std.error["z_fd_lead2"]),    df_test_STATA[df_test_STATA[1] == "z_fd_lead2",][[2*i+1]],      tolerance = epsilon)
-        expect_equal(unname(reg$std.error["eta_m"]),         df_test_STATA[df_test_STATA[1] ==      "eta_m",][[2*i+1]],      tolerance = epsilon)
-        expect_equal(unname(reg$std.error["z_fd_lag1"]),     df_test_STATA[df_test_STATA[1] ==  "z_fd_lag1",][[2*i+1]],      tolerance = epsilon)
-        expect_equal(unname(reg$std.error["z_fd_lag2"]),     df_test_STATA[df_test_STATA[1] ==  "z_fd_lag2",][[2*i+1]],      tolerance = epsilon)
-        expect_equal(unname(reg$std.error["z_lead3"]),       df_test_STATA[df_test_STATA[1] ==    "z_lead3",][[2*i+1]],      tolerance = epsilon)
-        expect_equal(unname(reg$std.error["z_lag3"]),        df_test_STATA[df_test_STATA[1] ==     "z_lag3",][[2*i+1]],      tolerance = epsilon)
-        expect_equal(unname(reg$std.error["x_r"]),           df_test_STATA[df_test_STATA[1] ==        "x_r",][[2*i+1]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$coefficients["z_fd"]),       df_test_STATA[df_test_STATA[1] ==         "z_fd",][[2*i]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$coefficients["z_fd_lead2"]), df_test_STATA[df_test_STATA[1] ==   "z_fd_lead2",][[2*i]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$coefficients["eta_m"]),      df_test_STATA[df_test_STATA[1] ==        "eta_m",][[2*i]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$coefficients["z_fd_lag1"]),  df_test_STATA[df_test_STATA[1] ==    "z_fd_lag1",][[2*i]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$coefficients["z_fd_lag2"]),  df_test_STATA[df_test_STATA[1] ==    "z_fd_lag2",][[2*i]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$coefficients["z_lead3"]),    df_test_STATA[df_test_STATA[1] ==      "z_lead3",][[2*i]]*(-1), tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$coefficients["z_lag3"]),     df_test_STATA[df_test_STATA[1] ==       "z_lag3",][[2*i]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$coefficients["x_r"]),        df_test_STATA[df_test_STATA[1] ==          "x_r",][[2*i]],      tolerance = epsilon)
+
+        expect_equal(unname(reg_estimatr$std.error["z_fd"]),          df_test_STATA[df_test_STATA[1] ==       "z_fd",][[2*i+1]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$std.error["z_fd_lead2"]),    df_test_STATA[df_test_STATA[1] == "z_fd_lead2",][[2*i+1]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$std.error["eta_m"]),         df_test_STATA[df_test_STATA[1] ==      "eta_m",][[2*i+1]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$std.error["z_fd_lag1"]),     df_test_STATA[df_test_STATA[1] ==  "z_fd_lag1",][[2*i+1]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$std.error["z_fd_lag2"]),     df_test_STATA[df_test_STATA[1] ==  "z_fd_lag2",][[2*i+1]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$std.error["z_lead3"]),       df_test_STATA[df_test_STATA[1] ==    "z_lead3",][[2*i+1]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$std.error["z_lag3"]),        df_test_STATA[df_test_STATA[1] ==     "z_lag3",][[2*i+1]],      tolerance = epsilon)
+        expect_equal(unname(reg_estimatr$std.error["x_r"]),           df_test_STATA[df_test_STATA[1] ==        "x_r",][[2*i+1]],      tolerance = epsilon)
+
+        epsilon <- 10e-4
+
+        coef_mappings <- list(
+            "z_fd" = "z_fd",
+            "z_fd_lead2" = "z_fd_lead2",
+            "fit_eta_m" = "eta_m",
+            "z_fd_lag1" = "z_fd_lag1",
+            "z_fd_lag2" = "z_fd_lag2",
+            "z_lead3" = "z_lead3",
+            "z_lag3" = "z_lag3",
+            "x_r" = "x_r"
+        )
+
+        # Test coefficients
+        for (r_name in names(coef_mappings)) {
+            stata_name <- coef_mappings[[r_name]]
+            expected <- df_test_STATA[df_test_STATA[1] == stata_name,][[2*i]]
+            if (r_name == "z_lead3") expected <- expected * (-1)  # STATA sign convention
+            tolerance <- abs(expected) * epsilon
+            expect_equal(unname(coef(reg_fixest)[r_name]), expected, tolerance = tolerance)
+        }
+
+        # Test standard errors
+        for (r_name in names(coef_mappings)) {
+            stata_name <- coef_mappings[[r_name]]
+            expected <- df_test_STATA[df_test_STATA[1] == stata_name,][[2*i+1]]
+            tolerance <- abs(expected) * epsilon
+            expect_equal(unname(fixest::se(reg_fixest)[r_name]), expected, tolerance = tolerance)
+        }
     }
 })
 
diff --git a/tests/testthat/test-EventStudyOLS.R b/tests/testthat/test-EventStudyOLS.R
index 265496fd..0aa893bc 100644
--- a/tests/testthat/test-EventStudyOLS.R
+++ b/tests/testthat/test-EventStudyOLS.R
@@ -333,7 +333,7 @@ test_that("Coefficients and Standard Errors agree with base STATA", {
     str_policy_vars <- c("z_lead3", "z_fd_lead3", "z_fd_lead2", "z_fd", "z_fd_lag1", "z_fd_lag2", "z_lag3")
     controls        <- "x_r"
 
-    event_study_formula <- PrepareModelFormula(estimator, outcomevar, str_policy_vars,
+    event_study_formula_estimatr <- PrepareModelFormula(estimator, outcomevar, str_policy_vars,
                                                static = F, controls = controls)
 
     idvar   <- "id"
@@ -343,8 +343,8 @@ test_that("Coefficients and Standard Errors agree with base STATA", {
     TFE     <- TRUE
     cluster <- TRUE
 
-    reg <- EventStudyOLS(
-            prepared_model_formula = event_study_formula,
+    reg_estimatr <- EventStudyOLS(
+            prepared_model_formula = event_study_formula_estimatr,
             prepared_data = df_test_EventStudyOLS,
             idvar = idvar,
             timevar = timevar,
@@ -357,21 +357,79 @@ test_that("Coefficients and Standard Errors agree with base STATA", {
 
     epsilon <- 10e-7
 
-    expect_equal(unname(reg$coefficients["z_fd"]), df_test_STATA[df_test_STATA["term"] == "zfd",][[2]], tolerance = epsilon)
-    expect_equal(unname(reg$coefficients["z_fd_lead2"]), df_test_STATA[df_test_STATA["term"] == "F2.zfd",][[2]], tolerance = epsilon)
-    expect_equal(unname(reg$coefficients["z_fd_lead3"]), df_test_STATA[df_test_STATA["term"] == "F3.zfd",][["coef"]][1], tolerance = epsilon)
-    expect_equal(unname(reg$coefficients["z_fd_lag1"]), df_test_STATA[df_test_STATA["term"] == "L.zfd",][[2]], tolerance = epsilon)
-    expect_equal(unname(reg$coefficients["z_fd_lag2"]), df_test_STATA[df_test_STATA["term"] == "L2.zfd",][[2]], tolerance = epsilon)
-    expect_equal(unname(reg$coefficients["z_lead3"]), -1 * df_test_STATA[df_test_STATA["term"] == "F3.z",][["coef"]], tolerance = epsilon)
-    expect_equal(unname(reg$coefficients["z_lag3"]), df_test_STATA[df_test_STATA["term"] == "L3.z",][[2]], tolerance = epsilon)
-    expect_equal(unname(reg$coefficients["x_r"]), df_test_STATA[df_test_STATA["term"] == "x_r",][[2]], tolerance = epsilon)
-
-    expect_equal(unname(reg$std.error["z_fd"]), df_test_STATA[df_test_STATA["term"] == "zfd",][[3]], tolerance = epsilon)
-    expect_equal(unname(reg$std.error["z_fd_lead2"]), df_test_STATA[df_test_STATA["term"] == "F2.zfd",][[3]], tolerance = epsilon)
-    expect_equal(unname(reg$std.error["z_fd_lead3"]), df_test_STATA[df_test_STATA["term"] == "F3.zfd",][["std_error"]][1], tolerance = epsilon)
-    expect_equal(unname(reg$std.error["z_fd_lag1"]), df_test_STATA[df_test_STATA["term"] == "L.zfd",][[3]], tolerance = epsilon)
-    expect_equal(unname(reg$std.error["z_fd_lag2"]),df_test_STATA[df_test_STATA["term"] == "L2.zfd",][[3]], tolerance = epsilon)
-    expect_equal(unname(reg$std.error["z_lead3"]), df_test_STATA[df_test_STATA["term"] == "F3.z",][["std_error"]], tolerance = epsilon)
-    expect_equal(unname(reg$std.error["z_lag3"]),df_test_STATA[df_test_STATA["term"] == "L3.z",][[3]], tolerance = epsilon)
-    expect_equal(unname(reg$std.error["x_r"]), df_test_STATA[df_test_STATA["term"] == "x_r",][[3]], tolerance = epsilon * 10)
+    expect_equal(unname(reg_estimatr$coefficients["z_fd"]), df_test_STATA[df_test_STATA["term"] == "zfd",][[2]], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$coefficients["z_fd_lead2"]), df_test_STATA[df_test_STATA["term"] == "F2.zfd",][[2]], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$coefficients["z_fd_lead3"]), df_test_STATA[df_test_STATA["term"] == "F3.zfd",][["coef"]][1], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$coefficients["z_fd_lag1"]), df_test_STATA[df_test_STATA["term"] == "L.zfd",][[2]], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$coefficients["z_fd_lag2"]), df_test_STATA[df_test_STATA["term"] == "L2.zfd",][[2]], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$coefficients["z_lead3"]), -1 * df_test_STATA[df_test_STATA["term"] == "F3.z",][["coef"]], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$coefficients["z_lag3"]), df_test_STATA[df_test_STATA["term"] == "L3.z",][[2]], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$coefficients["x_r"]), df_test_STATA[df_test_STATA["term"] == "x_r",][[2]], tolerance = epsilon)
+
+    expect_equal(unname(reg_estimatr$std.error["z_fd"]), df_test_STATA[df_test_STATA["term"] == "zfd",][[3]], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$std.error["z_fd_lead2"]), df_test_STATA[df_test_STATA["term"] == "F2.zfd",][[3]], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$std.error["z_fd_lead3"]), df_test_STATA[df_test_STATA["term"] == "F3.zfd",][["std_error"]][1], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$std.error["z_fd_lag1"]), df_test_STATA[df_test_STATA["term"] == "L.zfd",][[3]], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$std.error["z_fd_lag2"]),df_test_STATA[df_test_STATA["term"] == "L2.zfd",][[3]], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$std.error["z_lead3"]), df_test_STATA[df_test_STATA["term"] == "F3.z",][["std_error"]], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$std.error["z_lag3"]),df_test_STATA[df_test_STATA["term"] == "L3.z",][[3]], tolerance = epsilon)
+    expect_equal(unname(reg_estimatr$std.error["x_r"]), df_test_STATA[df_test_STATA["term"] == "x_r",][[3]], tolerance = epsilon * 10)
+
+    event_study_formula_fixest <- PrepareModelFormula(
+        estimator = "OLS",
+        outcomevar = outcomevar,
+        str_policy_vars = str_policy_vars,
+        static = FALSE,
+        controls = controls,
+        idvar = idvar,
+        timevar = timevar,
+        FE = FE,
+        TFE = TFE,
+        kernel = "fixest"
+    )
+
+    reg_fixest <- EventStudyFEOLS(
+        formula = event_study_formula_fixest,
+        prepared_data = df_test_EventStudyOLS,
+        idvar = idvar,
+        timevar = timevar,
+        FE = FE,
+        TFE = TFE,
+        cluster = cluster
+    )
+
+    coef_feols <- coef(reg_fixest)
+    se_feols <- fixest::se(reg_fixest)
+
+    epsilon <- 1e-4
+
+    coef_mappings <- list(
+        "z_fd" = "zfd",
+        "z_fd_lead2" = "F2.zfd",
+        "z_fd_lead3" = "F3.zfd",
+        "z_fd_lag1" = "L.zfd",
+        "z_fd_lag2" = "L2.zfd",
+        "z_lead3" = "F3.z",
+        "z_lag3" = "L3.z",
+        "x_r" = "x_r"
+    )
+
+    # Test coefficients
+    for (r_name in names(coef_mappings)) {
+        stata_term <- coef_mappings[[r_name]]
+        expected <- df_test_STATA[df_test_STATA["term"] == stata_term, "coef"]
+        if (r_name == "z_lead3") expected <- -1 * expected  # STATA sign convention
+        actual <- unname(coef_feols[r_name])
+        tolerance <- abs(expected) * epsilon
+        expect_equal(actual, expected, tolerance = tolerance)
+    }
+
+    # Test standard errors
+    for (r_name in names(coef_mappings)) {
+        stata_term <- coef_mappings[[r_name]]
+        expected <- df_test_STATA[df_test_STATA["term"] == stata_term, "std_error"]
+        actual <- unname(se_feols[r_name])
+        tolerance <- abs(expected) * epsilon
+        expect_equal(actual, expected, tolerance = tolerance)
+    }
 })
diff --git a/vignettes/documentation.Rmd b/vignettes/documentation.Rmd
index 8a06fe54..71002d40 100644
--- a/vignettes/documentation.Rmd
+++ b/vignettes/documentation.Rmd
@@ -55,8 +55,8 @@ results <- EventStudy(estimator = "OLS",
                       policyvar = "z",
                       idvar = "id",
                       timevar = "t",
-                      post = 3,
-                      pre = 0)
+                      post = 3, pre = 0,
+                      kernel = "fixest")
 ```
 ```{r Basic Eventstudy Example - Run Code, echo = FALSE}
 results <- EventStudy(estimator = "OLS",
@@ -65,8 +65,8 @@ results <- EventStudy(estimator = "OLS",
                       policyvar = "z",
                       idvar = "id",
                       timevar = "t",
-                      post = 3,
-                      pre = 0)
+                      post = 3, pre = 0,
+                      kernel = "fixest")
 ```
 ```{r Basic Eventstudy Example - Show Results 1, echo=TRUE, eval=TRUE}
 summary(results$output)
@@ -126,8 +126,8 @@ eventstudy_estimates_ols <- EventStudy(estimator = "OLS",
                                        policyvar = "z",
                                        idvar = "id",
                                        timevar = "t",
-                                       post = 3,
-                                       pre = 0)
+                                       post = 3, pre = 0,
+                                       kernel = "fixest")
  
 EventStudyPlot(estimates = eventstudy_estimates_ols,
                xtitle = "Event time",