diff --git a/R/class_clv_fitted.R b/R/class_clv_fitted.R
index 1a168c95..d3042107 100644
--- a/R/class_clv_fitted.R
+++ b/R/class_clv_fitted.R
@@ -49,3 +49,63 @@ setClass(Class = "clv.fitted", contains = "VIRTUAL",
            optimx.hessian           = matrix(data = numeric(0))))
 
 
+
+
+clv.fitted.get.LL <- function(clv.fitted){
+
+  # Calling the LL with the exact same inputs/specification as when the fitting
+  # is not trivial as there are plenty of options.
+  # To reproduce it, the exact steps of the estimation are repeated here.
+
+  # Start parameters are not really required, they are just stored as item
+  # `par` for `optimx()`
+  final.coefs <- drop(tail(coef(clv.fitted@optimx.estimation.output), n=1))
+
+  prepared.optimx.args <- clv.controlflow.estimate.prepare.optimx.args(
+    clv.fitted=clv.fitted,
+    start.params.all= final.coefs)
+
+  prepared.optimx.args <- clv.model.prepare.optimx.args(
+    clv.model=clv.fitted@clv.model,
+    clv.fitted=clv.fitted,
+    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. This also removes the LL itself (`fn`).
+  names.optimx.args <- setdiff(formalArgs(optimx), "...")
+  call.args <- prepared.optimx.args[!(names(prepared.optimx.args) %in% names.optimx.args)]
+
+  # Could save memory as the returned method is a closure and has the environment
+  # in which is was defined attached. Hence all variables in this method here
+  # which may be large. However, it can also be useful to have these objects.
+
+  # Wrapper to call the LL with the original args.
+  # It is preferred to return a method rather than calling it immediately
+  # because generating the call args may take time.
+  LL <- function(params){
+    req.names <- call.args$LL.param.names.to.optimx
+
+    # Ensure named same as original
+    if(!(identical(sort(names(params)), sort(req.names)))){
+      check_err_msg(paste0(
+        "'params' has to be named ",
+        paste(req.names, collapse = ", "),
+        ". Often, `drop(coef(model@optimx.estimation.output))` is useful."))
+    }
+
+    # Ensure name and position are the same as original order. This is required as
+    # any param input will be re-named by slapping `LL.param.names.to.optimx` on
+    # them (using names() <- ) in the first interlayer and then further accessed
+    # by name.
+    call.args$LL.params <- params[req.names] # bring to correct order
+    return(do.call(what = prepared.optimx.args$fn, args = call.args))
+  }
+
+  return(LL)
+}
diff --git a/R/f_generics_clvfitted.R b/R/f_generics_clvfitted.R
index 08e05676..7e2a1423 100644
--- a/R/f_generics_clvfitted.R
+++ b/R/f_generics_clvfitted.R
@@ -35,3 +35,4 @@ setMethod("clv.fitted.estimate.same.specification.on.new.data", signature = "clv
   new.fitted@call <- cl
   return(new.fitted)
 })
+
diff --git a/R/f_interface_hessian.R b/R/f_interface_hessian.R
index 9f9dcfab..08ae213d 100644
--- a/R/f_interface_hessian.R
+++ b/R/f_interface_hessian.R
@@ -25,52 +25,29 @@ 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.
-
+  # Get final parameters
   # To get coefficients at the same scale (log-transformed) and names as when
-  # estimating the model, get them from the optimx results.
+  # estimating the model, get them directly from the optimx output.
   # Cannot use coef() as the reported parameters are transformed back and named
-  # differently
+  # 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.
+  # Get LL
+  fn.LL <- clv.fitted.get.LL(object)
 
-  # 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))
-  }
+  # fn.call.LL.named <- function(par){
+  #   names(par) <- names(final.coefs)
+  #   return(fn.LL(par))
+  # }
 
   # 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,
+    func=fn.LL,
     x=final.coefs,
     method="Richardson",
     method.args = method.args)
diff --git a/tests/testthat/test_correctness_hessian.R b/tests/testthat/test_correctness_hessian.R
new file mode 100644
index 00000000..57975e07
--- /dev/null
+++ b/tests/testthat/test_correctness_hessian.R
@@ -0,0 +1,134 @@
+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)){
+    # Default specification
+    fn.compare.hessian(
+      fit.apparel.static(model=m, optimx.args = optimx.args)
+    )
+
+    # With constrained covs
+    fn.compare.hessian(
+      fit.apparel.static(
+        model=m,
+        names.cov.constr = "Gender",
+        optimx.args = optimx.args)
+    )
+
+    # With regularization
+    fn.compare.hessian(
+      fit.apparel.static(
+        model=m,
+        reg.lambdas = c(life = 10, trans=5),
+        optimx.args = optimx.args)
+    )
+
+    # With constrained covs & regularization
+    fn.compare.hessian(
+      fit.apparel.static(
+        model=m,
+        names.cov.constr = "Channel",
+        reg.lambdas = c(life = 10, trans=5),
+        optimx.args = optimx.args)
+    )
+  }
+
+
+  # PNBD only: 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()
+
+  # Default
+  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)
+  )
+
+  # With constrained covs
+  fn.compare.hessian(
+    fit.apparel.dyncov(model = pnbd, names.cov.constr = "Gender", optimx.args = optimx.args)
+  )
+
+  # With regularization
+  fn.compare.hessian(
+    fit.apparel.dyncov(model = pnbd, reg.lambdas = c(trans=10, life=5), optimx.args = optimx.args)
+  )
+
+})
+
+
+
+test_that("hessian() fails if parameters are non-finite",{
+  skip_on_cran()
+
+  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")
+})
+
+test_that("Internal clv.fitted.get.LL: Params position and order checked", {
+  # Indirectly tested for correctness by being used in hessian()
+  skip_on_cran()
+
+  p.reg.constr <- fit.apparel.static(
+    model = pnbd,
+    reg.lambdas = c(trans = 4, life = 9),
+    names.cov.constr = "Gender",
+    optimx.args= optimx.args)
+
+  LL.reg.constr <- clv.fitted.get.LL(p.reg.constr)
+  final.coefs <- drop(coef(p.reg.constr@optimx.estimation.output))
+
+  # Have to be named
+  expect_error(LL.reg.constr(setNames(final.coefs, NULL)), regexp = "has to be named")
+
+  # Does not work with extra coefs
+  expect_error(LL.reg.constr(coef(p.reg.constr)), regexp = "has to be named")
+
+  # Results are independent of order
+  expect_identical(
+    LL.reg.constr(sort(final.coefs, decreasing = FALSE)),
+    LL.reg.constr(sort(final.coefs, decreasing = TRUE))
+  )
+
+})
diff --git a/tests/testthat/tests_correctness_hessian.R b/tests/testthat/tests_correctness_hessian.R
deleted file mode 100644
index 9e09e7a2..00000000
--- a/tests/testthat/tests_correctness_hessian.R
+++ /dev/null
@@ -1,67 +0,0 @@
-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")
-})