diff --git a/DESCRIPTION b/DESCRIPTION
index 2f8dbbf..2c23f4f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -4,15 +4,15 @@ Title: Trend-Cycle Extraction with Linear Filters based on JDemetra+ v3.x
 Version: 2.1.1.9000
 Authors@R: c(
   person("Jean", "Palate", role = c("aut"),
-                email = "palatejean@gmail.com"),
-  person("Alain", "Quartier-la-Tente", role = c("aut", "cre"),
-                 email = "alain.quartier@yahoo.fr",
-                 comment = c(ORCID = "0000-0001-7890-3857")),
-  person("Tanguy", "Barthelemy", role = c("ctb", "art"),
-                 email ="tanguy.barthelemy@insee.fr"),
+         email = "palatejean@gmail.com"),
+  person("Alain", "Quartier-la-Tente", role = c("aut","cre"),
+         email = "alain.quartier@yahoo.fr",
+         comment = c(ORCID = "0000-0001-7890-3857")),
+  person("Tanguy", "Barthelemy", role = c("ctb"),
+         email ="tanguy.barthelemy@insee.fr"),
   person("Anna", "Smyk", role = c("ctb"),
-                 email ="anna.smyk@insee.fr")
-                 )
+         email ="anna.smyk@insee.fr")
+   )
 Description: This package provides functions to build and apply symmetric
   and asymmetric moving averages (= linear filters) for trend-cycle extraction.
   In particular, it implements several modern approaches for real-time estimates
@@ -29,17 +29,17 @@ Imports:
   MASS,
   graphics,
   stats,
-  rjd3toolkit (>= 3.2.2)
-Remotes: 
-    github::rjdverse/rjd3toolkit
+  rjd3toolkit (> 3.2.4)
+Remotes:
+  github::rjdverse/rjd3toolkit
 SystemRequirements: Java (>= 17)
-License: EUPL
+License: file LICENSE
 LazyData: TRUE
 URL: https://github.com/rjdverse/rjd3filters, https://rjdverse.github.io/rjd3filters/
 Suggests:
     knitr,
     rmarkdown
 VignetteBuilder: knitr
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 Roxygen: list(markdown = TRUE)
 Encoding: UTF-8
diff --git a/NAMESPACE b/NAMESPACE
index b4184cd..d99dda9 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -62,11 +62,13 @@ export(loocve)
 export(lower_bound)
 export(lp_filter)
 export(mirror)
+export(mmsre_filter)
 export(moving_average)
 export(mse)
 export(plot_coef)
 export(plot_gain)
 export(plot_phase)
+export(polynomial_matrix)
 export(rkhs_filter)
 export(rkhs_kernel)
 export(rkhs_optimal_bw)
diff --git a/NEWS.md b/NEWS.md
index 4a7cad1..bfaf50b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -7,8 +7,17 @@ to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
 
-## [2.1.1] - 2024-07-12
+## Added
 
+* New function `polynomial_matrix()` to create a matrix of polynomial terms.
+
+* New function `mmsre_filter()` to compute the general Proietti and Luati (2008) filter with extension for non symmetric filters and with Timeliness criterion.
+
+### Changed
+
+* `filter()` correction when the length of the series equals the length of the filter.
+
+## [2.1.1] - 2024-12-07
 
 ### Changed
 
diff --git a/R/2_finite_filters.R b/R/2_finite_filters.R
index 8ffe2cb..90bb0e0 100644
--- a/R/2_finite_filters.R
+++ b/R/2_finite_filters.R
@@ -46,7 +46,6 @@ finite_filters.moving_average <- function(sfilter,
     rfilters <- rev(lapply(lfilters, rev.moving_average))
   } else if (is.null(lfilters) && is.null(rfilters)) {
     rfilters <- lfilters <- list()
-
   }
   res <- new("finite_filters",
              sfilter = sfilter, lfilters = lfilters,
diff --git a/R/RKHS.R b/R/RKHS.R
index 91e531e..f38c385 100644
--- a/R/RKHS.R
+++ b/R/RKHS.R
@@ -89,13 +89,16 @@ rkhs_optimization_fun <- function(horizon = 6, leads = 0,  degree = 2,
                                smoothness = "Smoothness",
                                undefined = "Undefined")
   density <- match.arg(density)
-  optimalFunCriteria <- J("jdplus/filters/base/r/RKHSFilters")$optimalCriteria(
-    as.integer(horizon), as.integer(leads), as.integer(degree), kernel,
-    asymmetricCriterion, density=="rw", passband
-  )$applyAsDouble
-
+  jfun <-
+    .jcall(
+      "jdplus/filters/base/r/RKHSFilters",
+      "Ljava/util/function/DoubleUnaryOperator;",
+      "optimalCriteria",
+      as.integer(horizon), as.integer(leads), as.integer(degree), kernel,
+      asymmetricCriterion, density=="rw", passband
+    )
   Vectorize(function(x){
-    optimalFunCriteria(x)
+    .jcall(jfun, "D", "applyAsDouble", x)
   })
 }
 #' Optimal Bandwith of Reproducing Kernel Hilbert Space (RKHS) Filters
@@ -124,10 +127,14 @@ rkhs_optimal_bw <- function(horizon = 6,  degree = 2,
                                smoothness = "Smoothness",
                                undefined = "Undefined")
   density <- match.arg(density)
-  optimalBw <- J("jdplus/filters/base/r/RKHSFilters")$optimalBandwidth(
-    as.integer(horizon), as.integer(degree), kernel,
-    asymmetricCriterion, density=="rw", passband, optimal.minBandwidth, optimal.maxBandwidth
-  )
+  optimalBw <-
+    .jcall(
+      "jdplus/filters/base/r/RKHSFilters",
+      "[D",
+      "optimalBandwidth",
+      as.integer(horizon), as.integer(degree), kernel,
+      asymmetricCriterion, density=="rw", passband, optimal.minBandwidth, optimal.maxBandwidth
+    )
   names(optimalBw) <- sprintf("q=%i", 0:(horizon-1))
   optimalBw
 }
