diff --git a/DESCRIPTION b/DESCRIPTION
index fad0ad9c..128217a7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -122,6 +122,7 @@ Collate:
     'f_interface_clvdata.R'
     'f_interface_gg.R'
     'f_interface_ggomnbd.R'
+    'f_interface_hessian.R'
     'f_interface_latentattrition.R'
     'f_interface_lrtest.R'
     'f_interface_newcustomer.R'
diff --git a/NAMESPACE b/NAMESPACE
index 1b983919..6dc2c7d6 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -13,6 +13,7 @@ S3method(lmtest::lrtest,clv.fitted)
 S3method(logLik,clv.fitted)
 S3method(nobs,clv.data)
 S3method(nobs,clv.fitted)
+S3method(numDeriv::hessian,clv.fitted)
 S3method(plot,clv.data)
 S3method(plot,clv.fitted.spending)
 S3method(plot,clv.fitted.transactions)
@@ -53,6 +54,7 @@ exportMethods(bgbb)
 exportMethods(bgnbd)
 exportMethods(gg)
 exportMethods(ggomnbd)
+exportMethods(hessian)
 exportMethods(lrtest)
 exportMethods(plot)
 exportMethods(pmf)
diff --git a/R/f_interface_hessian.R b/R/f_interface_hessian.R
new file mode 100644
index 00000000..9f9dcfab
--- /dev/null
+++ b/R/f_interface_hessian.R
@@ -0,0 +1,106 @@
+
+#' @name hessian
+#' @title Calculate hessian for a fitted model
+#'
+#' @description Calculate a numerical approximation to the Hessian matrix at
+#' the final estimated parameters using \code{numDeriv::hessian}.
+#'
+#' @param object Fitted model
+#' @param method.args List of options forwarded to the numerical approximation
+#' method. See \link[numDeriv:hessian]{numDeriv::hessian}.
+#' @template template_param_dots
+#'
+#' @returns The hessian matrix, with column and row names set to the parameter
+#' names used to call the LL.
+NULL
+
+
+#' @rdname hessian
+#' @importFrom numDeriv hessian
+#' @importFrom utils tail
+#' @importFrom methods formalArgs
+#' @exportS3Method numDeriv::hessian
+hessian.clv.fitted <- function(object, method.args = list()){
+
+  # Register for dispatch on a method defined in another package by using
+  # @exportS3Method which adds `S3method(numDeriv::hessian,clv.fitted)` to NAMESPACE
+
+  # To calculate the Hessian, the LL needs to be called with the final parameters
+  # Calling the LL with the exact same inputs/specification as when the fitting
+  # is however not trivial as there are plenty of options.
+  # To reproduce it, the exact steps of the estimation are repeated here.
+
+  # To get coefficients at the same scale (log-transformed) and names as when
+  # estimating the model, get them from the optimx results.
+  # Cannot use coef() as the reported parameters are transformed back and named
+  # differently
+  final.coefs <- drop(tail(coef(object@optimx.estimation.output), n=1))
+
+  if(anyNA(final.coefs)){
+    check_err_msg("Cannot proceed because there are NAs in the estimated coefficients!")
+  }
+
+  prepared.optimx.args <- clv.controlflow.estimate.prepare.optimx.args(
+    clv.fitted=object,
+    start.params.all=final.coefs)
+
+  prepared.optimx.args <- clv.model.prepare.optimx.args(
+    clv.model=object@clv.model,
+    clv.fitted=object,
+    prepared.optimx.args=prepared.optimx.args)
+
+  prepared.optimx.args[["LL.param.names.to.optimx"]] <- names(prepared.optimx.args$par)
+
+  # In the estimation procedure, the user can also supply custom `optimx.args`
+  # which override `prepared.optimx.args` here. Because optimx is not called
+  # here, there is no need to add `optimx.args` here.
+
+  # The generated args also contain parameters for optimx. These are not required
+  # for the LL and need to be removed.
+  names.optimx.args <- setdiff(formalArgs(optimx), "...")
+  call.args <- prepared.optimx.args[!(names(prepared.optimx.args) %in% names.optimx.args)]
+
+
+  # Wrapper to call the LL with the original args and `par` given by numDeriv
+  fn.with.call.args <- function(par){
+    call.args$LL.params <- par
+    return(do.call(prepared.optimx.args$fn, call.args))
+  }
+
+  # Have to refer to `numDeriv` namespace directly (`::`) as `hessian()` would
+  # dispatch to `CLVTools::hessian` and fail if `numDeriv` is not attached
+  H <- numDeriv::hessian(
+    func=fn.with.call.args,
+    x=final.coefs,
+    method="Richardson",
+    method.args = method.args)
+  # Names as in optimx (log-scale etc)
+  colnames(H) <- rownames(H) <- names(final.coefs)
+
+  return(H)
+}
+
+
+# In order to be able to use `hessian()` without having `numDeriv` loaded or even
+# installed, define and export `hessian()` as a generic in CLVTools.
+# The S4 generic is NOT defined with the exact same signature as the
+# S3 generic `numDeriv::hessian <- function(func, x, method, method.args, ...){ ... }`.
+#
+# The numDeriv package exports an S3 generic `hessian()` what masks the generic (whether S3
+# or S4) exported by CLVTools if the numDeriv package is loaded after CLVTools.
+# Therefore, define and export also as a S3 method `CLVTools::hessian.clv.fitted`.
+#
+# ?Methods_for_Nongenerics on dispatching an S4 object to S3 generics method in
+# another package: Recommends to define both methods: The S3 method and also supply
+# the identical function as the definition of the S4 method.
+#' @rdname hessian
+#' @exportMethod hessian
+setGeneric(name = "hessian", def=function(object, ...)
+  standardGeneric("hessian"))
+
+
+#' @rdname hessian
+#' @include all_generics.R
+#' @exportMethod hessian
+setMethod("hessian", signature(object="clv.fitted"), definition = hessian.clv.fitted)
+
diff --git a/man/hessian.Rd b/man/hessian.Rd
new file mode 100644
index 00000000..16454dbc
--- /dev/null
+++ b/man/hessian.Rd
@@ -0,0 +1,30 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/f_interface_hessian.R
+\name{hessian}
+\alias{hessian}
+\alias{hessian.clv.fitted}
+\alias{hessian,clv.fitted-method}
+\title{Calculate hessian for a fitted model}
+\usage{
+\method{hessian}{clv.fitted}(object, method.args = list())
+
+hessian(object, ...)
+
+\S4method{hessian}{clv.fitted}(object, method.args = list())
+}
+\arguments{
+\item{object}{Fitted model}
+
+\item{method.args}{List of options forwarded to the numerical approximation
+method. See \link[numDeriv:hessian]{numDeriv::hessian}.}
+
+\item{...}{Ignored}
+}
+\value{
+The hessian matrix, with column and row names set to the parameter
+names used to call the LL.
+}
+\description{
+Calculate a numerical approximation to the Hessian matrix at
+the final estimated parameters using \code{numDeriv::hessian}.
+}
diff --git a/tests/testthat/helper_arrange.R b/tests/testthat/helper_arrange.R
index 938e3649..23c2b2a7 100644
--- a/tests/testthat/helper_arrange.R
+++ b/tests/testthat/helper_arrange.R
@@ -113,9 +113,9 @@ fit.cdnow <- function(
     estimation.split = 37,
     name.price = 'Price',
     model = pnbd,
-    start.params.model = c(),
     verbose = FALSE,
-    optimx.args = list()) {
+    ...
+    ) {
 
   clv.cdnow <- fct.helper.create.clvdata.cdnow(
     data.cdnow = data.cdnow,
@@ -127,9 +127,8 @@ fit.cdnow <- function(
     what = model,
     args = list(
       clv.data = clv.cdnow,
-      start.params.model = start.params.model,
-      optimx.args = optimx.args,
-      verbose = verbose
+      verbose = verbose,
+      ...
     )
   ))
 }
