diff --git a/DESCRIPTION b/DESCRIPTION
index c3b5994..81c7286 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -39,6 +39,7 @@ Depends:
     R (>= 3.4.0),
     fabletools (>= 0.3.0)
 Imports:
+    cli,
     Rcpp (>= 0.11.0),
     rlang (>= 0.4.6),
     stats,
diff --git a/R/00_specials.R b/R/00_specials.R
index a9376ae..8d0f98c 100644
--- a/R/00_specials.R
+++ b/R/00_specials.R
@@ -13,7 +13,7 @@ model_xreg <- function(...) {
 }
 
 no_xreg <- function(...) {
-  abort("Exogenous regressors are not supported for this model type.")
+  cli::cli_abort("Exogenous regressors are not supported for this model type.")
 }
 
 trend <- function(x, knots = NULL, origin = NULL) {
@@ -91,10 +91,10 @@ fourier.tbl_ts <- function(x, period, K, origin = NULL) {
 
 fourier.numeric <- function(x, period, K, origin = NULL) {
   if (length(period) != length(K)) {
-    abort("Number of periods does not match number of orders")
+    cli::cli_abort("Number of periods does not match number of orders")
   }
   if (any(2 * K > period)) {
-    abort("K must be not be greater than period/2")
+    cli::cli_abort("{.arg K} must be not be greater than period/2")
   }
 
   fourier_exprs <- map2(
diff --git a/R/VARIMA.R b/R/VARIMA.R
index b8bb3dd..ed8e36e 100644
--- a/R/VARIMA.R
+++ b/R/VARIMA.R
@@ -3,7 +3,7 @@ train_varima <- function(.data, specials, identification = NULL, ...) {
   y <- invoke(cbind, lapply(unclass(.data)[measured_vars(.data)], as.double))
   
   if(any(colnames(specials$xreg[[1]]) != "(Intercept)")) {
-    stop("Exogenous regressors for VARIMA are not yet supported.")
+    cli::cli_abort("Exogenous regressors for VARIMA are not yet supported.")
   }
   
   p <- specials$pdq[[1]]$p
@@ -20,7 +20,7 @@ train_varima <- function(.data, specials, identification = NULL, ...) {
     }
   }
   
-  require_package("MTS")
+  check_installed("MTS")
   utils::capture.output(
     fit <- if (identification == "kronecker_indices") {
       MTS::Kronfit(
@@ -39,7 +39,7 @@ train_varima <- function(.data, specials, identification = NULL, ...) {
       )
     } else {
       if(length(p) != 1 || length(q) != 1) {
-        stop("Model selection is not yet supported, please specify `p` and `q` exactly.")
+        cli::cli_abort("Model selection is not yet supported, please specify {.arg p} and {.arg q} exactly.")
       }
       MTS::VARMA(
         yd,
@@ -72,7 +72,7 @@ specials_varima <- new_specials(
     as.list(environment())
   },
   PDQ = function(P, D, Q, period = NULL) {
-    stop("Seasonal VARIMA models are not yet supported.")
+    cli::cli_abort("Seasonal VARIMA models are not yet supported.")
   },
   common_xregs,
   xreg = special_xreg(default_intercept = TRUE),
@@ -580,4 +580,4 @@ IRF.VARIMA <- function(x, new_data, specials, impulse = NULL, orthogonal = FALSE
   irf[colnames(x$data)] <- split(irf$.sim, col(irf$.sim))
   irf$.innov <- irf$.sim <- NULL
   irf
-}
\ No newline at end of file
+}
diff --git a/R/ar.R b/R/ar.R
index ad160fc..b62acfb 100644
--- a/R/ar.R
+++ b/R/ar.R
@@ -72,7 +72,7 @@ AR <- function(formula, ic = c("aicc", "aic", "bic"), ...) {
 specials_ar <- new_specials(
   order = function(p = 0:15, fixed = list()) {
     if (any(p < 0)) {
-      warn("The AR order must be non-negative. Only non-negative orders will be considered.")
+      cli::cli_warn("The AR order must be non-negative. Only non-negative orders will be considered.")
       p <- p[p >= 0]
     }
     list(p = p, fixed = fixed)
@@ -158,9 +158,9 @@ estimate_ar <- function(x, p, xreg, constant, fixed) {
   XX <- t(X_est) %*% X_est
   rank <- qr(XX)$rank
   if (rank != nrow(XX)) {
-    warning(paste("model order: ", p, "singularities in the computation of the projection matrix", 
-                  "results are only valid up to model order", 
-                  p - 1L), domain = NA)
+    cli::cli_warn(
+      "model order: {p} singularities in the computation of the projection matrix results are only valid up to model order {p - 1L}"
+    )
     return(NULL)
   }
   P <- if (ncol(XX) > 0) 
diff --git a/R/arima.R b/R/arima.R
index 398b93c..9694e44 100644
--- a/R/arima.R
+++ b/R/arima.R
@@ -8,12 +8,12 @@ train_arima <- function(.data, specials,
                         unitroot_spec = unitroot_options(), trace = FALSE,
                         fixed = NULL, method = NULL, ...) {
   if (length(measured_vars(.data)) > 1) {
-    abort("Only univariate responses are supported by ARIMA.")
+    cli::cli_abort("Only univariate responses are supported by ARIMA.")
   }
   
   # Get args
   if(length(specials$pdq) > 1 || length(specials$PDQ) > 1){
-    warn("Only one special for `pdq()` and `PDQ()` is allowed, defaulting to the first usage")
+    cli::cli_warn("Only one special for {.fn pdq} and {.fn PDQ} is allowed, defaulting to the first usage")
   }
   pdq <- specials$pdq[[1]]
   PDQ <- specials$PDQ[[1]]
@@ -24,7 +24,7 @@ train_arima <- function(.data, specials,
   y <- x <- ts(unclass(.data)[[measured_vars(.data)]], frequency = period)
 
   if (all(is.na(y))) {
-    abort("All observations are missing, a model cannot be estimated without data.")
+    cli::cli_abort("All observations are missing, a model cannot be estimated without data.")
   }
 
   # Get xreg
@@ -47,7 +47,7 @@ train_arima <- function(.data, specials,
       
       # Remove deficient regressors
       if(!is_empty(bad_regressors)){
-        warn(sprintf(
+        cli::cli_warn(sprintf(
           "Provided exogenous regressors are rank deficient, removing regressors: %s",
           paste("`", colnames(xreg)[bad_regressors], "`", sep = "", collapse = ", ")
         ))
@@ -72,8 +72,8 @@ train_arima <- function(.data, specials,
   if (NROW(model_opts) > 1) {
     model_opts <- filter(model_opts, !!enexpr(order_constraint))
     if (NROW(model_opts) == 0) {
-      if (mostly_specified) warn(mostly_specified_msg)
-      abort("There are no ARIMA models to choose from after imposing the `order_constraint`, please consider allowing more models.")
+      if (mostly_specified) cli::cli_warn(mostly_specified_msg)
+      cli::cli_abort("There are no ARIMA models to choose from after imposing the {.arg order_constraint}, please consider allowing more models.")
     }
     wrap_arima <- possibly(quietly(stats::arima), NULL)
   }
@@ -90,7 +90,7 @@ train_arima <- function(.data, specials,
 
   # Choose seasonal differencing
   if (length(seas_D <- unique(model_opts$D)) > 1) {
-    require_package("feasts")
+    check_installed("feasts")
     # Valid xregs
 
     if (!is.null(xreg)) {
@@ -111,7 +111,7 @@ train_arima <- function(.data, specials,
   x <- diff(x, lag = period, differences = seas_D)
   diff_xreg <- diff(xreg, lag = period, differences = seas_D)
   if (length(seas_d <- unique(model_opts$d)) > 1) {
-    require_package("feasts")
+    check_installed("feasts")
 
     # Valid xregs
     if (!is.null(xreg)) {
@@ -132,12 +132,12 @@ train_arima <- function(.data, specials,
   }
 
   # Check number of differences selected
-  if (length(seas_D) != 1) abort("Could not find appropriate number of seasonal differences.")
-  if (length(seas_d) != 1) abort("Could not find appropriate number of non-seasonal differences.")
+  if (length(seas_D) != 1) cli::cli_abort("Could not find appropriate number of seasonal differences.")
+  if (length(seas_d) != 1) cli::cli_abort("Could not find appropriate number of non-seasonal differences.")
   if (seas_D >= 2) {
-    warn("Having more than one seasonal difference is not recommended. Please consider using only one seasonal difference.")
+    cli::cli_warn("Having more than one seasonal difference is not recommended. Please consider using only one seasonal difference.")
   } else if (seas_D + seas_d > 2) {
-    warn("Having 3 or more differencing operations is not recommended. Please consider reducing the total number of differences.")
+    cli::cli_warn("Having 3 or more differencing operations is not recommended. Please consider reducing the total number of differences.")
   }
 
   # Find best model
@@ -236,7 +236,7 @@ train_arima <- function(.data, specials,
       method <- "CSS-ML"
     }
   } else {
-    if(isTRUE(approximation)) warn("Estimating ARIMA models with approximation is not supported when `method` is specified.")
+    if(isTRUE(approximation)) cli::cli_warn("Estimating ARIMA models with approximation is not supported when {.arg method} is specified.")
     approximation <- FALSE
   }
   
@@ -252,7 +252,7 @@ train_arima <- function(.data, specials,
   }
   
   if (any((model_opts$d + model_opts$D > 1) & model_opts$constant)) {
-    warn("Model specification induces a quadratic or higher order polynomial trend. 
+    cli::cli_warn("Model specification induces a quadratic or higher order polynomial trend. 
 This is generally discouraged, consider removing the constant or reducing the number of differences.")
   }
   constant <- unique(model_opts$constant)
@@ -337,8 +337,12 @@ This is generally discouraged, consider removing the constant or reducing the nu
   }
 
   if (is.null(best)) {
-    if (mostly_specified) warn(mostly_specified_msg)
-    abort("Could not find an appropriate ARIMA model.\nThis is likely because automatic selection does not select models with characteristic roots that may be numerically unstable.\nFor more details, refer to https://otexts.com/fpp3/arima-r.html#plotting-the-characteristic-roots")
+    if (mostly_specified) cli::cli_warn(mostly_specified_msg)
+    cli::cli_abort(c(
+      "Could not find an appropriate ARIMA model.",
+      "This is likely because automatic selection does not select models with characteristic roots that may be numerically unstable.",
+      "For more details, refer to {.url https://otexts.com/fpp3/arima-r.html#plotting-the-characteristic-roots}"
+    ))
   }
 
   # Compute ARMA roots
@@ -420,7 +424,7 @@ specials_arima <- new_specials(
     p_init <- p[which.min(abs(p - p_init))]
     q_init <- q[which.min(abs(q - q_init))]
     if(!all(grepl("^(ma|ar)\\d+", names(fixed)))){
-      abort("The 'fixed' coefficients for pdq() must begin with ar or ma, followed by a lag number.")
+      cli::cli_abort("The {.arg fixed} coefficients for {.fn pdq} must begin with ar or ma, followed by a lag number.")
     }
     as.list(environment())
   },
@@ -429,7 +433,7 @@ specials_arima <- new_specials(
                  fixed = list()) {
     period <- get_frequencies(period, self$data, .auto = "smallest")
     if (period < 1) {
-      abort("The seasonal period must be greater than or equal to 1.")
+      cli::cli_abort("The seasonal period must be greater than or equal to 1.")
     } else if (period == 1) {
       # Not seasonal
       P <- 0
@@ -440,13 +444,13 @@ specials_arima <- new_specials(
       P <- P[P <= floor(NROW(self$data) / 3 / period)]
       Q <- Q[Q <= floor(NROW(self$data) / 3 / period)]
       if(length(P) == 0 || length(Q) == 0) {
-        abort("Not enough data to estimate a model with those options of P and Q. Consider allowing smaller values of P and Q to be selected.")
+        cli::cli_abort("Not enough data to estimate a model with those options of P and Q. Consider allowing smaller values of P and Q to be selected.")
       }
     }
     P_init <- P[which.min(abs(P - P_init))]
     Q_init <- Q[which.min(abs(Q - Q_init))]
     if(!all(grepl("^(sma|sar)\\d+", names(fixed)))){
-      abort("The 'fixed' coefficients for PDQ() must begin with sar or sma, followed by a lag number.")
+      cli::cli_abort("The {.arg fixed} coefficients for {.fn PDQ} must begin with sar or sma, followed by a lag number.")
     }
     as.list(environment())
   },
@@ -716,7 +720,7 @@ residuals.ARIMA <- function(object, type = c("innovation", "regression"), ...) {
     object$est[[".regression_resid"]]
   }
   else {
-    abort(sprintf('Residuals of `type = "%s"` are not supported by ARIMA models', type))
+    cli::cli_abort("Residuals of {.code type = {.val type}} are not supported by ARIMA models")
   }
 }
 
diff --git a/R/checks.R b/R/checks.R
index f9b19fa..f5d7463 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -1,19 +1,19 @@
 check_gaps <- function(x) {
   if (any(tsibble::has_gaps(x)[[".gaps"]])) {
-    abort(sprintf("%s contains implicit gaps in time. You should check your data and convert implicit gaps into explicit missing values using `tsibble::fill_gaps()` if required.", deparse(substitute(x))))
+    cli::cli_abort(sprintf("%s contains implicit gaps in time. You should check your data and convert implicit gaps into explicit missing values using {.fn tsibble::fill_gaps} if required.", deparse(substitute(x))))
   }
 }
 
 check_regular <- function(x) {
   if (!is_regular(x)) {
-    abort(sprintf("%s is an irregular time series, which this model does not support. You should consider if your data can be made regular, and use `tsibble::update_tsibble(%s, regular = TRUE)` if appropriate.", deparse(substitute(x)), deparse(substitute(x))))
+    cli::cli_abort(sprintf("%s is an irregular time series, which this model does not support. You should consider if your data can be made regular, and use {.fn tsibble::update_tsibble(%s, regular = TRUE)} if appropriate.", deparse(substitute(x)), deparse(substitute(x))))
   }
 }
 
 check_ordered <- function(x) {
   if (!is_ordered(x)) {
-    abort(sprintf(
-      "%s is an unordered time series. To use this model, you first must sort the data in time order using `dplyr::arrange(%s, %s)`",
+    cli::cli_abort(sprintf(
+      "%s is an unordered time series. To use this model, you first must sort the data in time order using {.fn dplyr::arrange(%s, %s)}",
       deparse(substitute(x)), paste(c(deparse(substitute(x)), key_vars(x)), collapse = ", "), index_var(x)
     ))
   }
@@ -24,6 +24,6 @@ all_tsbl_checks <- function(.data) {
   check_regular(.data)
   check_ordered(.data)
   if (NROW(.data) == 0) {
-    abort("There is no data to model. Please provide a dataset with at least one observation.")
+    cli::cli_abort("There is no data to model. Please provide a dataset with at least one observation.")
   }
 }
diff --git a/R/croston.R b/R/croston.R
index 6a550f9..9e83f3d 100644
--- a/R/croston.R
+++ b/R/croston.R
@@ -99,10 +99,10 @@ CROSTON <- function(
 specials_croston <- new_specials(
   demand = function(initial = NULL, param = NULL, param_range = c(0, 1)) {
     if (!is.null(initial) && initial < 0) {
-      abort("The initial demand for Croston's method must be non-negative")
+      cli::cli_abort("The initial demand for Croston's method must be non-negative")
     }
     if (param_range[1] > param_range[2]) {
-      rlang::abort("Lower param limits must be less than upper limits")
+      cli::cli_abort("Lower param limits must be less than upper limits")
     }
 
     as.list(environment())
@@ -111,11 +111,11 @@ specials_croston <- new_specials(
     method <- match.arg(method)
 
     if (!is.null(initial) && initial < 1) {
-      abort("The initial interval for Croston's method must be greater than (or equal to) 1.")
+      cli::cli_abort("The initial interval for Croston's method must be greater than (or equal to) 1.")
     }
 
     if (param_range[1] > param_range[2]) {
-      rlang::abort("Lower param limits must be less than upper limits")
+      cli::cli_abort("Lower param limits must be less than upper limits")
     }
 
     as.list(environment())
@@ -125,7 +125,7 @@ specials_croston <- new_specials(
 
 train_croston <- function(.data, specials, opt_crit = "mse", type = "croston", ...) {
   if (length(measured_vars(.data)) > 1) {
-    abort("Only univariate responses are supported by Croston's method.")
+    cli::cli_abort("Only univariate responses are supported by Croston's method.")
   }
 
   # Get response
@@ -133,13 +133,13 @@ train_croston <- function(.data, specials, opt_crit = "mse", type = "croston", .
 
   # Check data
   if (any(y < 0)) {
-    abort("All observations must be non-negative for Croston's method.")
+    cli::cli_abort("All observations must be non-negative for Croston's method.")
   }
 
   non_zero <- which(y != 0)
 
   if (length(non_zero) < 2) {
-    abort("At least two non-zero values are required to use Croston's method.")
+    cli::cli_abort("At least two non-zero values are required to use Croston's method.")
   }
 
   # Get specials
diff --git a/R/ets.R b/R/ets.R
index a3a5dd8..a85c83f 100644
--- a/R/ets.R
+++ b/R/ets.R
@@ -1,14 +1,14 @@
 train_ets <- function(.data, specials, opt_crit,
                       nmse, bounds, ic, restrict = TRUE, ...) {
   if (length(measured_vars(.data)) > 1) {
-    abort("Only univariate responses are supported by ETS.")
+    cli::cli_abort("Only univariate responses are supported by ETS.")
   }
 
   # Rebuild `ets` arguments
   ets_spec <- specials[c("error", "trend", "season")]
   map(ets_spec, function(.x) {
     if (length(.x) > 1) {
-      abort("Only one special of each type is allowed for ETS.")
+      cli::cli_abort("Only one special of each type is allowed for ETS.")
     }
   })
   ets_spec <- unlist(ets_spec, recursive = FALSE)
@@ -18,7 +18,7 @@ train_ets <- function(.data, specials, opt_crit,
   idx <- unclass(.data)[[index_var(.data)]]
 
   if (any(is.na(y))) {
-    abort("ETS does not support missing values.")
+    cli::cli_abort("ETS does not support missing values.")
   }
 
   # Build possible models
@@ -47,7 +47,7 @@ train_ets <- function(.data, specials, opt_crit,
   }
 
   if (NROW(model_opts) == 0) {
-    abort("No valid ETS models have been allowed. Consider allowing different (more stable) models, or enabling the restricted models with `restrict = FALSE`.")
+    cli::cli_abort("No valid ETS models have been allowed. Consider allowing different (more stable) models, or enabling the restricted models with {.code restrict = FALSE}.")
   }
 
   # Find best model
@@ -78,7 +78,7 @@ train_ets <- function(.data, specials, opt_crit,
   ic <- pmap_dbl(model_opts, compare_ets)
 
   if (is.null(best)) {
-    abort(last_error$message)
+    cli::cli_abort(last_error$message)
   }
 
   best_spec <- model_opts[which.min(ic), ]
@@ -114,7 +114,7 @@ train_ets <- function(.data, specials, opt_crit,
 specials_ets <- new_specials(
   error = function(method = c("A", "M")) {
     if (!all(is.element(method, c("A", "M")))) {
-      stop("Invalid error type")
+      cli::cli_abort("Invalid error type")
     }
     list(method = method)
   },
@@ -123,16 +123,16 @@ specials_ets <- new_specials(
                    beta = NULL, beta_range = c(1e-04, 0.9999),
                    phi = NULL, phi_range = c(0.8, 0.98)) {
     if (!all(is.element(method, c("N", "A", "Ad", "M", "Md")))) {
-      stop("Invalid trend type")
+      cli::cli_abort("Invalid trend type")
     }
     if (alpha_range[1] > alpha_range[2]) {
-      abort("Lower alpha limits must be less than upper limits")
+      cli::cli_abort("Lower alpha limits must be less than upper limits")
     }
     if (beta_range[1] > beta_range[2]) {
-      abort("Lower beta limits must be less than upper limits")
+      cli::cli_abort("Lower beta limits must be less than upper limits")
     }
     if (phi_range[1] > phi_range[2]) {
-      abort("Lower phi limits must be less than upper limits")
+      cli::cli_abort("Lower phi limits must be less than upper limits")
     }
     if(!is.null(alpha)) alpha_range <- rep(alpha, 2)
     if(!is.null(beta)) beta_range <- rep(beta, 2)
@@ -147,10 +147,10 @@ specials_ets <- new_specials(
   season = function(method = c("N", "A", "M"), period = NULL,
                     gamma = NULL, gamma_range = c(1e-04, 0.9999)) {
     if (!all(is.element(method, c("N", "A", "M")))) {
-      abort("Invalid season type")
+      cli::cli_abort("Invalid season type")
     }
     if (gamma_range[1] > gamma_range[2]) {
-      abort("Lower gamma limits must be less than upper limits")
+      cli::cli_abort("Lower gamma limits must be less than upper limits")
     }
     if(!is.null(gamma)) gamma_range <- rep(gamma, 2)
 
@@ -160,14 +160,14 @@ specials_ets <- new_specials(
     }
     if (m > 24) {
       if (!is.element("N", method)) {
-        abort("Seasonal periods (`period`) of length greather than 24 are not supported by ETS.")
+        cli::cli_abort("Seasonal periods ({.arg period}) of length greather than 24 are not supported by ETS.")
       } else if (length(method) > 1) {
-        warn("Seasonal periods (`period`) of length greather than 24 are not supported by ETS. Seasonality will be ignored.")
+        cli::cli_warn("Seasonal periods ({.arg period}) of length greather than 24 are not supported by ETS. Seasonality will be ignored.")
         method <- "N"
       }
     }
     if (is_empty(method)) {
-      abort("A seasonal ETS model cannot be used for this data.")
+      cli::cli_abort("A seasonal ETS model cannot be used for this data.")
     }
     list(method = method, gamma = gamma, gamma_range = gamma_range, period = m)
   },
@@ -362,14 +362,14 @@ forecast.ETS <- function(object, new_data, specials = NULL, simulate = FALSE, bo
 #' @export
 generate.ETS <- function(x, new_data, specials, bootstrap = FALSE, ...) {
   if (!is_regular(new_data)) {
-    abort("Simulation new_data must be regularly spaced")
+    cli::cli_abort("Simulation new_data must be regularly spaced")
   }
 
   start_idx <- min(new_data[[index_var(new_data)]])
   start_pos <- match(start_idx - default_time_units(interval(new_data)), x$states[[index_var(x$states)]])
 
   if (is.na(start_pos)) {
-    abort("The first observation index of simulation data must be within the model's training set.")
+    cli::cli_abort("The first observation index of simulation data must be within the model's training set.")
   }
 
   initstate <- as.numeric(x$states[start_pos, measured_vars(x$states)])
@@ -414,7 +414,7 @@ generate.ETS <- function(x, new_data, specials, bootstrap = FALSE, ...) {
     )[[11]])
   
   if (is.na(result[[".sim"]][1])) {
-    stop("Problem with multiplicative damped trend")
+    cli::cli_abort("Problem with multiplicative damped trend")
   }
 
   result
@@ -720,4 +720,4 @@ initial_ets_states <- function(object) {
   )
   colnames(states_init) <- unsplit(states_names, states_type)
   states_init
-}
\ No newline at end of file
+}
diff --git a/R/etsmodel.R b/R/etsmodel.R
index fea5b25..cc08432 100644
--- a/R/etsmodel.R
+++ b/R/etsmodel.R
@@ -46,7 +46,7 @@ etsmodel <- function(y, m, errortype, trendtype, seasontype, damped,
   }
 
   if (!check.param(alpha, beta, gamma, phi, lower, upper, bounds, m)) {
-    abort(sprintf(
+    cli::cli_abort(sprintf(
       "Parameters out of range for ETS(%s,%s%s,%s) model",
       errortype, trendtype, ifelse(damped, "d", ""), seasontype
     ))
@@ -61,7 +61,7 @@ etsmodel <- function(y, m, errortype, trendtype, seasontype, damped,
 
   np <- length(par)
   if (np >= length(y) - 1) { # Not enough data to continue
-    abort("Not enough data to estimate this ETS model.")
+    cli::cli_abort("Not enough data to estimate this ETS model.")
   }
 
   env <- etsTargetFunctionInit(
@@ -161,12 +161,12 @@ etsTargetFunctionInit <- function(par, y, nstate, errortype, trendtype, seasonty
   names(par.noopt) <- pnames2
   alpha <- c(par["alpha"], par.noopt["alpha"])["alpha"]
   if (is.na(alpha)) {
-    stop("alpha problem!")
+    cli::cli_abort("alpha problem!")
   }
   if (trendtype != "N") {
     beta <- c(par["beta"], par.noopt["beta"])["beta"]
     if (is.na(beta)) {
-      stop("beta Problem!")
+      cli::cli_abort("beta Problem!")
     }
   }
   else {
@@ -175,7 +175,7 @@ etsTargetFunctionInit <- function(par, y, nstate, errortype, trendtype, seasonty
   if (seasontype != "N") {
     gamma <- c(par["gamma"], par.noopt["gamma"])["gamma"]
     if (is.na(gamma)) {
-      stop("gamma Problem!")
+      cli::cli_abort("gamma Problem!")
     }
   }
   else {
@@ -185,7 +185,7 @@ etsTargetFunctionInit <- function(par, y, nstate, errortype, trendtype, seasonty
   if (damped) {
     phi <- c(par["phi"], par.noopt["phi"])["phi"]
     if (is.na(phi)) {
-      stop("phi Problem!")
+      cli::cli_abort("phi Problem!")
     }
   }
   else {
@@ -259,7 +259,7 @@ initparam <- function(alpha, beta, gamma, phi, trendtype, seasontype, damped, lo
     lower[1L:3L] <- 0
     upper[1L:3L] <- 1e-3
   } else if (any(lower > upper)) {
-    stop("Inconsistent parameter boundaries")
+    cli::cli_abort("Inconsistent parameter boundaries")
   }
 
   # Select alpha
@@ -344,7 +344,7 @@ initstate <- function(y, m, trendtype, seasontype) {
     # Do decomposition
     n <- length(y)
     if (n < 4) {
-      stop("You've got to be joking (not enough data).")
+      cli::cli_abort("You've got to be joking (not enough data).")
     } else if (n < 3 * m) # Fit simple Fourier model.
       {
         fouriery <- as.matrix(fourier(seq_along(y), m, 1))
diff --git a/R/lagwalk.R b/R/lagwalk.R
index 4ddac67..d777d0f 100644
--- a/R/lagwalk.R
+++ b/R/lagwalk.R
@@ -1,13 +1,13 @@
 train_lagwalk <- function(.data, specials, ...) {
   if (length(measured_vars(.data)) > 1) {
-    abort("Only univariate responses are supported by lagwalks.")
+    cli::cli_abort("Only univariate responses are supported by lagwalks.")
   }
   
   y <- unclass(.data)[[measured_vars(.data)]]
   n <- length(y)
 
   if (all(is.na(y))) {
-    abort("All observations are missing, a model cannot be estimated without data.")
+    cli::cli_abort("All observations are missing, a model cannot be estimated without data.")
   }
 
   drift <- specials$drift[[1]][[1]] %||% FALSE
@@ -128,7 +128,7 @@ RW <- function(formula, ...) {
           lag <- 1
         }
         if (!rlang::is_integerish(lag)) {
-          warn("Non-integer lag orders for random walk models are not supported. Rounding to the nearest integer.")
+          cli::cli_warn("Non-integer lag orders for random walk models are not supported. Rounding to the nearest integer.")
           lag <- round(lag)
         }
         get_frequencies(lag, self$data, .auto = "smallest")
@@ -167,10 +167,10 @@ SNAIVE <- function(formula, ...) {
       lag = function(lag = NULL) {
         lag <- get_frequencies(lag, self$data, .auto = "smallest")
         if (lag == 1) {
-          abort("Non-seasonal model specification provided, use RW() or provide a different lag specification.")
+          cli::cli_abort("Non-seasonal model specification provided, use {.fn RW} or provide a different lag specification.")
         }
         if (!rlang::is_integerish(lag)) {
-          warn("Non-integer lag orders for random walk models are not supported. Rounding to the nearest integer.")
+          cli::cli_warn("Non-integer lag orders for random walk models are not supported. Rounding to the nearest integer.")
           lag <- round(lag)
         }
         lag
@@ -246,7 +246,7 @@ forecast.RW <- function(object, new_data, specials = NULL, simulate = FALSE, boo
 #' @export
 generate.RW <- function(x, new_data, bootstrap = FALSE, ...) {
   if (!is_regular(new_data)) {
-    abort("Simulation new_data must be regularly spaced")
+    cli::cli_abort("Simulation new_data must be regularly spaced")
   }
 
   lag <- x$lag
@@ -263,7 +263,7 @@ generate.RW <- function(x, new_data, bootstrap = FALSE, ...) {
   future <- fits[start_pos + seq_len(lag) - 1]
 
   if (any(is.na(future))) {
-    abort("The first lag window for simulation must be within the model's training set.")
+    cli::cli_abort("The first lag window for simulation must be within the model's training set.")
   }
 
   if (!(".innov" %in% names(new_data))) {
diff --git a/R/lm.R b/R/lm.R
index 77870ef..c4c802d 100644
--- a/R/lm.R
+++ b/R/lm.R
@@ -272,7 +272,7 @@ forecast.TSLM <- function(object, new_data, specials = NULL, bootstrap = FALSE,
   xreg <- specials$xreg[[1]]
 
   if (rank < ncol(xreg)) {
-    warn("prediction from a rank-deficient fit may be misleading")
+    cli::cli_warn("prediction from a rank-deficient fit may be misleading")
   }
 
   # Intervals
diff --git a/R/mean.R b/R/mean.R
index 6467dad..e112511 100644
--- a/R/mean.R
+++ b/R/mean.R
@@ -1,13 +1,13 @@
 #' @importFrom stats sd
 train_mean <- function(.data, specials, ...) {
   if (length(measured_vars(.data)) > 1) {
-    abort("Only univariate responses are supported by MEAN.")
+    cli::cli_abort("Only univariate responses are supported by MEAN.")
   }
 
   y <- unclass(.data)[[measured_vars(.data)]]
 
   if (all(is.na(y))) {
-    abort("All observations are missing, a model cannot be estimated without data.")
+    cli::cli_abort("All observations are missing, a model cannot be estimated without data.")
   }
 
   n <- length(y)
@@ -287,10 +287,10 @@ refit.model_mean <- function(object, new_data, specials = NULL, reestimate = FAL
   y <- unclass(new_data)[[measured_vars(new_data)]]
   
   if (all(is.na(y))) {
-    abort("All new observations are missing, model cannot be applied.")
+    cli::cli_abort("All new observations are missing, model cannot be applied.")
   }
 
-  if (!is_null(specials$window)) warn("A rolling mean model cannot be refitted, the most recent mean from the fitted model will be used as a fixed estimate of the mean.")
+  if (!is_null(specials$window)) cli::cli_warn("A rolling mean model cannot be refitted, the most recent mean from the fitted model will be used as a fixed estimate of the mean.")
   
   n <- length(y)
 
diff --git a/R/nnetar.R b/R/nnetar.R
index 5d55c5b..2e5097d 100644
--- a/R/nnetar.R
+++ b/R/nnetar.R
@@ -1,21 +1,21 @@
 #' @importFrom stats ar complete.cases
 train_nnetar <- function(.data, specials, n_nodes, n_networks, scale_inputs, wts = NULL,...) {
-  require_package("nnet")
+  check_installed("nnet")
 
   if (length(measured_vars(.data)) > 1) {
-    abort("Only univariate responses are supported by NNETAR.")
+    cli::cli_abort("Only univariate responses are supported by NNETAR.")
   }
 
   y <- x <- unclass(.data)[[measured_vars(.data)]]
 
   if (all(is.na(y))) {
-    abort("All observations are missing, a model cannot be estimated without data.")
+    cli::cli_abort("All observations are missing, a model cannot be estimated without data.")
   }
 
   n <- length(x)
 
   if (n < 3) {
-    stop("Not enough data to fit a model")
+    cli::cli_abort("Not enough data to fit a model")
   }
 
   # Get args
@@ -26,21 +26,21 @@ train_nnetar <- function(.data, specials, n_nodes, n_networks, scale_inputs, wts
   # Check for constant data in time series
   constant_data <- is.constant(x)
   if (constant_data) {
-    warn("Constant data, setting `AR(p=1, P=0)`, and `scale_inputs=FALSE`")
+    cli::cli_warn("Constant data, setting {.code AR(p=1, P=0)}, and {.code scale_inputs=FALSE}")
     scale_inputs <- FALSE
     p <- 1
     P <- 0
   }
   # Check for insufficient data for seasonal lags
   if (P > 0 && n < period * P + 2) {
-    warn("Series too short for seasonal lags")
+    cli::cli_warn("Series too short for seasonal lags")
     P <- 0
   }
   # Check for constant data in xreg
   if (!is.null(xreg)) {
     xreg <- as.matrix(xreg)
     if (any(apply(xreg, 2, is.constant))) {
-      warn("Constant xreg column, setting `scale_inputs=FALSE`")
+      cli::cli_warn("Constant xreg column, setting {.code scale_inputs=FALSE}")
       scale_inputs <- FALSE
     }
   }
@@ -87,7 +87,7 @@ train_nnetar <- function(.data, specials, n_nodes, n_networks, scale_inputs, wts
     p <- max(length(ar(y_sa, na.action=stats::na.pass)$ar), 1)
   }
   if (p >= n) {
-    warn("Reducing number of lagged inputs due to short series")
+    cli::cli_warn("Reducing number of lagged inputs due to short series")
     p <- n - 1
   }
   lags <- 1:p
@@ -115,7 +115,7 @@ train_nnetar <- function(.data, specials, n_nodes, n_networks, scale_inputs, wts
   x <- x[j]
   ## Stop if there's no data to fit
   if (NROW(x_lags) == 0) {
-    abort("No data to fit (possibly due to missing values)")
+    cli::cli_abort("No data to fit (possibly due to missing values)")
   }
 
   # Fit the nnet and consider the Wts argument for nnet::nnet() if provided: 
@@ -168,7 +168,7 @@ specials_nnetar <- new_specials(
     period <- get_frequencies(period, self$data, .auto = "smallest")
     if (period == 1) {
       if (!missing(P) && P > 0) {
-        warn("Non-seasonal data, ignoring seasonal lags")
+        cli::cli_warn("Non-seasonal data, ignoring seasonal lags")
       }
       P <- 0
     }
@@ -274,7 +274,7 @@ NNETAR <- function(formula, n_nodes = NULL, n_networks = 20, scale_inputs = TRUE
 #'   forecast(times = 10)
 #' @export
 forecast.NNETAR <- function(object, new_data, specials = NULL, simulate = TRUE, bootstrap = FALSE, times = 5000, ...) {
-  require_package("nnet")
+  check_installed("nnet")
 
   # Prepare xreg
   xreg <- specials$xreg[[1]]
@@ -293,7 +293,7 @@ forecast.NNETAR <- function(object, new_data, specials = NULL, simulate = TRUE,
 
   # Compute forecast intervals
   if (!simulate) {
-    warn("Analytical forecast distributions are not available for NNETAR.")
+    cli::cli_warn("Analytical forecast distributions are not available for NNETAR.")
     times <- 0
   }
   sim <- map(seq_len(times), function(x) {
@@ -313,7 +313,7 @@ forecast.NNETAR <- function(object, new_data, specials = NULL, simulate = TRUE,
     {
       fcdata <- c(future_lags[lags], xreg[i, ])
       if (any(is.na(fcdata))) {
-        abort("I can't use NNETAR to forecast with missing values near the end of the series.")
+        cli::cli_abort("I can't use NNETAR to forecast with missing values near the end of the series.")
       }
       fc[i] <- mean(map_dbl(object$model, predict, newdata = fcdata))
       future_lags <- c(fc[i], future_lags[-maxlag])
@@ -379,7 +379,7 @@ generate.NNETAR <- function(x, new_data, specials = NULL, bootstrap = FALSE, ...
     {
       fcdata <- c(future_lags[lags], xreg[i, ])
       if (any(is.na(fcdata))) {
-        abort("I can't use NNETAR to forecast with missing values near the end of the series.")
+        cli::cli_abort("I can't use NNETAR to forecast with missing values near the end of the series.")
       }
       path[i] <- mean(map_dbl(x$model, predict, newdata = fcdata)) + e[i]
       future_lags <- c(path[i], future_lags[-maxlag])
diff --git a/R/theta.R b/R/theta.R
index 8345bdf..0ace4d7 100644
--- a/R/theta.R
+++ b/R/theta.R
@@ -1,14 +1,14 @@
 #' @importFrom stats sd
 train_theta <- function(.data, specials, ...) {
   if (length(measured_vars(.data)) > 1) {
-    abort("Only univariate responses are supported by the Theta method")
+    cli::cli_abort("Only univariate responses are supported by the Theta method")
   }
   
   y <- unclass(.data)[[measured_vars(.data)]]
   n <- length(y)
   
   if (all(is.na(y))) {
-    abort("All observations are missing, a model cannot be estimated without data.")
+    cli::cli_abort("All observations are missing, a model cannot be estimated without data.")
   }
   
   # Check seasonality
@@ -28,7 +28,7 @@ train_theta <- function(.data, specials, ...) {
   if (m > 1L) {
     dcmp <- stats::decompose(y, type = specials$season[[1]]$method)
     if (any(abs(dcmp$seasonal) < 1e-4)) {
-      warning("Seasonal indexes equal to zero. Using non-seasonal Theta method")
+      cli::cli_warn("Seasonal indexes equal to zero. Using non-seasonal Theta method")
     } else {
       y_sa <- if(dcmp$type == "additive") dcmp$x - dcmp$seasonal else dcmp$x / dcmp$seasonal
     }
@@ -146,7 +146,7 @@ THETA <- function(formula, ...) {
 #' @export
 forecast.fable_theta <- function(object, new_data, specials = NULL, bootstrap = FALSE, times = 5000, ...) {
   if (bootstrap) {
-    abort("Bootstrapped forecasts are not yet supported for the Theta method.")
+    cli::cli_abort("Bootstrapped forecasts are not yet supported for the Theta method.")
   }
   h <- NROW(new_data)
   
diff --git a/R/utils.R b/R/utils.R
index 0c2adc2..363d9ae 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -6,17 +6,9 @@ is.constant <- function(x) {
 
 assignSpecials <- function(x, env = caller_env()) {
   imap(x, function(.x, nm) {
-    if (length(.x) > 1) warn(sprintf("Only one special for `%s` is allowed, defaulting to the first usage", nm))
+    if (length(.x) > 1) cli::cli_warn("Only one special for {.code {nm}} is allowed, defaulting to the first usage")
     imap(.x[[1]], function(.x, .y) assign(.y, .x, envir = env))
   })
 }
 
-require_package <- function(pkg) {
-  if (!requireNamespace(pkg, quietly = TRUE)) {
-    abort(
-      sprintf('The `%s` package must be installed to use this functionality. It can be installed with install.packages("%s")', pkg, pkg)
-    )
-  }
-}
-
 `%||%` <- function(x, y) if (is_null(x)) y else x
diff --git a/R/var.R b/R/var.R
index 615ce57..b8ae899 100644
--- a/R/var.R
+++ b/R/var.R
@@ -84,7 +84,7 @@ estimate_var <- function(y, p, xreg, constant) {
 specials_var <- new_specials(
   AR = function(p = 0:5) {
     if (any(p < 0)) {
-      warn("The AR order must be non-negative. Only non-negative orders will be considered.")
+      cli::cli_warn("The AR order must be non-negative. Only non-negative orders will be considered.")
       p <- p[p >= 0]
     }
     list(p = p)
@@ -200,7 +200,7 @@ VAR <- function(formula, ic = c("aicc", "aic", "bic"), ...) {
 forecast.VAR <- function(object, new_data = NULL, specials = NULL,
                          bootstrap = FALSE, times = 5000, ...) {
   if (bootstrap) {
-    abort("Bootstrapped forecasts for VARs are not yet implemented.")
+    cli::cli_abort("Bootstrapped forecasts for VARs are not yet implemented.")
   }
 
   h <- NROW(new_data)
@@ -468,4 +468,4 @@ IRF.VAR <- function(x, new_data, specials, impulse = NULL, orthogonal = FALSE, .
   irf[colnames(x$coef)] <- split(irf$.sim, col(irf$.sim))
   irf$.innov <- irf$.sim <- NULL
   irf
-}
\ No newline at end of file
+}
diff --git a/R/vecm.R b/R/vecm.R
index 2d56a63..e292362 100644
--- a/R/vecm.R
+++ b/R/vecm.R
@@ -130,7 +130,7 @@ estimate_vecm <- function(y, p, sr_xreg, constant, r, ...) {
 specials_vecm <- new_specials(
   AR = function(p = 0:5) {
     if (any(p < 0)) {
-      warn("The AR order must be non-negative. Only non-negative orders will be considered.")
+      cli::cli_warn("The AR order must be non-negative. Only non-negative orders will be considered.")
       p <- p[p >= 0]
     }
     list(p = p)
@@ -446,4 +446,4 @@ IRF.VECM <- function(x, new_data, specials, impulse = NULL, orthogonal = FALSE,
   irf[colnames(x$coef)] <- split(irf$.sim, col(irf$.sim))
   irf$.innov <- irf$.sim <- NULL
   irf
-}
\ No newline at end of file
+}