@@ -148,11 +155,14 @@ rkhs_kernel <- function(kernel = c("Biweight", "Henderson", "Epanechnikov", "Tri
     "epanechnikov" = "Epanechnikov",
     "henderson" = "Henderson"
   )
-  kernel_fun <- J("jdplus/filters/base/r/RKHSFilters")$kernel(
-    kernel, as.integer(degree), as.integer(horizon)
-  )$applyAsDouble
-
+  jfun <-
+    .jcall(
+      "jdplus/filters/base/r/RKHSFilters",
+      "Ljava/util/function/DoubleUnaryOperator;",
+      "kernel",
+      kernel, as.integer(degree), as.integer(horizon)
+    )
   Vectorize(function(x){
-    kernel_fun(x)
+    .jcall(jfun, "D", "applyAsDouble", x)
   })
 }
diff --git a/R/dfa.R b/R/dfa.R
index 50b24c6..64dc525 100644
--- a/R/dfa.R
+++ b/R/dfa.R
@@ -21,7 +21,7 @@
 #' dfa_filter(horizon = 6, degree = 2)
 dfa_filter <- function(horizon = 6, degree = 0,
                        density = c("uniform", "rw"),
-                       targetfilter = lp_filter(horizon = horizon)[,1],
+                       targetfilter = lp_filter(horizon = horizon)@sfilter,
                        passband = 2*pi/12,
                        accuracy.weight = 1/3,
                        smoothness.weight = 1/3,
@@ -37,7 +37,10 @@ dfa_filter <- function(horizon = 6, degree = 0,
       targetfilter <- coef(targetfilter)
     }
   }
