diff --git a/.gitignore b/.gitignore
index 8735c5f..7065f16 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,8 +7,6 @@ inst/include
 src/*.o
 src/*.so
 src/*.dll
-R/RcppExports.R
-src/RcppExports.cpp
 *.html
 *.RData
 *.bak
diff --git a/R/RcppExports.R b/R/RcppExports.R
new file mode 100644
index 0000000..71634cb
--- /dev/null
+++ b/R/RcppExports.R
@@ -0,0 +1,635 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#' @name kf
+#' @rdname kf
+#' @export
+kf_cpp <- function(A, C, Q, R, S, y_t, P1, a1) {
+    .Call(`_RLDM_kf_cpp`, A, C, Q, R, S, y_t, P1, a1)
+}
+
+#' @name kf
+#' @rdname kf
+#' @export
+kf2_cpp <- function(A, C, H_t, y_t, P1_R, a1) {
+    .Call(`_RLDM_kf2_cpp`, A, C, H_t, y_t, P1_R, a1)
+}
+
+#' @name ll_kf
+#' @rdname ll_kf
+#' @export
+ll_kf_cpp <- function(A, C, Q, R, S, y_t, P1, a1, tol) {
+    .Call(`_RLDM_ll_kf_cpp`, A, C, Q, R, S, y_t, P1, a1, tol)
+}
+
+#' @name ll_kf
+#' @rdname ll_kf
+#' @export
+ll_kf2_cpp <- function(A, C, H_t, y_t, P1_R, a1, tol) {
+    .Call(`_RLDM_ll_kf2_cpp`, A, C, H_t, y_t, P1_R, a1, tol)
+}
+
+#' Compute the log likelihood for a statespace system
+#' described by a model template.
+#'
+#' This is an internal helper function, used by the function factory \code{\link{ll_FUN}}. For a more detailed
+#' documentation of the log Likelihood, see \code{\link{ll_kf}}.
+#'
+#' @param theta \eqn{(K)} dimensional vector of "deep" parameters.
+#' @param y \eqn{(m,N)} matrix with the observed outputs:
+#'        \eqn{(y_1,y_2,\ldots,y_N)}{(y[1],y[2],...,y[N])}.
+#' @param SYS \eqn{(m+s,m+s)} matrix, is overwritten with the system matrix
+#'        \eqn{[A,B | C,D]}.
+#' @param H_SYS \eqn{(m+s)^2, K)} matrix.
+#' @param h_SYS \eqn{((m+s)^2)}-dimensional vector. Note that \code{vec(SYS) = H_SYS*theta + h_SYS}.
+#' @param sigma_L \eqn{(m,m)} matrix, is overwritten with the left square root of the
+#'        noise covariance matrix.
+#' @param H_sigma_L \eqn{(m^2, K)} matrix.
+#' @param h_sigma_L \eqn{(m^2)}-dimensional vector. Note that
+#'        \code{vec(sigma_L) = H_sigma_L*theta + h_sigma_L}.
+#' @param VAR \eqn{(m+s,m+s)} matrix, is overwritten with the
+#'        covariance matrix \eqn{[Q,S | S',R] = [B | C] sigma_L sigma_L' [B', C']}
+#' @param P1 \eqn{(s,s)} matrix, is overwritten with the
+#'        initial state covariance matrix (computed via a Lyapunov equation).
+#' @param tol (double) tolerance used by ll_kf_cpp.
+#' @param err (double) return err, if the computation of P1 fails.
+#' @export
+ll_kf_theta_cpp <- function(theta, y, SYS, H_SYS, h_SYS, sigma_L, H_sigma_L, h_sigma_L, VAR, P1, tol, err) {
+    .Call(`_RLDM_ll_kf_theta_cpp`, theta, y, SYS, H_SYS, h_SYS, sigma_L, H_sigma_L, h_sigma_L, VAR, P1, tol, err)
+}
+
+#' Outputs of a statespace system
+#'
+#' @description
+#' This internal helper function computes the outputs and states
+#' for a statespace system of the form
+#' \deqn{a_{t+1} = A a_t + B u_t, \; y_t   = C a_t + D u_t}{
+#'       a[t+1] = A a[t] + B u[t] and y[t]   = C a[t] + D u[t]}
+#'
+#' @note
+#' Use this procedure with care!
+#'
+#' * The procedure does \bold{not} check the input arguments.
+#' * The procedure \bold{overwrites} the input arguments \code{a} and \code{u}.
+#' * The data matrices are organized columnwise (to avoid memory shuffling)!
+#'
+#' @param A \eqn{(s,s)} matrix.
+#' @param B \eqn{(s,n)} matrix.
+#' @param C \eqn{(m,s)} matrix.
+#' @param D \eqn{(m,n)} matrix.
+#' @param u \eqn{(n,N)} matrix with the inputs/disturbances:
+#'          \eqn{(u_1,u_2,\ldots,u_N)}{(u[1],u[2],...,u[N])}.
+#' @param a \eqn{(s,N+1)} matrix. This matrix is overwritten with the (computed) states:
+#'            \eqn{(a_1,a_2,\ldots,a_N,a_{N+1})}{(a[1],a[2],\ldots,a[N],a[N+1])}.
+#'            On input \code{a[,1]} must hold the initial state \eqn{a_1}{a[1]}.
+#' @param y \eqn{(m,N)} matrix. This matrix is overwritten with (computed) outputs:
+#'           \eqn{(y_1,y_2,\ldots,y_N)}{(y[1],y[2],...,y[N])}.
+#'
+#' @return This RcppArmadillo routine returns \code{NULL} but \bold{overwrites} the input
+#'         arguments \code{a} and \code{u}!
+#' @export
+#'
+#' @rdname outputs_STSP_cpp
+#' @name outputs_STSP_cpp
+#'
+#' @seealso \code{\link{outputs_ARMA_cpp}}, \code{\link{residuals_ARMA_cpp}},  \code{\link{cll_theta_ARMA_cpp}},
+#'    \code{\link{outputs_STSP_cpp}}, \code{\link{residuals_STSP_cpp}},  \code{\link{cll_theta_STSP_cpp}} and
+#'    \code{\link{solve_de}}, \code{\link{solve_inverse_de}} and \code{\link{ll}}.
+#'
+#' @examples
+#' # generate a random statespace model (3 outputs, 2 inputs and 4 states)
+#' m = 3
+#' n = 2
+#' s = 4
+#' model = test_stspmod(dim = c(m, n), s = s, digits = 2)
+#'
+#' # generate random noise sequence (sample size N = 10)
+#' n.obs = 10
+#' u = matrix(rnorm(n.obs*n), nrow = n, ncol = n.obs)
+#' print(u)
+#'
+#' # generate matrix for the state sequence
+#' a = matrix(0, nrow = s, ncol = n.obs+1)
+#' a[,1] = rnorm(s) # random initial state a[1]
+#' print(a)
+#'
+#' # generate matrix for the outputs
+#' y = matrix(0, nrow = m, ncol = n.obs)
+#'
+#' # call outputs_STSP_cpp()
+#' outputs_STSP_cpp(model$sys$A, model$sys$B, model$sys$C, model$sys$D, u, a, y)
+#' print(u)
+#' print(a)  # a is overwritten with the computed states
+#' print(y)  # y is overwritten with the computed outputs
+NULL
+
+#' Forward-backward solution of statespace systems
+#'
+#' @description
+#' **DEPRECATED**?
+#' This internal helper function computes the outputs of an, in general **unstable**, statespace system
+#' \deqn{a_{t+1} = A a_t + B u_t, \; y_t = C a_t + D u_t}{
+#'       a[t+1] = A a[t] + B u[t] and y[t] = C a[t] + D u[t]}
+#' by forward and backward recursion. The procedure assumes that the state transition matrix
+#' \eqn{A} is block upper triangular, where the upper block \eqn{A_{11}}{A[11]} is stable (i.e.
+#' all eigenvalues have moduli less than one) and the
+#' lower block \eqn{A_{22}}{A[22]} is unstable (i.e. all eigenvalues have moduli larger than one).
+#' This function is mainly used in the routine \code{\link{innovation_form}}.
+#'
+#' @note
+#' Use this procedure with care!
+#'
+#' * The procedure does \bold{not} check the input arguments. We require \eqn{m > 0}, \eqn{n > 0}.
+#'   Furthermore it is assumed that the state transition matrix \eqn{A} is block upper triangular,
+#'   as explained above.
+#' * The procedure \bold{overwrites} the input arguments \code{y}, \code{as} and \code{au}.
+#' * The data matrices are organized columnwise (to avoid memory shuffling)!
+#'
+#' @param A \eqn{(s, s)} matrix.
+#' @param B \eqn{(s, n)} matrix.
+#' @param C \eqn{(m, s)} matrix.
+#' @param D \eqn{(m, n)} matrix.
+#' @param u \eqn{(n, N)} matrix with the inputs \eqn{(u_1,...,u_N}{(u[1],...,u[N])}.
+#' @param au \eqn{(su,N+1)} matrix. This matrix is **overwritten** with the (computed)
+#'           states of the unstable part of the system.
+#'            \eqn{(a_{u1},a_{u2},\ldots,a_{uN},a_{u,N+1})}{(au[1],au[2],\ldots,au[N],au[N+1])}.
+#'            On input \code{au[,N+1]} must hold the "initial" state \eqn{a_{u,N+1}}{au[N+1]}.
+#' @param as \eqn{(ss,N+1)} matrix. This matrix is **overwritten** with the (computed)
+#'           states of the stable part of the system.
+#'            \eqn{(a_{s1},a_{s2},\ldots,a_{sN},a_{s,N+1})}{(as[1],as[2],\ldots,as[N],as[N+1])}.
+#'            On input \code{as[,1]} must hold the "initial" state \eqn{a_{s1}}{as[1]}.
+#' @param y \eqn{(m,N)} matrix. This matrix is **overwritten** with (computed) outputs:
+#'           \eqn{(y_1,y_2,\ldots,y_N)}{(y[1],y[2],...,y[N])}.
+#'
+#' @return This RcppArmadillo routine returns \code{NULL} but \bold{overwrites}
+#'         the input argument \code{y}, \code{au} and \code{as} with the computed outputs
+#'         and states!
+#'
+#' @seealso \code{\link{outputs_STSP_cpp}} and \code{\link{innovation_form}}.
+#'
+#' @export
+#' @name fbsolve_STSP_cpp
+#' @rdname fbsolve_STSP_cpp
+#'
+NULL
+
+#' Simulating Output from an RMFD Model (or obtain residuals)
+#'
+#' This RcppArmadillo function calculates for given inputs \code{data_in}
+#' of dimension \eqn{(n \times nobs)}{(n x nobs)},  where \eqn{nobs} is the sample size,
+#' the outputs of dimension \eqn{m}.
+#' Note that data matrices are "transposed" in the sense that every column
+#' corresponds to one observation because of memory management.
+#' \code{data_out} is thus of dimension \eqn{(m x nobs)}.
+#' This function is intended for internal use and thus arguments are not checked.
+#'
+#' @param poly_inv Matrix of dimension \eqn{(n \times n p)}{(n x n p)},
+#'     representing a square matrix polynomial \eqn{c(z)} with \eqn{c_0}{c[0]} equal to the identity matrix (and therefore not stored).
+#'     The coefficients need to be in \strong{reverse} direction, i.e. \eqn{(c_p, ... , c_1)}{(c[p], ... , c[1])},
+#'     where \eqn{p} denotes the degree of \eqn{c(z)}.
+#' @param poly_fwd Matrix of dimensions \eqn{(m \times n(q+1))}{(m x n(q+1))},
+#'    representing a (possibly tall) matrix polynomial \eqn{d(z)} of dimension
+#'    \eqn{(m \times n)}{(m x n)}, where \eqn{m \geq n}{m \ge n}.
+#'    The coefficient are stored "as usual" and including \eqn{d_0}{d[0]}, i.e. \eqn{(d_0, d_1, ... , d_{q-1}, d_{q})}{(d[0], d[1], ... , d[q-1], d[q])},
+#'    where \eqn{q} denotes the degree of \eqn{d(z)}.
+#' @param data_in Matrix of dimension \eqn{(n \times n_obs)}{(n x n_obs)}, i.e. \eqn{(u_1, ..., u_T)}{(u[1], ..., u[T])}.
+#'     Inputs to the RMFD system.
+#' @param data_out Matrix of dimension \eqn{(m \times n_obs)}{(m x n_obs)}, i.e. \eqn{(y_1, ..., y_T)}{(y[1], ..., t[T])}.
+#'     Outputs of the RMFD system. Initially zero and will be \strong{overwritten}.
+#' @param t0 Integer. Time index from which we should start calculating a solution.
+#'     Usually equal to 1.
+#'
+#' @return \code{data_out} is overwritten with the outputs of the RMFD system.
+#'
+#' @export
+#'
+#' @name solve_rmfd_cpp
+#' @rdname solve_rmfd_cpp
+NULL
+
+#' Residuals of an ARMA system
+#'
+#' This internal helper function computes the residuals and the directional derivatives of the
+#' residuals of an ARMA system of the form
+#' \deqn{a_0 y_t + a_1 y_{t-1} + \cdots + a_p y_{t-p} = b_0 u_t + \cdots + b_q u_{t-q}}{
+#'       a[0] y[t] + a[1] y[t-1] + ... + a[p] y[t-p] = b[0] u[t] + ... + b[q] u[t-q]}
+#'
+#' Values \eqn{y_t}{y[t]}, \eqn{u_t}{u[t]} for \eqn{t\leq 0}{t\le 0} are implicitly set to be zero.
+#' However, by starting the iteration with some \eqn{t_0>1}{t0>1} we can enforce non-zero
+#' initial values.
+#'
+#' @note
+#' Use this procedure with care!
+#'
+#' * The procedure does \bold{not} check the input arguments. We require \eqn{m = n > 0},
+#'   \eqn{p,q \geq 0}{p,q \ge 0} and \eqn{1 \leq t_0 \leq N}{1 \le t0  \le N}.
+#' * The procedure \bold{overwrites} the input argument \code{u} (and \code{dU}).
+#' * The data matrices are organized columnwise (to avoid memory shuffling)!
+#' * Note also the non standard representation of the coefficient matrices.
+#'
+#' @param ib0 \eqn{(m, m)} matrix, **inverse** of the coefficient matrix \eqn{b[0]}{b[0]}.
+#' @param B1 \eqn{(m, mq)} matrix, \eqn{-b_0^{-1}(b_q,...,b_1)}{-b[0]^(-1)(b[q],...,b[1])}.
+#' @param A \eqn{(m, n(q+1))} matrix \eqn{b_0^{-1}(a_0,...,a_p}{b[0]^(-1)(a[0],...,a[p])}.
+#' @param t0 integer, start iteration at t = t0.
+#' @param y \eqn{(m, N)} matrix with the observed outputs \eqn{(y_1,...,y_N}{(y[1],...,y[N])}.
+#' @param u \eqn{(m, N)} matrix. This matrix is **overwritten** with the computed
+#'        residuals \eqn{(u_1,...,u_N}{(u[1],...,u[N])}.
+#' @param dU \eqn{(mN, m^2(p+q+2))} matrix or an empty matrix. If non empty then this
+#'        matrix is **overwritten** with the directional derivatives of the vectorized residuals.
+#'        The \eqn{j}-th column of \code{dU} is the derivative of \eqn{vec(u)} with respect
+#'        to the \eqn{j}-th entry of
+#'        \eqn{\mathrm{vec}(a_0,a_1,\ldots,a_p,b_0,\ldots,b_q)}{vec(a[0],a[1],...,a[p],b[0],...,b[q])}
+#'
+#' @return This RcppArmadillo routine returns \code{NULL} but \bold{overwrites}
+#'         the input arguments \code{u} (and \code{dU})!
+#'
+#' @export
+#'
+#' @rdname residuals_ARMA_cpp
+#' @name residuals_ARMA_cpp
+#'
+#' @seealso \code{\link{outputs_ARMA_cpp}}, \code{\link{residuals_ARMA_cpp}},  \code{\link{cll_theta_ARMA_cpp}},
+#'    \code{\link{outputs_STSP_cpp}}, \code{\link{residuals_STSP_cpp}},  \code{\link{cll_theta_STSP_cpp}} and
+#'    \code{\link{solve_de}}, \code{\link{solve_inverse_de}} and \code{\link{ll}}.
+#'
+#' @examples
+#' # generate a random ARMA(2,1) model (3 outputs, 2 inputs)
+#' p = 2
+#' q = 1
+#' m = 2
+#' model = test_armamod(dim = c(m, m), degrees = c(p,q), digits = 2)
+#'
+#' # prepare parameters for "outputs_ARMA_cpp"
+#' A = unclass(model$sys$a)
+#' a0 = A[,,1]
+#' A1 = -A[,,(p+1):2]
+#' dim(A1) = c(m, m*p)
+#' A1 = solve(a0, A1)
+#'
+#' B = unclass(model$sys$b)
+#' dim(B) = c(m, m*(q+1))
+#' B = solve(a0, B)
+#'
+#' # generate random noise sequence (sample size N = 10)
+#' n.obs = 10
+#' u = matrix(rnorm(n.obs*m), nrow = m, ncol = n.obs)
+#'
+#' # generate matrix for the outputs
+#' y = matrix(0, nrow = m, ncol = n.obs)
+#'
+#' # compute outputs
+#' t0 = 2   # start iterations from t>=t0=2
+#' outputs_ARMA_cpp(A1, B, t0, u, y)
+#'
+#' # recompute the disturbances/residuals from the given outputs:
+#' B = unclass(model$sys$b)
+#' ib0 = B[,,1]
+#' B1 = -B[,,(q+1):2]
+#' dim(B1) = c(m, m*q)
+#' B1 = solve(ib0, B1)
+#'
+#' A = unclass(model$sys$a)
+#' dim(A) = c(m, m*(p+1))
+#' A = solve(ib0, A)
+#'
+#' ib0 = solve(ib0)
+#'
+#' uu = u + 0 # "deep copy" of the disturbances
+#' uu[, t0:(n.obs)] = 0 # clear values for t >= t0
+#' residuals_ARMA_cpp(ib0, B1, A, t0 = 2, y, uu, diag(0))
+#' all.equal(u, uu) # check
+#'
+#' # compute directional derivatives of residuals
+#' dU = matrix(0, nrow = n.obs*m, ncol = (m^2)*(p+q+2))
+#' residuals_ARMA_cpp(ib0, B1, A, t0 = 2, y, uu, dU)
+NULL
+
+#' Compute the (concentrated) conditional log likelihood for a statespace system
+#' described by a model template.
+#'
+#' This is an internal helper function, used by the function factory \code{\link{ll_FUN}}. For a more detailed
+#' documentation of the conditional log Likelihood, see \code{\link{ll}}.
+#' The conditional likelihood is computed for the initial state \eqn{a_1} given in the first column `a[,1]` of the
+#' matrix `a`.
+#'
+#' @param th \eqn{(K)} dimensional vector of "deep" parameters.
+#' @param y \eqn{(m,N)} matrix with the observed outputs:
+#'        \eqn{(y_1,y_2,\ldots,y_N)}{(y[1],y[2],...,y[N])}.
+#' @param skip (integer), skip the first residuals, when computing the sample covariance of the
+#'        residuals.
+#' @param concentrated (bool), if TRUE then the *concentrated*, conditional log Likelihood is computed
+#' @param pi \eqn{(m+s,m+s)} matrix, is overwritten with the system matrix
+#'        \eqn{[A,B | C,D]}.
+#' @param H_pi \eqn{(m+s)^2, K)} matrix.
+#' @param h_pi \eqn{((m+s)^2)}-dimensional vector. Note that \code{vec(pi) = H_pi*th + h_pi}.
+#' @param L \eqn{(m,m)} matrix. If (concentrated==FALSE) then L is overwritten with
+#'        the left square of the noise covariance matrix L corresponding
+#'        to the deep parameters th. However, if (concentrated==TRUE) then
+#'        L is overwritten with sample covariance matrix of the computed residuals!
+#' @param H_L \eqn{(m^2, K)} matrix.
+#' @param h_L \eqn{(m^2)}-dimensional vector. Note that
+#'        \code{vec(L) = H_L*th + h_L}.
+#' @param a \eqn{(s,N+1)} matrix. This matrix is overwritten with the (computed) states:
+#'            \eqn{(a_1,a_2,\ldots,a_N,a_{N+1})}{(a[1],a[2],\ldots,a[N],a[N+1])}.
+#'            On input \code{a[,1]} must hold the initial state \eqn{a_1}{a[1]}.
+#' @param u \eqn{(m,N)} matrix. This matrix is overwritten with (computed) residuals:
+#'           \eqn{(u_1,u_2,\ldots,u_N)}{(u[1],u[2],...,u[N])}.
+#' @param dU \eqn{(mN,K)} matrix or \eqn{(0,0)} matrix. This matrix is overwritten with the
+#'           directional derivatives of the residuals. However, if
+#'           the matrix is empty then no derivatives are computed.
+#'
+#' @return (double) log Likelihood
+#'
+#' @export
+#'
+#' @seealso \code{\link{outputs_ARMA_cpp}}, \code{\link{residuals_ARMA_cpp}},  \code{\link{cll_theta_ARMA_cpp}},
+#'    \code{\link{outputs_STSP_cpp}}, \code{\link{residuals_STSP_cpp}},  \code{\link{cll_theta_STSP_cpp}} and
+#'    \code{\link{solve_de}}, \code{\link{solve_inverse_de}} and \code{\link{ll}}.
+#'
+#' @rdname cll_theta_STSP_cpp
+#' @name cll_theta_STSP_cpp
+NULL
+
+#' Outputs of an ARMA systems
+#'
+#' This internal helper function computes the outputs of an ARMA system
+#' \deqn{a_0 y_t + a_1 y_{t-1} + \cdots + a_p y_{t-p} = b_0 u_t + \cdots + b_q u_{t-q}}{
+#'       a[0] y[t] + a[1] y[t-1] + ... + a[p] y[t-p] = b[0] u[t] + ... + b[q] u[t-q]}
+#'
+#' Values \eqn{y_t}{y[t]}, \eqn{u_t}{u[t]} for \eqn{t\leq 0}{t\le 0} are implicitly set to be zero.
+#' However, by starting the iteration with some \eqn{t_0>1}{t0>1} we can enforce non-zero
+#' initial values.
+#'
+#' @note
+#' Use this procedure with care!
+#'
+#' * The procedure does \bold{not} check the input arguments. We require \eqn{m > 0},
+#'   \eqn{p \geq 0}{p \ge 0}, \eqn{n(q+1) \geq 0}{n(q+1)\ge 0} and
+#'   \eqn{1 \leq t_0 \leq N}{1 \le t_0  \le N}.
+#' * The procedure \bold{overwrites} the input argument \code{y}.
+#' * The data matrices are organized columnwise (to avoid memory shuffling)!
+#' * Note also the non standard representation of the coefficient matrices.
+#'
+#' @param A1 \eqn{(m, mp)} matrix \eqn{-a_0^{-1}(a_p,...,a_1)}{-a[0]^{-1}(a[p],...,a[1])}.
+#' @param B \eqn{(m, n(q+1))} matrix \eqn{a_0^{-1}(b_0,...,b_q}{a[0]^{-1}(b[0],...,b[q])}.
+#' @param u \eqn{(n, N)} matrix with the inputs \eqn{(u_1,...,u_N}{(u[1],...,u[N])}.
+#' @param y \eqn{(m, N)} matrix with the outputs \eqn{(y_1,...,y_N}{(y[1],...,y[N])}.
+#' @param t0 integer, start iteration at t = t0.
+#'
+#' @return This RcppArmadillo routine returns \code{NULL} but \bold{overwrites}
+#'         the input argument \code{y} with the computed outputs!
+#'
+#' @export
+#'
+#' @rdname outputs_ARMA_cpp
+#' @name outputs_ARMA_cpp
+#'
+#' @seealso \code{\link{outputs_ARMA_cpp}}, \code{\link{residuals_ARMA_cpp}},  \code{\link{cll_theta_ARMA_cpp}},
+#'    \code{\link{outputs_STSP_cpp}}, \code{\link{residuals_STSP_cpp}},  \code{\link{cll_theta_STSP_cpp}} and
+#'    \code{\link{solve_de}}, \code{\link{solve_inverse_de}} and \code{\link{ll}}.
+#'
+#' @examples
+#' # generate a random ARMA(2,1) model (3 outputs, 2 inputs)
+#' p = 2
+#' q = 1
+#' m = 3
+#' n = 2
+#' model = test_armamod(dim = c(m, n), degrees = c(p,q), digits = 2)
+#' A = unclass(model$sys$a)
+#' a0 = A[,,1]
+#' A1 = -A[,,(p+1):2]
+#' dim(A1) = c(m, m*p)
+#' A1 = solve(a0, A1)
+#' B = unclass(model$sys$b)
+#' dim(B) = c(m, n*(q+1))
+#' B = solve(a0, B)
+#'
+#' # generate random noise sequence (sample size N = 10)
+#' n.obs = 10
+#' u = matrix(rnorm(n.obs*n), nrow = n, ncol = n.obs)
+#' print(u)
+#'
+#' # generate matrix for the outputs
+#' y = matrix(0, nrow = m, ncol = n.obs)
+#'
+#' # call outputs_ARMA_cpp()
+#' outputs_ARMA_cpp(A1, B, t0 = 2, u, y) # start with t>=2
+#' print(u)
+#' print(y)  # y is overwritten with the computed outputs
+outputs_ARMA_cpp <- function(A1, B, t0, u, y) {
+    invisible(.Call(`_RLDM_outputs_ARMA_cpp`, A1, B, t0, u, y))
+}
+
+outputs_STSP_cpp <- function(A, B, C, D, u, a, y) {
+    invisible(.Call(`_RLDM_outputs_STSP_cpp`, A, B, C, D, u, a, y))
+}
+
+fbsolve_STSP_cpp <- function(A, B, C, D, u, au, as, y) {
+    invisible(.Call(`_RLDM_fbsolve_STSP_cpp`, A, B, C, D, u, au, as, y))
+}
+
+solve_rmfd_cpp <- function(poly_inv, poly_fwd, data_in, data_out, t0) {
+    invisible(.Call(`_RLDM_solve_rmfd_cpp`, poly_inv, poly_fwd, data_in, data_out, t0))
+}
+
+residuals_ARMA_cpp <- function(ib0, B1, A, t0, y, u, dU) {
+    invisible(.Call(`_RLDM_residuals_ARMA_cpp`, ib0, B1, A, t0, y, u, dU))
+}
+
+#' Residuals of a statespace system
+#'
+#' This internal helper function computes the residuals
+#' (and the directional derivatives of the residuals)
+#' for a statespace system of the form
+#' \deqn{a_{t+1} = A a_t + B u_t, \; y_t = C a_t + D u_t}{
+#'       a[t+1] = A a[t] + B u[t] and y[t] = C a[t] + D u[t]}
+#' The system must be square and non-empty, i.e. \eqn{m=n>0}.
+#'
+#'
+#' @note
+#' Use this procedure with care!
+#'
+#' * The procedure does \bold{not} check the input arguments.
+#' * The procedure \bold{overwrites} the input arguments
+#'    \code{a}, \code{u} and \code{dU}.
+#' * The data matrices are organized columnwise (to avoid memory shuffling)!
+#'
+#' @param A \eqn{(s,s)} matrix.
+#' @param B \eqn{(s,m)} matrix.
+#' @param C \eqn{(m,s)} matrix.
+#' @param D \eqn{(m,m)} matrix, must be regular.
+#' @param y \eqn{(m,N)} matrix with the outputs: \eqn{(y_1,y_2,\ldots,y_N)}{(y[1],y[2],...,y[N])}.
+#' @param a \eqn{(s,N+1)} matrix. This matrix is overwritten with the (computed) states:
+#'            \eqn{(a_1,a_2,\ldots,a_N,a_{N+1})}{(a[1],a[2],\ldots,a[N],a[N+1])}.
+#'            On input \code{a[,1]} must hold the initial state \eqn{a_1}{a[1]}.
+#' @param u \eqn{(m,N)} matrix. This matrix is overwritten with (computed) residuals:
+#'           \eqn{(u_1,u_2,\ldots,u_N)}{(u[1],u[2],...,u[N])}.
+#' @param dPI \eqn{((m+s)^2,K)} matrix.
+#' @param dU \eqn{(mN,K)} matrix or \eqn{(0,0)} matrix. This matrix is overwritten with the
+#'           directional derivatives of the residuals. However, if
+#'           the matrix is empty then no derivatives are computed.
+#'
+#' @return This RcppArmadillo implementation returns \code{NULL} but \bold{overwrites} the input
+#'         arguments \code{a}, \code{u} and \code{dU}!
+#'
+#' @export
+#'
+#' @rdname residuals_STSP_cpp
+#' @name residuals_STSP_cpp
+#'
+#' @seealso \code{\link{outputs_ARMA_cpp}}, \code{\link{residuals_ARMA_cpp}},  \code{\link{cll_theta_ARMA_cpp}},
+#'    \code{\link{outputs_STSP_cpp}}, \code{\link{residuals_STSP_cpp}},  \code{\link{cll_theta_STSP_cpp}} and
+#'    \code{\link{solve_de}}, \code{\link{solve_inverse_de}} and \code{\link{ll}}.
+#'
+#' @examples
+#' # generate a random statespace model (3 outputs, 3 inputs and 4 states)
+#' m = 2
+#' s = 3
+#' model = test_stspmod(dim = c(m, m), s = s, digits = 2)
+#'
+#' # generate random noise sequence (sample size N = 10)
+#' n.obs = 10
+#' u = matrix(rnorm(n.obs*m), nrow = m, ncol = n.obs)
+#'
+#' # generate matrix for the state sequence
+#' a = matrix(0, nrow = s, ncol = n.obs+1)
+#' a[,1] = rnorm(s) # random initial state a[1]
+#'
+#' # generate matrix for the outputs
+#' y = matrix(0, nrow = m, ncol = n.obs)
+#'
+#' # compute outputs and states
+#' outputs_STSP_cpp(model$sys$A, model$sys$B, model$sys$C, model$sys$D, u, a, y)
+#'
+#' # recompute the states and disturbances/residuals from the given outputs:
+#' uu = u + 0 # "deep copy" of the disturbances
+#' aa = a + 0 # and the states
+#' aa[, 2:(n.obs+1)] = 0 # clear all states a[t], t > 1
+#' residuals_STSP_cpp(model$sys$A, model$sys$B, model$sys$C, model$sys$D,
+#'                    y, aa, uu, diag(0), diag(0))
+#' all.equal(u, uu) # check
+#' all.equal(a, aa) # check
+#'
+#' # compute directional derivatives of residuals
+#' dPI = diag((m+s)^2)
+#' dU = matrix(0, nrow = n.obs*m, ncol = ncol(dPI))
+#' residuals_STSP_cpp(model$sys$A, model$sys$B, model$sys$C, model$sys$D,
+#'                    y, aa, uu, dPI, dU)
+#'
+#' # check the directional derivatives
+#' eps = 1e-8
+#' dU_num = matrix(0, nrow = m*n.obs, ncol = (m+s)^2)
+#' dPI = matrix(0, nrow = (m+s), ncol = (m+s))
+#' for (k in (1:((m+s)^2))) {
+#'   dPI[] = 0
+#'   dPI[k] = eps
+#'   uu = u + 0 # "deep copy" of the disturbances
+#'   aa = a + 0 # and the states
+#'   residuals_STSP_cpp(model$sys$A + dPI[1:s,1:s],
+#'                      model$sys$B + dPI[1:s,(s+1):(s+m)],
+#'                      model$sys$C + dPI[(s+1):(s+m),1:s],
+#'                      model$sys$D + dPI[(s+1):(s+m),(s+1):(s+m)],
+#'                      y, aa, uu, diag(0), diag(0))
+#'   dU_num[, k] = c(uu - u )/eps  # num. approx. of the derivative in direction "dPI"
+#' }
+#' # relative error of the numerical approximation
+#' junk = (abs(dU)+abs(dU_num))
+#' junk[junk == 0] = 1
+#' 2*abs(dU_num - dU)/junk
+#'
+residuals_STSP_cpp <- function(A, B, C, D, y, a, u, dPI, dU) {
+    invisible(.Call(`_RLDM_residuals_STSP_cpp`, A, B, C, D, y, a, u, dPI, dU))
+}
+
+#' Compute the (concentrated) conditional log likelihood for ARMA models
+#' described by a model template.
+#'
+#' This internal helper function computes the (concentrated) conditional log Likelihood
+#' of ARMA systems of the form
+#' \deqn{a_0 y_t + a_1 y_{t-1} + \cdots + a_p y_{t-p} = b_0 u_t + \cdots + b_q u_{t-q}}{
+#'       a[0] y[t] + a[1] y[t-1] + ... + a[p] y[t-p] = b[0] u[t] + ... + b[q] u[t-q]}
+#' The conditional likelihood is computed for **zero** initial values \eqn{u_s=y_s=0}{u[s]=y[s]=0} for
+#' \eqn{s\leq 0}{s\le 0}.
+#'
+#' This function is mainly used by the function factory \code{\link{ll_FUN}}. For a more detailed
+#' documentation of the (concentrated) conditional log Likelihood, see \code{\link{ll}}.
+#'
+#' The procedure first constructs the ARMA parameter matrices from the given vector
+#' \code{th} of "deep" parameters.
+#'
+#' \itemize{
+#' \item  AR parameters \code{vec((a[0],a[1],...,a[p])) = h_A + H_A * th}.
+#' \item  MA parameters \code{vec(b[0]) = h_b + H_b * th} and
+#'        \code{vec(-(b[q],...,b[1])) = h_B + H_B * th}
+#' \item  Left square root of noise covariance matrix
+#'        \eqn{\Sigma = LL'} and \code{vec(L) = h_L + H_L * th}.
+#' }
+#'
+#' The residuals (and their directional derivatives) are computed with
+#' \code{\link{residuals_ARMA_cpp}}.
+#'
+#' @note
+#' Use this procedure with care!
+#'
+#' * The procedure does \bold{not} check the input arguments.
+#' * The procedure \bold{overwrites} some of the input arguments
+#' * The data matrices are organized columnwise (to avoid memory shuffling)!
+#' * Note also the non standard representation of the coefficient matrices.
+#'
+#' @param th \eqn{(K)} dimensional vector of "deep" parameters.
+#' @param y \eqn{(m,N)} matrix with the observed outputs:
+#'        \eqn{(y_1,y_2,\ldots,y_N)}{(y[1],y[2],...,y[N])}.
+#' @param skip (integer), omit the first "skip" residuals, when computing the likelihood.
+#' @param concentrated (bool), if TRUE then the *concentrated*, conditional log Likelihood is computed
+#' @param ib0 \eqn{(m, m)} matrix, is **overwritten** with the matrix \eqn{b_0^{-1}a_0}{b[0]^(-1)a[0]}.
+#' @param H_b \eqn{(m^2, K)} matrix.
+#' @param h_b \eqn{((m^2)}-dimensional vector. Note that
+#'         \code{vec(b[0]) = H_b*th + h_b}.
+#' @param B1 \eqn{(m, mq)} matrix, is **overwritten** with
+#'        \eqn{-b_0^{-1}(b_q,...,b_1)}{-b[0]^(-1)(b[q],...,b[1])}.
+#' @param H_B \eqn{((m^2)*q, K)} matrix.
+#' @param h_B \eqn{((m^2)*q)}-dimensional vector. Note that
+#'        \code{vec(-(b[q],...,b[1])) = H_B*th + h_B}.
+#' @param a0 \eqn{(m, m)} matrix, is **overwritten** with
+#'        \eqn{a_0}{a[0]}.
+#' @param A \eqn{(m, m(q+1))} matrix, is **overwritten** with
+#'        \eqn{b_0^{-1}(a_0,...,a_p}{b[0]^(-1)(a[0],...,a[p])}.
+#' @param H_A \eqn{((m^2)*(p+1), K)} matrix.
+#' @param h_A \eqn{((m^2)*(p+1))}-dimensional vector. Note that
+#'        \code{vec((a[0],a[1],...,a[p])) = H_A*th + h_A}.
+#' @param L \eqn{(m,m)} matrix. If (concentrated==FALSE) then \code{L} is **overwritten** with
+#'        the left square \eqn{L} of the noise covariance matrix \eqn{\Sigma=LL'} corresponding
+#'        to the deep parameters th. However, if (concentrated==TRUE) then
+#'        L is **overwritten** with sample covariance matrix of the computed residuals!
+#' @param H_L \eqn{(m^2, K)} matrix.
+#' @param h_L \eqn{(m^2)}-dimensional vector. Note that
+#'        \code{vec(L) = H_L*th + h_L}.
+#' @param u \eqn{(m,N)} matrix. This matrix is **overwritten** with (computed) residuals:
+#'           \eqn{(u_1,u_2,\ldots,u_N)}{(u[1],u[2],...,u[N])}.
+#' @param dU \eqn{(mN,(m^2)(p+q+2))} matrix or \eqn{(0,0)} matrix. If non empty this
+#'           matrix is **overwritten** with the
+#'           directional derivatives of the residuals. However, if
+#'           the matrix is empty then no derivatives are computed.
+#'
+#' @seealso \code{\link{outputs_ARMA_cpp}}, \code{\link{residuals_ARMA_cpp}},  \code{\link{cll_theta_ARMA_cpp}},
+#'    \code{\link{outputs_STSP_cpp}}, \code{\link{residuals_STSP_cpp}},  \code{\link{cll_theta_STSP_cpp}} and
+#'    \code{\link{solve_de}}, \code{\link{solve_inverse_de}} and \code{\link{ll}}.
+#'
+#' @return (double) log Likelihood
+#'
+#' @export
+#'
+#' @rdname cll_theta_ARMA_cpp
+#' @name cll_theta_ARMA_cpp
+cll_theta_ARMA_cpp <- function(th, y, skip, concentrated, ib0, H_b, h_b, B1, H_B, h_B, a0, A, H_A, h_A, L, H_L, h_L, u, dU) {
+    .Call(`_RLDM_cll_theta_ARMA_cpp`, th, y, skip, concentrated, ib0, H_b, h_b, B1, H_B, h_B, a0, A, H_A, h_A, L, H_L, h_L, u, dU)
+}
+
+cll_theta_STSP_cpp <- function(th, y, skip, concentrated, pi, H_pi, h_pi, L, H_L, h_L, a, u, dU) {
+    .Call(`_RLDM_cll_theta_STSP_cpp`, th, y, skip, concentrated, pi, H_pi, h_pi, L, H_L, h_L, a, u, dU)
+}
+
+# Register entry points for exported C++ functions
+methods::setLoadAction(function(ns) {
+    .Call(`_RLDM_RcppExport_registerCCallable`)
+})
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
new file mode 100644
index 0000000..646472b
--- /dev/null
+++ b/src/RcppExports.cpp
@@ -0,0 +1,617 @@
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include "../inst/include/RLDM.h"
+#include <RcppArmadillo.h>
+#include <Rcpp.h>
+#include <string>
+#include <set>
+
+using namespace Rcpp;
+
+#ifdef RCPP_USE_GLOBAL_ROSTREAM
+Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
+Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
+#endif
+
+// kf_cpp
+Rcpp::List kf_cpp(const arma::mat& A, const arma::mat& C, const arma::mat& Q, const arma::mat& R, const arma::mat& S, const arma::mat& y_t, const arma::mat& P1, const arma::colvec& a1);
+static SEXP _RLDM_kf_cpp_try(SEXP ASEXP, SEXP CSEXP, SEXP QSEXP, SEXP RSEXP, SEXP SSEXP, SEXP y_tSEXP, SEXP P1SEXP, SEXP a1SEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type C(CSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type Q(QSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type R(RSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type S(SSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type y_t(y_tSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type P1(P1SEXP);
+    Rcpp::traits::input_parameter< const arma::colvec& >::type a1(a1SEXP);
+    rcpp_result_gen = Rcpp::wrap(kf_cpp(A, C, Q, R, S, y_t, P1, a1));
+    return rcpp_result_gen;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_kf_cpp(SEXP ASEXP, SEXP CSEXP, SEXP QSEXP, SEXP RSEXP, SEXP SSEXP, SEXP y_tSEXP, SEXP P1SEXP, SEXP a1SEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_kf_cpp_try(ASEXP, CSEXP, QSEXP, RSEXP, SSEXP, y_tSEXP, P1SEXP, a1SEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+// kf2_cpp
+Rcpp::List kf2_cpp(const arma::mat& A, const arma::mat& C, const arma::mat& H_t, const arma::mat& y_t, arma::mat& P1_R, const arma::colvec& a1);
+static SEXP _RLDM_kf2_cpp_try(SEXP ASEXP, SEXP CSEXP, SEXP H_tSEXP, SEXP y_tSEXP, SEXP P1_RSEXP, SEXP a1SEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type C(CSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type H_t(H_tSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type y_t(y_tSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type P1_R(P1_RSEXP);
+    Rcpp::traits::input_parameter< const arma::colvec& >::type a1(a1SEXP);
+    rcpp_result_gen = Rcpp::wrap(kf2_cpp(A, C, H_t, y_t, P1_R, a1));
+    return rcpp_result_gen;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_kf2_cpp(SEXP ASEXP, SEXP CSEXP, SEXP H_tSEXP, SEXP y_tSEXP, SEXP P1_RSEXP, SEXP a1SEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_kf2_cpp_try(ASEXP, CSEXP, H_tSEXP, y_tSEXP, P1_RSEXP, a1SEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+// ll_kf_cpp
+double ll_kf_cpp(const arma::mat& A, const arma::mat& C, const arma::mat& Q, const arma::mat& R, const arma::mat& S, const arma::mat& y_t, const arma::mat& P1, const arma::colvec& a1, double tol);
+static SEXP _RLDM_ll_kf_cpp_try(SEXP ASEXP, SEXP CSEXP, SEXP QSEXP, SEXP RSEXP, SEXP SSEXP, SEXP y_tSEXP, SEXP P1SEXP, SEXP a1SEXP, SEXP tolSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type C(CSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type Q(QSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type R(RSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type S(SSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type y_t(y_tSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type P1(P1SEXP);
+    Rcpp::traits::input_parameter< const arma::colvec& >::type a1(a1SEXP);
+    Rcpp::traits::input_parameter< double >::type tol(tolSEXP);
+    rcpp_result_gen = Rcpp::wrap(ll_kf_cpp(A, C, Q, R, S, y_t, P1, a1, tol));
+    return rcpp_result_gen;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_ll_kf_cpp(SEXP ASEXP, SEXP CSEXP, SEXP QSEXP, SEXP RSEXP, SEXP SSEXP, SEXP y_tSEXP, SEXP P1SEXP, SEXP a1SEXP, SEXP tolSEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_ll_kf_cpp_try(ASEXP, CSEXP, QSEXP, RSEXP, SSEXP, y_tSEXP, P1SEXP, a1SEXP, tolSEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+// ll_kf2_cpp
+double ll_kf2_cpp(arma::mat& A, arma::mat& C, arma::mat& H_t, arma::mat& y_t, arma::mat& P1_R, arma::colvec& a1, double tol);
+static SEXP _RLDM_ll_kf2_cpp_try(SEXP ASEXP, SEXP CSEXP, SEXP H_tSEXP, SEXP y_tSEXP, SEXP P1_RSEXP, SEXP a1SEXP, SEXP tolSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::traits::input_parameter< arma::mat& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type C(CSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type H_t(H_tSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type y_t(y_tSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type P1_R(P1_RSEXP);
+    Rcpp::traits::input_parameter< arma::colvec& >::type a1(a1SEXP);
+    Rcpp::traits::input_parameter< double >::type tol(tolSEXP);
+    rcpp_result_gen = Rcpp::wrap(ll_kf2_cpp(A, C, H_t, y_t, P1_R, a1, tol));
+    return rcpp_result_gen;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_ll_kf2_cpp(SEXP ASEXP, SEXP CSEXP, SEXP H_tSEXP, SEXP y_tSEXP, SEXP P1_RSEXP, SEXP a1SEXP, SEXP tolSEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_ll_kf2_cpp_try(ASEXP, CSEXP, H_tSEXP, y_tSEXP, P1_RSEXP, a1SEXP, tolSEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+// ll_kf_theta_cpp
+double ll_kf_theta_cpp(const arma::vec& theta, const arma::mat& y, arma::mat& SYS, const arma::mat& H_SYS, const arma::vec& h_SYS, arma::mat& sigma_L, const arma::mat& H_sigma_L, const arma::vec& h_sigma_L, arma::mat& VAR, arma::mat& P1, double tol, double err);
+static SEXP _RLDM_ll_kf_theta_cpp_try(SEXP thetaSEXP, SEXP ySEXP, SEXP SYSSEXP, SEXP H_SYSSEXP, SEXP h_SYSSEXP, SEXP sigma_LSEXP, SEXP H_sigma_LSEXP, SEXP h_sigma_LSEXP, SEXP VARSEXP, SEXP P1SEXP, SEXP tolSEXP, SEXP errSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::traits::input_parameter< const arma::vec& >::type theta(thetaSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type SYS(SYSSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type H_SYS(H_SYSSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type h_SYS(h_SYSSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type sigma_L(sigma_LSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type H_sigma_L(H_sigma_LSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type h_sigma_L(h_sigma_LSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type VAR(VARSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type P1(P1SEXP);
+    Rcpp::traits::input_parameter< double >::type tol(tolSEXP);
+    Rcpp::traits::input_parameter< double >::type err(errSEXP);
+    rcpp_result_gen = Rcpp::wrap(ll_kf_theta_cpp(theta, y, SYS, H_SYS, h_SYS, sigma_L, H_sigma_L, h_sigma_L, VAR, P1, tol, err));
+    return rcpp_result_gen;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_ll_kf_theta_cpp(SEXP thetaSEXP, SEXP ySEXP, SEXP SYSSEXP, SEXP H_SYSSEXP, SEXP h_SYSSEXP, SEXP sigma_LSEXP, SEXP H_sigma_LSEXP, SEXP h_sigma_LSEXP, SEXP VARSEXP, SEXP P1SEXP, SEXP tolSEXP, SEXP errSEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_ll_kf_theta_cpp_try(thetaSEXP, ySEXP, SYSSEXP, H_SYSSEXP, h_SYSSEXP, sigma_LSEXP, H_sigma_LSEXP, h_sigma_LSEXP, VARSEXP, P1SEXP, tolSEXP, errSEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+// outputs_ARMA_cpp
+void outputs_ARMA_cpp(const arma::mat& A1, const arma::mat& B, int t0, const arma::mat& u, arma::mat& y);
+static SEXP _RLDM_outputs_ARMA_cpp_try(SEXP A1SEXP, SEXP BSEXP, SEXP t0SEXP, SEXP uSEXP, SEXP ySEXP) {
+BEGIN_RCPP
+    Rcpp::traits::input_parameter< const arma::mat& >::type A1(A1SEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type B(BSEXP);
+    Rcpp::traits::input_parameter< int >::type t0(t0SEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type u(uSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type y(ySEXP);
+    outputs_ARMA_cpp(A1, B, t0, u, y);
+    return R_NilValue;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_outputs_ARMA_cpp(SEXP A1SEXP, SEXP BSEXP, SEXP t0SEXP, SEXP uSEXP, SEXP ySEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_outputs_ARMA_cpp_try(A1SEXP, BSEXP, t0SEXP, uSEXP, ySEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+// outputs_STSP_cpp
+void outputs_STSP_cpp(const arma::mat& A, const arma::mat& B, const arma::mat& C, const arma::mat& D, const arma::mat& u, arma::mat& a, arma::mat& y);
+static SEXP _RLDM_outputs_STSP_cpp_try(SEXP ASEXP, SEXP BSEXP, SEXP CSEXP, SEXP DSEXP, SEXP uSEXP, SEXP aSEXP, SEXP ySEXP) {
+BEGIN_RCPP
+    Rcpp::traits::input_parameter< const arma::mat& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type B(BSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type C(CSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type D(DSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type u(uSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type a(aSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type y(ySEXP);
+    outputs_STSP_cpp(A, B, C, D, u, a, y);
+    return R_NilValue;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_outputs_STSP_cpp(SEXP ASEXP, SEXP BSEXP, SEXP CSEXP, SEXP DSEXP, SEXP uSEXP, SEXP aSEXP, SEXP ySEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_outputs_STSP_cpp_try(ASEXP, BSEXP, CSEXP, DSEXP, uSEXP, aSEXP, ySEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+// fbsolve_STSP_cpp
+void fbsolve_STSP_cpp(const arma::mat& A, const arma::mat& B, const arma::mat& C, const arma::mat& D, const arma::mat& u, arma::mat& au, arma::mat& as, arma::mat& y);
+static SEXP _RLDM_fbsolve_STSP_cpp_try(SEXP ASEXP, SEXP BSEXP, SEXP CSEXP, SEXP DSEXP, SEXP uSEXP, SEXP auSEXP, SEXP asSEXP, SEXP ySEXP) {
+BEGIN_RCPP
+    Rcpp::traits::input_parameter< const arma::mat& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type B(BSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type C(CSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type D(DSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type u(uSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type au(auSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type as(asSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type y(ySEXP);
+    fbsolve_STSP_cpp(A, B, C, D, u, au, as, y);
+    return R_NilValue;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_fbsolve_STSP_cpp(SEXP ASEXP, SEXP BSEXP, SEXP CSEXP, SEXP DSEXP, SEXP uSEXP, SEXP auSEXP, SEXP asSEXP, SEXP ySEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_fbsolve_STSP_cpp_try(ASEXP, BSEXP, CSEXP, DSEXP, uSEXP, auSEXP, asSEXP, ySEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+// solve_rmfd_cpp
+void solve_rmfd_cpp(const arma::mat& poly_inv, const arma::mat& poly_fwd, arma::mat& data_in, arma::mat& data_out, int t0);
+static SEXP _RLDM_solve_rmfd_cpp_try(SEXP poly_invSEXP, SEXP poly_fwdSEXP, SEXP data_inSEXP, SEXP data_outSEXP, SEXP t0SEXP) {
+BEGIN_RCPP
+    Rcpp::traits::input_parameter< const arma::mat& >::type poly_inv(poly_invSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type poly_fwd(poly_fwdSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type data_in(data_inSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type data_out(data_outSEXP);
+    Rcpp::traits::input_parameter< int >::type t0(t0SEXP);
+    solve_rmfd_cpp(poly_inv, poly_fwd, data_in, data_out, t0);
+    return R_NilValue;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_solve_rmfd_cpp(SEXP poly_invSEXP, SEXP poly_fwdSEXP, SEXP data_inSEXP, SEXP data_outSEXP, SEXP t0SEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_solve_rmfd_cpp_try(poly_invSEXP, poly_fwdSEXP, data_inSEXP, data_outSEXP, t0SEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+// residuals_ARMA_cpp
+void residuals_ARMA_cpp(const arma::mat& ib0, const arma::mat& B1, const arma::mat& A, int t0, const arma::mat& y, arma::mat& u, arma::mat& dU);
+static SEXP _RLDM_residuals_ARMA_cpp_try(SEXP ib0SEXP, SEXP B1SEXP, SEXP ASEXP, SEXP t0SEXP, SEXP ySEXP, SEXP uSEXP, SEXP dUSEXP) {
+BEGIN_RCPP
+    Rcpp::traits::input_parameter< const arma::mat& >::type ib0(ib0SEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type B1(B1SEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< int >::type t0(t0SEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type u(uSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type dU(dUSEXP);
+    residuals_ARMA_cpp(ib0, B1, A, t0, y, u, dU);
+    return R_NilValue;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_residuals_ARMA_cpp(SEXP ib0SEXP, SEXP B1SEXP, SEXP ASEXP, SEXP t0SEXP, SEXP ySEXP, SEXP uSEXP, SEXP dUSEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_residuals_ARMA_cpp_try(ib0SEXP, B1SEXP, ASEXP, t0SEXP, ySEXP, uSEXP, dUSEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+// residuals_STSP_cpp
+void residuals_STSP_cpp(const arma::mat& A, const arma::mat& B, const arma::mat& C, const arma::mat& D, const arma::mat& y, arma::mat& a, arma::mat& u, const arma::mat& dPI, arma::mat& dU);
+static SEXP _RLDM_residuals_STSP_cpp_try(SEXP ASEXP, SEXP BSEXP, SEXP CSEXP, SEXP DSEXP, SEXP ySEXP, SEXP aSEXP, SEXP uSEXP, SEXP dPISEXP, SEXP dUSEXP) {
+BEGIN_RCPP
+    Rcpp::traits::input_parameter< const arma::mat& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type B(BSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type C(CSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type D(DSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type a(aSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type u(uSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type dPI(dPISEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type dU(dUSEXP);
+    residuals_STSP_cpp(A, B, C, D, y, a, u, dPI, dU);
+    return R_NilValue;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_residuals_STSP_cpp(SEXP ASEXP, SEXP BSEXP, SEXP CSEXP, SEXP DSEXP, SEXP ySEXP, SEXP aSEXP, SEXP uSEXP, SEXP dPISEXP, SEXP dUSEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_residuals_STSP_cpp_try(ASEXP, BSEXP, CSEXP, DSEXP, ySEXP, aSEXP, uSEXP, dPISEXP, dUSEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+// cll_theta_ARMA_cpp
+double cll_theta_ARMA_cpp(const arma::vec& th, const arma::mat& y, unsigned long int skip, bool concentrated, arma::mat& ib0, const arma::mat& H_b, const arma::vec& h_b, arma::mat& B1, const arma::mat& H_B, const arma::vec& h_B, arma::mat& a0, arma::mat& A, const arma::mat& H_A, const arma::vec& h_A, arma::mat& L, const arma::mat& H_L, const arma::vec& h_L, arma::mat& u, arma::mat& dU);
+static SEXP _RLDM_cll_theta_ARMA_cpp_try(SEXP thSEXP, SEXP ySEXP, SEXP skipSEXP, SEXP concentratedSEXP, SEXP ib0SEXP, SEXP H_bSEXP, SEXP h_bSEXP, SEXP B1SEXP, SEXP H_BSEXP, SEXP h_BSEXP, SEXP a0SEXP, SEXP ASEXP, SEXP H_ASEXP, SEXP h_ASEXP, SEXP LSEXP, SEXP H_LSEXP, SEXP h_LSEXP, SEXP uSEXP, SEXP dUSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::traits::input_parameter< const arma::vec& >::type th(thSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< unsigned long int >::type skip(skipSEXP);
+    Rcpp::traits::input_parameter< bool >::type concentrated(concentratedSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type ib0(ib0SEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type H_b(H_bSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type h_b(h_bSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type B1(B1SEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type H_B(H_BSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type h_B(h_BSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type a0(a0SEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type H_A(H_ASEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type h_A(h_ASEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type L(LSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type H_L(H_LSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type h_L(h_LSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type u(uSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type dU(dUSEXP);
+    rcpp_result_gen = Rcpp::wrap(cll_theta_ARMA_cpp(th, y, skip, concentrated, ib0, H_b, h_b, B1, H_B, h_B, a0, A, H_A, h_A, L, H_L, h_L, u, dU));
+    return rcpp_result_gen;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_cll_theta_ARMA_cpp(SEXP thSEXP, SEXP ySEXP, SEXP skipSEXP, SEXP concentratedSEXP, SEXP ib0SEXP, SEXP H_bSEXP, SEXP h_bSEXP, SEXP B1SEXP, SEXP H_BSEXP, SEXP h_BSEXP, SEXP a0SEXP, SEXP ASEXP, SEXP H_ASEXP, SEXP h_ASEXP, SEXP LSEXP, SEXP H_LSEXP, SEXP h_LSEXP, SEXP uSEXP, SEXP dUSEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_cll_theta_ARMA_cpp_try(thSEXP, ySEXP, skipSEXP, concentratedSEXP, ib0SEXP, H_bSEXP, h_bSEXP, B1SEXP, H_BSEXP, h_BSEXP, a0SEXP, ASEXP, H_ASEXP, h_ASEXP, LSEXP, H_LSEXP, h_LSEXP, uSEXP, dUSEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+// cll_theta_STSP_cpp
+double cll_theta_STSP_cpp(const arma::vec& th, const arma::mat& y, unsigned long int skip, bool concentrated, arma::mat& pi, const arma::mat& H_pi, const arma::vec& h_pi, arma::mat& L, const arma::mat& H_L, const arma::vec& h_L, arma::mat& a, arma::mat& u, arma::mat& dU);
+static SEXP _RLDM_cll_theta_STSP_cpp_try(SEXP thSEXP, SEXP ySEXP, SEXP skipSEXP, SEXP concentratedSEXP, SEXP piSEXP, SEXP H_piSEXP, SEXP h_piSEXP, SEXP LSEXP, SEXP H_LSEXP, SEXP h_LSEXP, SEXP aSEXP, SEXP uSEXP, SEXP dUSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::traits::input_parameter< const arma::vec& >::type th(thSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< unsigned long int >::type skip(skipSEXP);
+    Rcpp::traits::input_parameter< bool >::type concentrated(concentratedSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type pi(piSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type H_pi(H_piSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type h_pi(h_piSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type L(LSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type H_L(H_LSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type h_L(h_LSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type a(aSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type u(uSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type dU(dUSEXP);
+    rcpp_result_gen = Rcpp::wrap(cll_theta_STSP_cpp(th, y, skip, concentrated, pi, H_pi, h_pi, L, H_L, h_L, a, u, dU));
+    return rcpp_result_gen;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _RLDM_cll_theta_STSP_cpp(SEXP thSEXP, SEXP ySEXP, SEXP skipSEXP, SEXP concentratedSEXP, SEXP piSEXP, SEXP H_piSEXP, SEXP h_piSEXP, SEXP LSEXP, SEXP H_LSEXP, SEXP h_LSEXP, SEXP aSEXP, SEXP uSEXP, SEXP dUSEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_RLDM_cll_theta_STSP_cpp_try(thSEXP, ySEXP, skipSEXP, concentratedSEXP, piSEXP, H_piSEXP, h_piSEXP, LSEXP, H_LSEXP, h_LSEXP, aSEXP, uSEXP, dUSEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
+
+// validate (ensure exported C++ functions exist before calling them)
+static int _RLDM_RcppExport_validate(const char* sig) { 
+    static std::set<std::string> signatures;
+    if (signatures.empty()) {
+        signatures.insert("Rcpp::List(*kf_cpp)(const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,const arma::colvec&)");
+        signatures.insert("Rcpp::List(*kf2_cpp)(const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,arma::mat&,const arma::colvec&)");
+        signatures.insert("double(*ll_kf_cpp)(const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,const arma::colvec&,double)");
+        signatures.insert("double(*ll_kf2_cpp)(arma::mat&,arma::mat&,arma::mat&,arma::mat&,arma::mat&,arma::colvec&,double)");
+        signatures.insert("double(*ll_kf_theta_cpp)(const arma::vec&,const arma::mat&,arma::mat&,const arma::mat&,const arma::vec&,arma::mat&,const arma::mat&,const arma::vec&,arma::mat&,arma::mat&,double,double)");
+        signatures.insert("void(*outputs_ARMA_cpp)(const arma::mat&,const arma::mat&,int,const arma::mat&,arma::mat&)");
+        signatures.insert("void(*outputs_STSP_cpp)(const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,arma::mat&,arma::mat&)");
+        signatures.insert("void(*fbsolve_STSP_cpp)(const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,arma::mat&,arma::mat&,arma::mat&)");
+        signatures.insert("void(*solve_rmfd_cpp)(const arma::mat&,const arma::mat&,arma::mat&,arma::mat&,int)");
+        signatures.insert("void(*residuals_ARMA_cpp)(const arma::mat&,const arma::mat&,const arma::mat&,int,const arma::mat&,arma::mat&,arma::mat&)");
+        signatures.insert("void(*residuals_STSP_cpp)(const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,const arma::mat&,arma::mat&,arma::mat&,const arma::mat&,arma::mat&)");
+        signatures.insert("double(*cll_theta_ARMA_cpp)(const arma::vec&,const arma::mat&,unsigned long int,bool,arma::mat&,const arma::mat&,const arma::vec&,arma::mat&,const arma::mat&,const arma::vec&,arma::mat&,arma::mat&,const arma::mat&,const arma::vec&,arma::mat&,const arma::mat&,const arma::vec&,arma::mat&,arma::mat&)");
+        signatures.insert("double(*cll_theta_STSP_cpp)(const arma::vec&,const arma::mat&,unsigned long int,bool,arma::mat&,const arma::mat&,const arma::vec&,arma::mat&,const arma::mat&,const arma::vec&,arma::mat&,arma::mat&,arma::mat&)");
+    }
+    return signatures.find(sig) != signatures.end();
+}
+
+// registerCCallable (register entry points for exported C++ functions)
+RcppExport SEXP _RLDM_RcppExport_registerCCallable() { 
+    R_RegisterCCallable("RLDM", "_RLDM_kf_cpp", (DL_FUNC)_RLDM_kf_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_kf2_cpp", (DL_FUNC)_RLDM_kf2_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_ll_kf_cpp", (DL_FUNC)_RLDM_ll_kf_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_ll_kf2_cpp", (DL_FUNC)_RLDM_ll_kf2_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_ll_kf_theta_cpp", (DL_FUNC)_RLDM_ll_kf_theta_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_outputs_ARMA_cpp", (DL_FUNC)_RLDM_outputs_ARMA_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_outputs_STSP_cpp", (DL_FUNC)_RLDM_outputs_STSP_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_fbsolve_STSP_cpp", (DL_FUNC)_RLDM_fbsolve_STSP_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_solve_rmfd_cpp", (DL_FUNC)_RLDM_solve_rmfd_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_residuals_ARMA_cpp", (DL_FUNC)_RLDM_residuals_ARMA_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_residuals_STSP_cpp", (DL_FUNC)_RLDM_residuals_STSP_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_cll_theta_ARMA_cpp", (DL_FUNC)_RLDM_cll_theta_ARMA_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_cll_theta_STSP_cpp", (DL_FUNC)_RLDM_cll_theta_STSP_cpp_try);
+    R_RegisterCCallable("RLDM", "_RLDM_RcppExport_validate", (DL_FUNC)_RLDM_RcppExport_validate);
+    return R_NilValue;
+}
+
+static const R_CallMethodDef CallEntries[] = {
+    {"_RLDM_kf_cpp", (DL_FUNC) &_RLDM_kf_cpp, 8},
+    {"_RLDM_kf2_cpp", (DL_FUNC) &_RLDM_kf2_cpp, 6},
+    {"_RLDM_ll_kf_cpp", (DL_FUNC) &_RLDM_ll_kf_cpp, 9},
+    {"_RLDM_ll_kf2_cpp", (DL_FUNC) &_RLDM_ll_kf2_cpp, 7},
+    {"_RLDM_ll_kf_theta_cpp", (DL_FUNC) &_RLDM_ll_kf_theta_cpp, 12},
+    {"_RLDM_outputs_ARMA_cpp", (DL_FUNC) &_RLDM_outputs_ARMA_cpp, 5},
+    {"_RLDM_outputs_STSP_cpp", (DL_FUNC) &_RLDM_outputs_STSP_cpp, 7},
+    {"_RLDM_fbsolve_STSP_cpp", (DL_FUNC) &_RLDM_fbsolve_STSP_cpp, 8},
+    {"_RLDM_solve_rmfd_cpp", (DL_FUNC) &_RLDM_solve_rmfd_cpp, 5},
+    {"_RLDM_residuals_ARMA_cpp", (DL_FUNC) &_RLDM_residuals_ARMA_cpp, 7},
+    {"_RLDM_residuals_STSP_cpp", (DL_FUNC) &_RLDM_residuals_STSP_cpp, 9},
+    {"_RLDM_cll_theta_ARMA_cpp", (DL_FUNC) &_RLDM_cll_theta_ARMA_cpp, 19},
+    {"_RLDM_cll_theta_STSP_cpp", (DL_FUNC) &_RLDM_cll_theta_STSP_cpp, 13},
+    {"_RLDM_RcppExport_registerCCallable", (DL_FUNC) &_RLDM_RcppExport_registerCCallable, 0},
+    {NULL, NULL, 0}
+};
+
+RcppExport void R_init_RLDM(DllInfo *dll) {
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+    R_useDynamicSymbols(dll, FALSE);
+}