@@ -142,7 +141,6 @@ fit.apparel.nocov <- function(
     model = pnbd,
     verbose=FALSE,
     # start.params.model = c(),
-    # verbose = FALSE,
     # optimx.args = list()
     ...
 ) {
diff --git a/tests/testthat/tests_correctness_hessian.R b/tests/testthat/tests_correctness_hessian.R
new file mode 100644
index 00000000..9e09e7a2
--- /dev/null
+++ b/tests/testthat/tests_correctness_hessian.R
@@ -0,0 +1,67 @@
+skip_on_cran()
+
+optimx.args <- list(itnmax=100)
+
+fn.compare.hessian <- function(clv.fitted){
+  expect_equal(
+    hessian(clv.fitted),
+    clv.fitted@optimx.hessian
+    )
+}
+
+
+test_that("hessian() produces same result - no cov", {
+  skip_on_cran()
+
+  for(m in list(pnbd, bgnbd, ggomnbd, gg)){
+    fn.compare.hessian(
+      fit.cdnow(model=m, optimx.args = optimx.args)
+    )
+  }
+
+  # With cor
+  fn.compare.hessian(
+    fit.cdnow(model = pnbd, use.cor=TRUE, optimx.args = optimx.args)
+  )
+})
+
+
+
+test_that("hessian() produces same result - static cov", {
+  skip_on_cran()
+
+  for(m in list(pnbd, bgnbd, ggomnbd, gg)){
+    fn.compare.hessian(
+      fit.apparel.static(model=m, optimx.args = optimx.args)
+    )
+  }
+
+  # With cor
+  fn.compare.hessian(
+    fit.apparel.static(model = pnbd, use.cor=TRUE, optimx.args = optimx.args)
+  )
+})
+
+
+test_that("hessian() produces same result - dyn cov", {
+  skip_on_cran()
+
+  # No cor
+  fn.compare.hessian(
+    fit.apparel.dyncov(model = pnbd, optimx.args = optimx.args)
+  )
+
+  # With cor
+  fn.compare.hessian(
+    fit.apparel.dyncov(model = pnbd, use.cor=TRUE, optimx.args = optimx.args)
+  )
+})
+
+
+
+test_that("hessian() fails if parameters are non-finite",{
+  p.cdnow <- fit.cdnow(optimx.args=optimx.args)
+  p.cdnow@optimx.estimation.output[1, "log.r"] <- NA_real_
+
+  expect_error(hessian(p.cdnow), regexp = "Cannot proceed")
+})