-  dfa_filter <- J("jdplus/filters/base/r/DFAFilters")$filters(
+  dfa_filter <- .jcall(
+    "jdplus/filters/base/r/DFAFilters",
+    "Ljdplus/toolkit/base/core/math/linearfilters/ISymmetricFiltering;",
+    "filters",
     targetfilter,
     as.integer(horizon), as.integer(degree), density=="rw",
     passband,
diff --git a/R/filter.R b/R/filter.R
index dfcb2ca..223c675 100644
--- a/R/filter.R
+++ b/R/filter.R
@@ -75,15 +75,25 @@ filter_ma <- function(x, coefs){
   lb <- lower_bound(coefs)
   ub <- upper_bound(coefs)
 
-  if (length(x) <= length(coefs))
+  if (length(x) < length(coefs))
     return(x * NA)
 
-  DataBlock <- J("jdplus.toolkit.base.core.data.DataBlock")
-  jx <- DataBlock$of(as.numeric(x))
-  out <- DataBlock$of(as.numeric(rep(NA, length(x) - length(coefs)+1)))
-  .ma2jd(coefs)$apply(jx,
-                     out)
-  result <- out$toArray()
+  jx <- .jcall(
+    "jdplus/toolkit/base/core/data/DataBlock",
+    "Ljdplus/toolkit/base/core/data/DataBlock;",
+    "of",
+    as.numeric(x)
+    )
+  out <- .jcall(
+    "jdplus/toolkit/base/core/data/DataBlock",
+    "Ljdplus/toolkit/base/core/data/DataBlock;",
+    "of",
+    .jarray(as.numeric(rep(NA, length(x) - length(coefs)+1)))
+  )
+  jfilter <- .ma2jd(coefs)
+  .jcall(jfilter, "V", "apply",
+    jx, out)
+  result <- .jcall(out, "[D", "toArray")
   result <- c(rep(NA, abs(min(lb, 0))),
               result,
               rep(NA, abs(max(ub, 0))))
diff --git a/R/get_properties_function.R b/R/get_properties_function.R
index 0d949c4..1710ef2 100644
--- a/R/get_properties_function.R
+++ b/R/get_properties_function.R
@@ -21,22 +21,29 @@ get_properties_function <- function(x,
 }
 
 get_gain_function <- function(x){
-  jgain <- x$gainFunction()$applyAsDouble
+  jgain <- .jcall(x, "Ljava/util/function/DoubleUnaryOperator;",
+                  "gainFunction")
   Vectorize(function(x){
-    jgain(x)
+    .jcall(jgain, "D", "applyAsDouble", x)
   })
 }
 get_phase_function <- function(x){
-  jphase <- x$phaseFunction()$applyAsDouble
+  jphase <- .jcall(x, "Ljava/util/function/DoubleUnaryOperator;",
+                   "phaseFunction")
   Vectorize(function(x){
-    jphase(x)
+    .jcall(jphase, "D", "applyAsDouble", x)
   })
 }
 get_frequency_response_function <- function(x){
-  jfrf <- x$frequencyResponseFunction()$apply
+  jfrf <- .jcall(x,
+                 "Ljava/lang/Object;",
+                 "frequencyResponseFunction")
+
   Vectorize(function(x){
-    res <- jfrf(x)
-    complex(real = res$getRe(), imaginary = res$getIm())
+    res <- .jcall(jfrf, "Ljava/lang/Object;", "apply", x)
+
+    complex(real = .jcall(res, "D", "getRe"),
+            imaginary = .jcall(res, "D", "getIm"))
   })
 }
 
diff --git a/R/kernels.R b/R/kernels.R
index a561ca6..f151a34 100644
--- a/R/kernels.R
+++ b/R/kernels.R
@@ -16,11 +16,28 @@ get_kernel <- function(kernel = c("Henderson","Uniform", "Triangular",
                                   "Trapezoidal", "Gaussian"),
                        horizon,
                        sd_gauss = 0.25){
+  jkernel <- .r2jd_kernel(kernel, horizon, sd_gauss)
+  coef <- sapply(as.integer(seq.int(from = 0, to = horizon, by = 1)),
+                function(x) .jcall(jkernel, "D", "applyAsDouble", x))
+  m <- horizon
+  result <- list(coef = coef, m = m)
+  attr(result, "name") <- kernel
+  attr(result, "class") <- "tskernel"
+  result
+}
+.r2jd_kernel <- function(
+    kernel = c("Henderson","Uniform", "Triangular",
+               "Epanechnikov","Parabolic","BiWeight", "TriWeight","Tricube",
+               "Trapezoidal", "Gaussian"),
+    horizon, sd_gauss = 0.25){
+
+  if (is.null(kernel) || kernel[1]=="")
+    return(.jnull("java/util/function/IntToDoubleFunction"))
 
   kernel <- match.arg(tolower(kernel)[1],
-                     choices = c("henderson", "uniform", "triangular", "epanechnikov", "parabolic",
-                                 "biweight", "triweight", "tricube", "trapezoidal", "gaussian"
-                     ))
+                      choices = c("henderson", "uniform", "triangular", "epanechnikov", "parabolic",
+                                  "biweight", "triweight", "tricube", "trapezoidal", "gaussian"
+                      ))
   if (kernel == "parabolic")
     kernel <- "epanechnikov"
   h <- as.integer(horizon)
@@ -33,12 +50,5 @@ get_kernel <- function(kernel = c("Henderson","Uniform", "Triangular",
                       "Ljava/util/function/IntToDoubleFunction;",
                       tolower(kernel), h)
   }
-
-  coef <- sapply(as.integer(seq.int(from = 0, to = horizon, by = 1)),
-                jkernel$applyAsDouble)
-  m <- horizon
-  result <- list(coef = coef, m = m)
-  attr(result, "name") <- kernel
-  attr(result, "class") <- "tskernel"
-  result
+  jkernel
 }
diff --git a/R/lp_filters.R b/R/lp_filters.R
index 12a5922..3d23126 100644
--- a/R/lp_filters.R
+++ b/R/lp_filters.R
@@ -7,7 +7,7 @@ NULL
 #' @param horizon horizon (bandwidth) of the symmetric filter.
 #' @param degree degree of polynomial.
 #' @param kernel kernel uses.
-#' @param endpoints methode for endpoints.
+#' @param endpoints method for endpoints.
 #' @param tweight timeliness weight.
 #' @param passband passband threshold.
 #' @param ic ic ratio.
@@ -67,7 +67,7 @@ localpolynomials<-function(x,
 #' * "CN": Cut and Normalized Filter
 #'
 #' @return a [finite_filters()] object.
-#' @seealso [localpolynomials()].
+#' @seealso [mmsre_filter()] [localpolynomials()].
 #' @examples
 #' henderson_f <- lp_filter(horizon = 6, kernel = "Henderson")
 #' plot_coef(henderson_f)
diff --git a/R/mmsreFilter.R b/R/mmsreFilter.R
new file mode 100644
index 0000000..b34d302
--- /dev/null
+++ b/R/mmsreFilter.R
@@ -0,0 +1,107 @@
+#' Mean Square Revision Error (mmsre) filter
+#'
+#' Provides an asymmetric filter based on the given reference
+#' filter (usually symmetric) minimizing the mean square revision error.
+#'
+#' @param ref_filter The reference filter (a [moving_average()] object).
+#' @param q The horizon of the asymmetric filter.
+#' @param U Matrix of the constraints.
+#' @param Z Matrix of the bias (can be `NULL`).
+#' @param delta Coefficients of the linear model.
+#' @param kernel The kernel used for weighting factors, by default, no weight is used.
+#' See [lp_filter()] for the available kernels.
+#'
+#' @details
+#' The asymmetric filter \eqn{\boldsymbol v=(v_{-h},\dots,v{q})'} minimizes the mean square revision error
+#' (mmsre) relative to the reference filter \eqn{\boldsymbol \theta=(\theta_{-h},\dots,\theta_{h'})'}.
+#' The series follows the model
+#' \deqn{
+#' \boldsymbol y=\boldsymbol U \boldsymbol \gamma +
+#' \boldsymbol Z \boldsymbol \delta + \boldsymbol \varepsilon, \quad
+#' \boldsymbol \varepsilon \sim \mathcal N(0,\sigma^2 \boldsymbol K^{-1}).
+#' }
+#'
+#' With \eqn{K} a set of weights (kernel), by default (`kernel = NULL`) no weight is used.
+#' The matrix \eqn{U} represents the constraints of the symmetric filter (usually polynomials preservations), \eqn{\boldsymbol \theta},
+#' imposed to the asymmetric filter, \eqn{\boldsymbol v}.
+#' Partitionning the matrix \eqn{\boldsymbol U=\begin{pmatrix} \boldsymbol U_p' & \boldsymbol U_f'\end{pmatrix}'}
+#' with \eqn{\boldsymbol U_p} the first \eqn{h+q+1} rows and \eqn{\boldsymbol U_f} the remaining rows, the constraints are
+#' \eqn{\boldsymbol U_p'\boldsymbol v=\boldsymbol U'\boldsymbol \theta}.
+#'
+#' The matrix \eqn{\boldsymbol Z} represents the bias of the asymmetric filter: usually constraints imposed to the symmetric filter but not to the asymmetric filter.
+#'
+#'
+#'
+#'
+#' @inheritParams lp_filter
+#' @examples
+#' QL <- lp_filter(endpoints = "QL", ic = 3.5)
+#' LC <- lp_filter(endpoints = "LC", ic = 3.5)
+#' DAF <- lp_filter(endpoints = "DAF")
+#' h6 <- QL[, "q=6"]
+#' # To reproduce DAF filter
+#' mmsre_filter(
+#'   ref_filter = h6, q = 0,
+#'   U = polynomial_matrix(l = - 6, d0 = 0, d1 = 3),
+#'   kernel = "Henderson"
+#' )
+#' DAF[, "q=0"]
+#' # To reproduce QL filter
+#' mmsre_filter(
+#'   ref_filter = h6, q = 1,
+#'   delta = 2 / (sqrt(pi) * 3.5),
+#'   U = polynomial_matrix(l = -6, d0 = 0, d1 = 1),
+#'   Z = polynomial_matrix(l = -6, d0 = 2, d1 = 2)
+#' )
+#' QL[, "q=1"]
+#'
+#' # Or using the Uniform kernel
+#' mmsre_filter(
+#'   ref_filter = h6, q = 2,
+#'   # we multiply by the square root of the inverse of weights (1/13)
+#'   # to get the same result as the QL filter
+#'   delta = 2 / (sqrt(pi) * 3.5) * (sqrt(13)),
+#'   U = polynomial_matrix(l = -6, d0 = 0, d1 = 0),
+#'   Z = polynomial_matrix(l = -6, d0 = 1, d1 = 1),
+#'   kernel = "Uniform"
+#' )
+#' LC[, "q=2"]
+#' @seealso [lp_filter()].
+#' @references Proietti, Tommaso and Alessandra Luati (2008). “Real time estimation in local polynomial regression, with application to trend-cycle analysis”.
+#' @export
+mmsre_filter <- function(
+    ref_filter, q, U, Z = NULL, delta = NULL,
+    kernel = NULL,
+    tweight = 0, passband = pi/12){
+  jref <- .jcast(.ma2jd(ref_filter), "jdplus/toolkit/base/core/math/linearfilters/IFiniteFilter")
+  if (is.null(delta))
+    delta <- numeric()
+
+  jkernel <- .r2jd_kernel(kernel, abs(.jcall(jref, "I", "getLowerBound")))
+  jf <- .jcall(
+    "jdplus/toolkit/base/core/math/linearfilters/AsymmetricFiltersFactory",
+    "Ljdplus/toolkit/base/core/math/linearfilters/IFiniteFilter;",
+    "mmsreFilter",
+    jref,
+    as.integer(q),
+    .r2jd_fast_matrix(U),
+    .r2jd_fast_matrix(Z),
+    .jarray(delta),
+    jkernel,
+    as.numeric(tweight),
+    as.numeric(passband)
+  )
+  return(.jd2ma(jf))
+}
+
+.r2jd_fast_matrix <- function(s){
+  if (is.null(s))
+    return(.jnull("jdplus/toolkit/base/core/math/matrices/FastMatrix"))
+
+  .jcall(
+    "jdplus/toolkit/base/core/math/matrices/FastMatrix",
+    "Ljdplus/toolkit/base/core/math/matrices/FastMatrix;",
+    "of",
+    rjd3toolkit::.r2jd_matrix(s))
+}
+
diff --git a/R/plots.R b/R/plots.R
index 1f5a36d..ecec798 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -13,6 +13,7 @@
 #' @param normalized boolean indicatif if the phase function is normalized by the frequency.
 #' @param zero_as_na boolean indicating if the trailing zero of the coefficients should be plotted (\code{FALSE}) or removed (\code{TRUE}).
 #' @param n number of points used to plot the functions.
+#' @param xlab,ylab labels of axis.
 #'
 #' @examples
 #' filter <- lp_filter(6, endpoints = "DAF", kernel = "Henderson")
@@ -23,14 +24,14 @@
 #' @rdname plot_filters
 #' @importFrom MASS fractions
 #' @export
-plot_coef <- function(x, nxlab = 7, add = FALSE, ...){
+plot_coef <- function(x, nxlab = 7, add = FALSE, ..., xlab = "", ylab = "coefficient"){
   UseMethod("plot_coef", x)
 }
 #' @rdname plot_filters
 #' @export
 plot_coef.default <- function(x, nxlab = 7, add = FALSE,
                               zero_as_na = TRUE, q = 0, legend = FALSE,
-                              legend.pos = "topright", ...){
+                              legend.pos = "topright", ..., xlab = "", ylab = "coefficient"){
   if (zero_as_na)
     x  <- apply(x,2, trailingzero_as_na)
   col_to_plot <- sprintf("q=%i",q)
@@ -38,16 +39,16 @@ plot_coef.default <- function(x, nxlab = 7, add = FALSE,
   horizon <- (nrow(x)-1)/2
   if (length(col_to_plot) == 0) {
     if (!add){
-      plot(1, type="n",xaxt = "n", xlab = "",
-           ylab = "coefficient", xlim=c(-horizon, horizon), ylim=c(0, 1),
+      plot(1, type="n",xaxt = "n", xlab = xlab,
+           ylab = ylab, xlim=c(-horizon, horizon), ylim=c(0, 1),
            ...)
       axis(1, at=seq(-horizon, horizon, by = 1), labels = rownames(x))
     }
     return(invisible(0))
   }
   matplot(seq(-horizon, horizon, by = 1),x[,col_to_plot],
-          xaxt = "n", xlab = "", type = "o", pch = 20,
-          ylab = "coefficient", add = add, ...)
+          xaxt = "n", xlab = xlab, type = "o", pch = 20,
+          ylab = ylab, add = add, ...)
   if (legend)
     legend(legend.pos,col_to_plot,
            col = seq_along(col_to_plot), lty=seq_along(col_to_plot), lwd=2)
@@ -57,11 +58,11 @@ plot_coef.default <- function(x, nxlab = 7, add = FALSE,
 
 #' @rdname plot_filters
 #' @export
-plot_coef.moving_average <- function(x, nxlab = 7, add = FALSE, ...){
+plot_coef.moving_average <- function(x, nxlab = 7, add = FALSE, ..., xlab = "", ylab = "coefficient"){
   x_plot <- coef(x)
   matplot(seq(lower_bound(x), upper_bound(x), by = 1), x_plot,
-          xaxt = "n", xlab = "", type = "o", pch = 20,
-          ylab = "coefficient", add = add, ...)
+          xaxt = "n", xlab = xlab, type = "o", pch = 20,
+          ylab = ylab, add = add, ...)
   if (!add)
     axis(1, at=seq(lower_bound(x), upper_bound(x), by = 1), labels = names(x_plot))
 }
@@ -70,27 +71,28 @@ plot_coef.moving_average <- function(x, nxlab = 7, add = FALSE, ...){
 #' @export
 plot_coef.finite_filters <- function(x, nxlab = 7, add = FALSE,
                                      zero_as_na = TRUE, q = 0, legend = length(q) > 1,
-                                     legend.pos = "topright", ...){
+                                     legend.pos = "topright", ..., xlab = "", ylab = "coefficient"){
   plot_coef(x = as.matrix(x, zero_as_na = zero_as_na),
             nxlab = nxlab, add = add,
             zero_as_na = FALSE, q = q, legend = legend,
-            legend.pos = legend.pos, ...)
+            legend.pos = legend.pos, ...,
+            xlab = xlab, ylab = ylab)
 }
 
 #' @rdname plot_filters
 #' @export
 plot_gain <- function(x, nxlab = 7, add = FALSE,
-                      xlim = c(0, pi), ...){
+                      xlim = c(0, pi), ..., xlab = "", ylab = "gain"){
   UseMethod("plot_gain", x)
 }
 #' @rdname plot_filters
 #' @export
 plot_gain.moving_average<- function(x, nxlab = 7, add = FALSE,
-                                    xlim = c(0, pi), ...){
+                                    xlim = c(0, pi), ..., xlab = "", ylab = "gain"){
   g <- get_properties_function(x, "Symmetric Gain")
   plot(g, type = "l",
-       xaxt = "n", xlab = "",
-       ylab = "gain", add = add, xlim = xlim, ...)
+       xaxt = "n", xlab = xlab,
+       ylab = ylab, add = add, xlim = xlim, ...)
   if (!add){
     x_lab_at <- seq(xlim[1]/pi, xlim[2]/pi, length.out = nxlab)
     axis(1, at = x_lab_at * pi, labels = xlabel(x_lab_at))
@@ -101,7 +103,7 @@ plot_gain.moving_average<- function(x, nxlab = 7, add = FALSE,
 plot_gain.finite_filters <- function(x, nxlab = 7, add = FALSE,
                                 xlim = c(0, pi), q = 0, legend = length(q) > 1,
                                 legend.pos = "topright",
-                                n = 101, ...){
+                                n = 101, ..., xlab = "", ylab = "gain"){
   x_values <- seq.int(xlim[1], xlim[2], length.out = n)
   gsym <- get_properties_function(x, "Symmetric Gain")
   gasym <- get_properties_function(x, "Asymmetric Gain")
@@ -113,8 +115,8 @@ plot_gain.finite_filters <- function(x, nxlab = 7, add = FALSE,
   y_val <- sapply(all_g_f, function(f) f(x_values))
   if (length(col_to_plot) == 0) {
     if (!add){
-      plot(1, type="n",xaxt = "n", xlab = "",
-           ylab = "gain", xlim=xlim, ylim=c(0, 1),
+      plot(1, type="n",xaxt = "n", xlab = xlab,
+           ylab = ylab, xlim=xlim, ylim=c(0, 1),
            ...)
       x_lab_at <- seq(xlim[1]/pi, xlim[2]/pi, length.out = nxlab)
       axis(1, at = x_lab_at * pi, labels = xlabel(x_lab_at))
@@ -122,8 +124,8 @@ plot_gain.finite_filters <- function(x, nxlab = 7, add = FALSE,
     return(invisible(0))
   }
   matplot(x_values, y_val[, col_to_plot], type = "l",
-          xaxt = "n", xlab = "",
-          ylab = "gain", add = add, xlim = xlim, ...)
+          xaxt = "n", xlab = xlab,
+          ylab = ylab, add = add, xlim = xlim, ...)
 
   if (legend)
     legend(legend.pos,col_to_plot,
@@ -137,13 +139,14 @@ plot_gain.finite_filters <- function(x, nxlab = 7, add = FALSE,
 #' @rdname plot_filters
 #' @export
 plot_phase <- function(x, nxlab = 7, add = FALSE,
-                       xlim = c(0, pi), normalized = FALSE, ...){
+                       xlim = c(0, pi), normalized = FALSE, ..., xlab = "", ylab = "phase"){
   UseMethod("plot_phase", x)
 }
 #' @rdname plot_filters
 #' @export
 plot_phase.moving_average<- function(x, nxlab = 7, add = FALSE,
-                                     xlim = c(0, pi), normalized = FALSE, ...){
+                                     xlim = c(0, pi), normalized = FALSE, ...,
+                                     xlab = "", ylab = "phase"){
   p <- get_properties_function(x, "Symmetric Phase")
 
   if (normalized) {
@@ -155,8 +158,8 @@ plot_phase.moving_average<- function(x, nxlab = 7, add = FALSE,
   }
 
   plot(p_plot, type = "l",
-       xaxt = "n", xlab = "",
-       ylab = "phase", add = add, xlim = xlim, ...)
+       xaxt = "n", xlab = xlab,
+       ylab = ylab, add = add, xlim = xlim, ...)
   if (!add){
     x_lab_at <- seq(xlim[1]/pi, xlim[2]/pi, length.out = nxlab)
     axis(1, at = x_lab_at * pi, labels = xlabel(x_lab_at))
@@ -169,7 +172,8 @@ plot_phase.moving_average<- function(x, nxlab = 7, add = FALSE,
 plot_phase.finite_filters <- function(x, nxlab = 7, add = FALSE,
                                       xlim = c(0, pi), normalized = FALSE,
                                       q = 0, legend = length(q) > 1, legend.pos = "topright",
-                                     n = 101, ...){
+                                     n = 101, ...,
+                                     xlab = "", ylab = "phase"){
   x_values <- seq.int(xlim[1], xlim[2], length.out = n)
   psym <- get_properties_function(x, "Symmetric Phase")
   pasym <- get_properties_function(x, "Asymmetric Phase")
@@ -185,8 +189,8 @@ plot_phase.finite_filters <- function(x, nxlab = 7, add = FALSE,
   }
   if (length(col_to_plot) == 0) {
     if (!add){
-      plot(1, type="n",xaxt = "n", xlab = "",
-           ylab = "phase", xlim=xlim, ylim=c(0, 1),
+      plot(1, type="n",xaxt = "n", xlab = xlab,
+           ylab = ylab, xlim=xlim, ylim=c(0, 1),
            ...)
       x_lab_at <- seq(xlim[1]/pi, xlim[2]/pi, length.out = nxlab)
       axis(1, at = x_lab_at * pi, labels = xlabel(x_lab_at))
@@ -194,8 +198,8 @@ plot_phase.finite_filters <- function(x, nxlab = 7, add = FALSE,
     return(invisible(0))
   }
   matplot(x_values, y_val[, col_to_plot], type = "l",
-          xaxt = "n", xlab = "",
-          ylab = "phase", add = add, xlim = xlim, ...)
+          xaxt = "n", xlab = xlab,
+          ylab = ylab, add = add, xlim = xlim, ...)
 
   if (legend)
     legend(legend.pos,col_to_plot,
diff --git a/R/polynomial_matrix.R b/R/polynomial_matrix.R
new file mode 100644
index 0000000..7ffcbbe
--- /dev/null
+++ b/R/polynomial_matrix.R
@@ -0,0 +1,40 @@
+#' Create polynomial matrix
+#'
+#' Create polynomial matrix used in local polynomial regression (see details).
+#'
+#' @param l,u lower bound (usually negative) and upper bound (usually positive) of the polynomial matrix.
+#' @param d0,d1 lower and polynomial degree of the polynomial matrix.
+#'
+#' @details
+#' `polynomial_matrix()` computes the following matrix
+#' \deqn{
+#' \begin{pmatrix}
+#' (l)^{d_0} & (l)^{d_0+1} & \cdots&(l)^{d_1}\\
+#' (l+1)^{d_0} & (l+1)^{d_0+1} & \cdots&(l+1)^{d_1} \\
+#' \vdots & \vdots & \cdots & \vdots \\
+#' (p)^{d_0} & (p)^{d_0+1} & \cdots&(p)^{d_1}
+#' \end{pmatrix}
+#' }
+#'
+#' @examples
+#' # For example to reproduce DAF filters
+#' daf <- lp_filter(horizon = 6, endpoints = "DAF")
+#' q <- 0
+#' X <- polynomial_matrix(l = -6, u = q, d0 = 0, d1 = 3)
+#' K <- diag(sapply(-6:q, function(i) get_kernel(horizon = 6)[i]))
+#' e_1 <- c(1, 0, 0, 0)
+#' q0 <- K %*% X %*% solve(t(X) %*% K %*% X, e_1)
+#' q0
+#' daf[, "q=0"]
+#' @export
+polynomial_matrix <- function(l, u = abs(l), d0 = 0, d1 = 3){
+  .jd2r_matrix(
+    .jcall(
+      "jdplus/toolkit/base/core/math/linearfilters/LocalPolynomialFilters",
+      "Ljdplus/toolkit/base/core/math/matrices/FastMatrix;",
+      "z",
+      .jnull("jdplus/toolkit/base/core/math/matrices/FastMatrix"),
+      as.integer(l), as.integer(u),
+      as.integer(d0), as.integer(d1))
+  )
+}
diff --git a/man/dfa_filter.Rd b/man/dfa_filter.Rd
index 6f68e82..269586f 100644
--- a/man/dfa_filter.Rd
+++ b/man/dfa_filter.Rd
@@ -8,7 +8,7 @@ dfa_filter(
   horizon = 6,
   degree = 0,
   density = c("uniform", "rw"),
-  targetfilter = lp_filter(horizon = horizon)[, 1],
+  targetfilter = lp_filter(horizon = horizon)@sfilter,
   passband = 2 * pi/12,
   accuracy.weight = 1/3,
   smoothness.weight = 1/3,
diff --git a/man/localpolynomials.Rd b/man/localpolynomials.Rd
index b4611f2..66c5497 100644
--- a/man/localpolynomials.Rd
+++ b/man/localpolynomials.Rd
@@ -25,7 +25,7 @@ localpolynomials(
 
 \item{kernel}{kernel uses.}
 
-\item{endpoints}{methode for endpoints.}
+\item{endpoints}{method for endpoints.}
 
 \item{ic}{ic ratio.}
 
diff --git a/man/lp_filter.Rd b/man/lp_filter.Rd
index 50ac419..ae076c8 100644
--- a/man/lp_filter.Rd
+++ b/man/lp_filter.Rd
@@ -22,7 +22,7 @@ lp_filter(
 
 \item{kernel}{kernel uses.}
 
-\item{endpoints}{methode for endpoints.}
+\item{endpoints}{method for endpoints.}
 
 \item{ic}{ic ratio.}
 
@@ -54,5 +54,5 @@ plot_coef(henderson_f)
 Proietti, Tommaso and Alessandra Luati (2008). “Real time estimation in local polynomial regression, with application to trend-cycle analysis”.
 }
 \seealso{
-\code{\link[=localpolynomials]{localpolynomials()}}.
+\code{\link[=mmsre_filter]{mmsre_filter()}} \code{\link[=localpolynomials]{localpolynomials()}}.
 }
diff --git a/man/mmsre_filter.Rd b/man/mmsre_filter.Rd
new file mode 100644
index 0000000..571ff5f
--- /dev/null
+++ b/man/mmsre_filter.Rd
@@ -0,0 +1,97 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mmsreFilter.R
+\name{mmsre_filter}
+\alias{mmsre_filter}
+\title{Mean Square Revision Error (mmsre) filter}
+\usage{
+mmsre_filter(
+  ref_filter,
+  q,
+  U,
+  Z = NULL,
+  delta = NULL,
+  kernel = NULL,
+  tweight = 0,
+  passband = pi/12
+)
+}
+\arguments{
+\item{ref_filter}{The reference filter (a \code{\link[=moving_average]{moving_average()}} object).}
+
+\item{q}{The horizon of the asymmetric filter.}
+
+\item{U}{Matrix of the constraints.}
+
+\item{Z}{Matrix of the bias (can be \code{NULL}).}
+
+\item{delta}{Coefficients of the linear model.}
+
+\item{kernel}{The kernel used for weighting factors, by default, no weight is used.
+See \code{\link[=lp_filter]{lp_filter()}} for the available kernels.}
+
+\item{tweight}{timeliness weight.}
+
+\item{passband}{passband threshold.}
+}
+\description{
+Provides an asymmetric filter based on the given reference
+filter (usually symmetric) minimizing the mean square revision error.
+}
+\details{
+The asymmetric filter \eqn{\boldsymbol v=(v_{-h},\dots,v{q})'} minimizes the mean square revision error
+(mmsre) relative to the reference filter \eqn{\boldsymbol \theta=(\theta_{-h},\dots,\theta_{h'})'}.
+The series follows the model
+\deqn{
+\boldsymbol y=\boldsymbol U \boldsymbol \gamma +
+\boldsymbol Z \boldsymbol \delta + \boldsymbol \varepsilon, \quad
+\boldsymbol \varepsilon \sim \mathcal N(0,\sigma^2 \boldsymbol K^{-1}).
+}
+
+With \eqn{K} a set of weights (kernel), by default (\code{kernel = NULL}) no weight is used.
+The matrix \eqn{U} represents the constraints of the symmetric filter (usually polynomials preservations), \eqn{\boldsymbol \theta},
+imposed to the asymmetric filter, \eqn{\boldsymbol v}.
+Partitionning the matrix \eqn{\boldsymbol U=\begin{pmatrix} \boldsymbol U_p' & \boldsymbol U_f'\end{pmatrix}'}
+with \eqn{\boldsymbol U_p} the first \eqn{h+q+1} rows and \eqn{\boldsymbol U_f} the remaining rows, the constraints are
+\eqn{\boldsymbol U_p'\boldsymbol v=\boldsymbol U'\boldsymbol \theta}.
+
+The matrix \eqn{\boldsymbol Z} represents the bias of the asymmetric filter: usually constraints imposed to the symmetric filter but not to the asymmetric filter.
+}
+\examples{
+QL <- lp_filter(endpoints = "QL", ic = 3.5)
+LC <- lp_filter(endpoints = "LC", ic = 3.5)
+DAF <- lp_filter(endpoints = "DAF")
+h6 <- QL[, "q=6"]
+# To reproduce DAF filter
+mmsre_filter(
+  ref_filter = h6, q = 0,
+  U = polynomial_matrix(l = - 6, d0 = 0, d1 = 3),
+  kernel = "Henderson"
+)
+DAF[, "q=0"]
+# To reproduce QL filter
+mmsre_filter(
+  ref_filter = h6, q = 1,
+  delta = 2 / (sqrt(pi) * 3.5),
+  U = polynomial_matrix(l = -6, d0 = 0, d1 = 1),
+  Z = polynomial_matrix(l = -6, d0 = 2, d1 = 2)
+)
+QL[, "q=1"]
+
+# Or using the Uniform kernel
+mmsre_filter(
+  ref_filter = h6, q = 2,
+  # we multiply by the square root of the inverse of weights (1/13)
+  # to get the same result as the QL filter
+  delta = 2 / (sqrt(pi) * 3.5) * (sqrt(13)),
+  U = polynomial_matrix(l = -6, d0 = 0, d1 = 0),
+  Z = polynomial_matrix(l = -6, d0 = 1, d1 = 1),
+  kernel = "Uniform"
+)
+LC[, "q=2"]
+}
+\references{
+Proietti, Tommaso and Alessandra Luati (2008). “Real time estimation in local polynomial regression, with application to trend-cycle analysis”.
+}
+\seealso{
+\code{\link[=lp_filter]{lp_filter()}}.
+}
diff --git a/man/plot_filters.Rd b/man/plot_filters.Rd
index 293f975..b4e2fa7 100644
--- a/man/plot_filters.Rd
+++ b/man/plot_filters.Rd
@@ -14,7 +14,7 @@
 \alias{plot_phase.finite_filters}
 \title{Plots filters properties}
 \usage{
-plot_coef(x, nxlab = 7, add = FALSE, ...)
+plot_coef(x, nxlab = 7, add = FALSE, ..., xlab = "", ylab = "coefficient")
 
 \method{plot_coef}{default}(
   x,
@@ -24,10 +24,12 @@ plot_coef(x, nxlab = 7, add = FALSE, ...)
   q = 0,
   legend = FALSE,
   legend.pos = "topright",
-  ...
+  ...,
+  xlab = "",
+  ylab = "coefficient"
 )
 
-\method{plot_coef}{moving_average}(x, nxlab = 7, add = FALSE, ...)
+\method{plot_coef}{moving_average}(x, nxlab = 7, add = FALSE, ..., xlab = "", ylab = "coefficient")
 
 \method{plot_coef}{finite_filters}(
   x,
@@ -37,12 +39,30 @@ plot_coef(x, nxlab = 7, add = FALSE, ...)
   q = 0,
   legend = length(q) > 1,
   legend.pos = "topright",
-  ...
+  ...,
+  xlab = "",
+  ylab = "coefficient"
 )
 
-plot_gain(x, nxlab = 7, add = FALSE, xlim = c(0, pi), ...)
+plot_gain(
+  x,
+  nxlab = 7,
+  add = FALSE,
+  xlim = c(0, pi),
+  ...,
+  xlab = "",
+  ylab = "gain"
+)
 
-\method{plot_gain}{moving_average}(x, nxlab = 7, add = FALSE, xlim = c(0, pi), ...)
+\method{plot_gain}{moving_average}(
+  x,
+  nxlab = 7,
+  add = FALSE,
+  xlim = c(0, pi),
+  ...,
+  xlab = "",
+  ylab = "gain"
+)
 
 \method{plot_gain}{finite_filters}(
   x,
@@ -53,12 +73,32 @@ plot_gain(x, nxlab = 7, add = FALSE, xlim = c(0, pi), ...)
   legend = length(q) > 1,
   legend.pos = "topright",
   n = 101,
-  ...
+  ...,
+  xlab = "",
+  ylab = "gain"
 )
 
-plot_phase(x, nxlab = 7, add = FALSE, xlim = c(0, pi), normalized = FALSE, ...)
+plot_phase(
+  x,
+  nxlab = 7,
+  add = FALSE,
+  xlim = c(0, pi),
+  normalized = FALSE,
+  ...,
+  xlab = "",
+  ylab = "phase"
+)
 
-\method{plot_phase}{moving_average}(x, nxlab = 7, add = FALSE, xlim = c(0, pi), normalized = FALSE, ...)
+\method{plot_phase}{moving_average}(
+  x,
+  nxlab = 7,
+  add = FALSE,
+  xlim = c(0, pi),
+  normalized = FALSE,
+  ...,
+  xlab = "",
+  ylab = "phase"
+)
 
 \method{plot_phase}{finite_filters}(
   x,
@@ -70,7 +110,9 @@ plot_phase(x, nxlab = 7, add = FALSE, xlim = c(0, pi), normalized = FALSE, ...)
   legend = length(q) > 1,
   legend.pos = "topright",
   n = 101,
-  ...
+  ...,
+  xlab = "",
+  ylab = "phase"
 )
 }
 \arguments{
@@ -82,6 +124,8 @@ plot_phase(x, nxlab = 7, add = FALSE, xlim = c(0, pi), normalized = FALSE, ...)
 
 \item{...}{other arguments to \code{matplot}.}
 
+\item{xlab, ylab}{labels of axis.}
+
 \item{zero_as_na}{boolean indicating if the trailing zero of the coefficients should be plotted (\code{FALSE}) or removed (\code{TRUE}).}
 
 \item{q}{q.}
diff --git a/man/polynomial_matrix.Rd b/man/polynomial_matrix.Rd
new file mode 100644
index 0000000..6acbbac
--- /dev/null
+++ b/man/polynomial_matrix.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/polynomial_matrix.R
+\name{polynomial_matrix}
+\alias{polynomial_matrix}
+\title{Create polynomial matrix}
+\usage{
+polynomial_matrix(l, u = abs(l), d0 = 0, d1 = 3)
+}
+\arguments{
+\item{l, u}{lower bound (usually negative) and upper bound (usually positive) of the polynomial matrix.}
+
+\item{d0, d1}{lower and polynomial degree of the polynomial matrix.}
+}
+\description{
+Create polynomial matrix used in local polynomial regression (see details).
+}
+\details{
+\code{polynomial_matrix()} computes the following matrix
+\deqn{
+\begin{pmatrix}
+(l)^{d_0} & (l)^{d_0+1} & \cdots&(l)^{d_1}\\
+(l+1)^{d_0} & (l+1)^{d_0+1} & \cdots&(l+1)^{d_1} \\
+\vdots & \vdots & \cdots & \vdots \\
+(p)^{d_0} & (p)^{d_0+1} & \cdots&(p)^{d_1}
+\end{pmatrix}
+}
+}
+\examples{
+# For example to reproduce DAF filters
+daf <- lp_filter(horizon = 6, endpoints = "DAF")
+q <- 0
+X <- polynomial_matrix(l = -6, u = q, d0 = 0, d1 = 3)
+K <- diag(sapply(-6:q, function(i) get_kernel(horizon = 6)[i]))
+e_1 <- c(1, 0, 0, 0)
+q0 <- K \%*\% X \%*\% solve(t(X) \%*\% K \%*\% X, e_1)
+q0
+daf[, "q=0"]
+}