Repository: alan-turing-institute/sktime
Pull request #3381 (2022-09-04 16:15:15)

[ENH] enable default `get_fitted_params` for forecasters and delegated estimators
This PR enables a default implementation of `get_fitted_params` for forecasters and delegated estimators:

* for forecasting, `BaseForecaster` accidentally overrode the recently introduced default implementation of `get_fitted_params`. The abstract method has been removed from `BaseForecaster` now, which should enable the default in all forecasters. (all other estimator types were unaffected)
* for delegators, `_get_fitted_params` now delegates to the delegate's `get_fitted_params`, in line with all other delegation defaults.
* the default `get_fitted_params` is robustified so that it queries fitted parameters (recursion) only from fitted components
* a conditional override is added to `BaseForecaster` and `BaseTransformer` for the vectorized case (to avoid unfitted estimators present in underscore attributes)

Also see https://github.com/alan-turing-institute/sktime/issues/1497.

This should solve #1497 for all compositors that are delegates.

Relies on:
* #3391 for solving vectorization/`get_fitted_params` clashes
* #3393 for solving vectorization/`get_fitted_params` clashes

fkiraly: Note to self: works, but the `BaseForecaster` override masked a bug where vectorization and the default interact.
In case of vectorization, components may remain uninitialized/unfitted, and the vectorized components do not get returned via `get_fitted_params`.

Patch:
diff --git a/docs/source/api_reference/transformations.rst b/docs/source/api_reference/transformations.rst
index e0e3174d7..0a0727261 100644
--- a/docs/source/api_reference/transformations.rst
+++ b/docs/source/api_reference/transformations.rst
@@ -382,8 +382,8 @@ Missing value imputation
 
     Imputer
 
-Datetime feature generation
-~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Seasonality and Date-Time Features
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 .. currentmodule:: sktime.transformations.series.date
 
@@ -393,6 +393,14 @@ Datetime feature generation
 
     DateTimeFeatures
 
+.. currentmodule:: sktime.transformations.series.fourier
+
+.. autosummary::
+    :toctree: auto_generated/
+    :template: class.rst
+
+    FourierFeatures
+
 Auto-correlation series
 ~~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/extension_templates/transformer.py b/extension_templates/transformer.py
index 3d08b09af..2f224060e 100644
--- a/extension_templates/transformer.py
+++ b/extension_templates/transformer.py
@@ -99,35 +99,150 @@ class MyTransformer(BaseTransformer):
     #   y_inner_mtype must be changed to one or a list of compatible sktime mtypes
     #  the other tags are "safe defaults" which can usually be left as-is
     _tags = {
+        # to list all valid tags with description, use sktime.registry.all_tags
+        #   all_tags(estimator_types="transformer", as_dataframe=True)
+        #
+        #
+        # behavioural tags: transformer type
+        # ----------------------------------
+        #
+        # scitype:transform-input, scitype:transform-output, scitype:transform-labels
+        # control the input/output type of transform, in terms of scitype
+        #
+        # scitype:transform-input, scitype:transform-output should be the
+        # simplest scitype that describes the mapping, taking into account vectorization
+        # a transform that produces Series when given Series, Panel when given Panel
+        #   should have both transform-input and transform-output as "Series"
+        # a transform that produces a tabular DataFrame (Table)
+        #   when given Series or Panel should have transform-input "Series"
+        #       and transform-output as "Primitives"
         "scitype:transform-input": "Series",
-        # what is the scitype of X: Series, or Panel
+        # valid values: "Series", "Panel"
         "scitype:transform-output": "Series",
-        # what scitype is returned: Primitives, Series, Panel
+        # valid values: "Series", "Panel", "Primitives"
+        #
+        # scitype:instancewise = is fit_transform an instance-wise operation?
+        # instance-wise = only values of a given series instance are used to transform
+        #   that instance. Example: Fourier transform; non-example: series PCA
+        "scitype:instancewise": True,
+        #
+        # scitype:transform-labels types the y used in transform
+        #   if y is not used in transform, this should be "None"
         "scitype:transform-labels": "None",
-        # what is the scitype of y: None (not needed), Primitives, Series, Panel
-        "scitype:instancewise": True,  # is this an instance-wise transform?
-        "capability:inverse_transform": False,  # can the transformer inverse transform?
-        "univariate-only": False,  # can the transformer handle multivariate X?
-        "X_inner_mtype": "pd.DataFrame",  # which mtypes do _fit/_predict support for X?
-        # this can be a Panel mtype even if transform-input is Series, vectorized
-        "y_inner_mtype": "None",  # which mtypes do _fit/_predict support for y?
-        "requires_y": False,  # does y need to be passed in fit?
-        "enforce_index_type": None,  # index type that needs to be enforced in X/y
-        "fit_is_empty": True,  # is fit empty and can be skipped? Yes = True
-        "X-y-must-have-same-index": False,  # can estimator handle different X/y index?
+        # valid values: "None" (not needed), "Primitives", "Series", "Panel"
+        #
+        #
+        # behavioural tags: internal type
+        # ----------------------------------
+        #
+        # X_inner_mtype, y_inner_mtype control which format X/y appears in
+        # in the inner functions _fit, _transform, etc
+        "X_inner_mtype": "pd.DataFrame",
+        "y_inner_mtype": "None",
+        # valid values: str and list of str
+        # if str, must be a valid mtype str, in sktime.datatypes.MTYPE_REGISTER
+        #   of scitype Series, Panel (panel data) or Hierarchical (hierarchical series)
+        #   y_inner_mtype can also be of scitype Table (one row/instance per series)
+        #   in that case, all inputs are converted to that one type
+        # if list of str, must be a list of valid str specifiers
+        #   in that case, X/y are passed through without conversion if on the list
+        #   if not on the list, converted to the first entry of the same scitype
+        #
+        # univariate-only controls whether internal X can be univariate/multivariate
+        # if True (only univariate), always applies vectorization over variables
+        "univariate-only": False,
+        # valid values: True = inner _fit, _transform receive only univariate serie
+        #   False = uni- and multivariate series are passed to inner methods
+        #
+        # requires_y = does y need to be passed in fit?
+        "requires_y": False,
+        # valid values: False (no), True = exception is raised if no y is seen in _fit
+        #   y can be passed or not in _transform for either value of requires_y
+        #
+        # remember_data = whether all data seen is remembered as self._X
+        "remember_data": False,
+        # valid vales: False (no), True = self._X is created/update in fit/update
+        #   self._X is all X passed via fit or update, updated via update_data
+        #   self._X is of mtype seen in fit, update adds more data to the same container
+        #   self._X can be used (readonly) by the estimator in _fit, _transform, _update
+        #   if set to True, fit-is-empty must be set to False
+        #
+        # capability tags: properties of the estimator
+        # --------------------------------------------
+        #
+        # fit_is_empty = is fit empty and can be skipped?
+        "fit_is_empty": True,
+        # valid values: True = _fit is considered empty and skipped, False = No
+        # CAUTION: default is "True", i.e., _fit will be skipped even if implemented
+        #
+        # X-y-must-have-same-index = can estimator handle different X/y index?
+        "X-y-must-have-same-index": False,
+        # valid values: boolean True (yes), False (no)
+        # if True, raises exception if X.index is not contained in y.index
+        #
+        # enforce_index_type = index type that needs to be enforced in X/y
+        "enforce_index_type": None,
+        # valid values: pd.Index subtype, or list of pd.Index subtype
+        # if not None, raises exception if X.index, y.index level -1 is not of that type
+        #
+        # transform-returns-same-time-index = does transform return same index as input?
         "transform-returns-same-time-index": False,
-        # does transform return have the same time index as input X
-        "skip-inverse-transform": False,  # is inverse-transform skipped when called?
+        # valid values: boolean True (yes), False (no)
+        # if True, transform and inverse_transform returns should have
+        #   same length and same index (if pandas) as inputs
+        # no exception is raised if this tag is incorrectly set
+        #
+        # capability:inverse_transform = is inverse_transform implemented?
+        "capability:inverse_transform": False,
+        # valid values: boolean True (yes), False (no)
+        # if True, _inverse_transform must be implemented
+        # if False, exception is raised if inverse_transform is called,
+        #   unless the skip-inverse-transform tag is set to True
+        #
+        # skip-inverse-transform = is inverse-transform skipped when called?
+        "skip-inverse-transform": False,
+        # if False, capability:inverse_transform tag behaviour is as per devault
+        # if True, inverse_transform is the identity transform and raises no exception
+        #   this is useful for transformers where inverse_transform
+        #   may be called but should behave as the identity, e.g., imputers
+        #
+        # capability:unequal_length = can the transformer handle unequal length panels,
+        #   i.e., when passed unequal length instances in Panel or Hierarchical data
         "capability:unequal_length": True,
-        # can the transformer handle unequal length time series (if passed Panel)?
+        # valid values: boolean True (yes), False (no)
+        # if False, may raise exception when passed unequal length Panel/Hierarchical
+        #
+        # capability:unequal_length:removes = if passed Panel/Hierarchical,
+        #   is transform result always guaranteed to be equal length (and series)?
         "capability:unequal_length:removes": False,
-        # is transform result always guaranteed to be equal length (and series)?
-        #   not relevant for transformers that return Primitives in transform-output
+        # valid values: boolean True (yes), False (no)
+        # applicable only if scitype:transform-output is not "Primitives"
+        # used for search index and validity checking, does not raise direct exception
+        #
+        # handles-missing-data = can the transformer handle missing data (np or pd.NA)?
         "handles-missing-data": False,  # can estimator handle missing data?
-        # todo: rename to capability:missing_values
+        # valid values: boolean True (yes), False (no)
+        # if False, may raise exception when passed time series with missing values
+        #
+        # capability:missing_values:removes = if passed time series
+        #   is transform result always guaranteed to contain no missing values?
         "capability:missing_values:removes": False,
-        # is transform result always guaranteed to contain no missing values?
-        "python_version": None,  # PEP 440 python version specifier to limit versions
+        # valid values: boolean True (yes), False (no)
+        # used for search index and validity checking, does not raise direct exception
+        #
+        #
+        # dependency tags: python version and soft dependencies
+        # -----------------------------------------------------
+        #
+        # python version requirement
+        "python_version": None,
+        # valid values: str, PEP 440 valid python version specifiers
+        # raises exception at construction if local python veresion is incompatible
+        #
+        # soft dependency requirement
+        "python_dependencies": None
+        # valid values: str or list of str
+        # raises exception at construction if modules at strings cannot be imported
     }
     # in case of inheritance, concrete class should typically set tags
     #  alternatively, descendants can set tags in __init__
diff --git a/extension_templates/transformer_simple.py b/extension_templates/transformer_simple.py
index d79cf511f..b52c03c64 100644
--- a/extension_templates/transformer_simple.py
+++ b/extension_templates/transformer_simple.py
@@ -87,33 +87,92 @@ class MyTransformer(BaseTransformer):
     #   y_inner_mtype must be changed to one or a list of compatible sktime mtypes
     #  the other tags are "safe defaults" which can usually be left as-is
     _tags = {
+        # to list all valid tags with description, use sktime.registry.all_tags
+        #   all_tags(estimator_types="transformer", as_dataframe=True)
+        #
+        #
+        # behavioural tags: transformer type
+        # ----------------------------------
+        #
+        # scitype:transform-input, scitype:transform-output, scitype:transform-labels
+        # control the input/output type of transform, in terms of scitype
+        #
+        # scitype:transform-input, scitype:transform-output should be the
+        # simplest scitype that describes the mapping, taking into account vectorization
+        # a transform that produces Series when given Series, Panel when given Panel
+        #   should have both transform-input and transform-output as "Series"
+        # a transform that produces a tabular DataFrame (Table)
+        #   when given Series or Panel should have transform-input "Series"
+        #       and transform-output as "Primitives"
         "scitype:transform-input": "Series",
-        # what is the scitype of X: Series, or Panel
+        # valid values: "Series", "Panel"
         "scitype:transform-output": "Series",
-        # what scitype is returned: Primitives, Series, Panel
+        # valid values: "Series", "Panel", "Primitives"
+        #
+        # scitype:instancewise = is fit_transform an instance-wise operation?
+        # instance-wise = only values of a given series instance are used to transform
+        #   that instance. Example: Fourier transform; non-example: series PCA
+        "scitype:instancewise": True,
+        #
+        # scitype:transform-labels types the y used in transform
+        #   if y is not used in transform, this should be "None"
         "scitype:transform-labels": "None",
-        # what is the scitype of y: None (not needed), Primitives, Series, Panel
-        "scitype:instancewise": True,  # is this an instance-wise transform?
-        "capability:inverse_transform": False,  # can the transformer inverse transform?
-        "univariate-only": False,  # can the transformer handle multivariate X?
-        "X_inner_mtype": "pd.DataFrame",  # which mtypes do _fit/_predict support for X?
-        # this can be a Panel mtype even if transform-input is Series, vectorized
-        "y_inner_mtype": "None",  # which mtypes do _fit/_predict support for y?
-        "requires_y": False,  # does y need to be passed in fit?
-        "enforce_index_type": None,  # index type that needs to be enforced in X/y
-        "fit_is_empty": True,  # is fit empty and can be skipped? Yes = True
-        "X-y-must-have-same-index": False,  # can estimator handle different X/y index?
-        "transform-returns-same-time-index": False,
-        # does transform return have the same time index as input X
-        "skip-inverse-transform": False,  # is inverse-transform skipped when called?
+        # valid values: "None" (not needed), "Primitives", "Series", "Panel"
+        #
+        #
+        # behavioural tags: internal type
+        # ----------------------------------
+        #
+        # X_inner_mtype, y_inner_mtype control which format X/y appears in
+        # in the inner functions _fit, _transform, etc
+        "X_inner_mtype": "pd.DataFrame",
+        "y_inner_mtype": "None",
+        # valid values: str and list of str
+        # if str, must be a valid mtype str, in sktime.datatypes.MTYPE_REGISTER
+        #   of scitype Series, Panel (panel data) or Hierarchical (hierarchical series)
+        #   y_inner_mtype can also be of scitype Table (one row/instance per series)
+        #   in that case, all inputs are converted to that one type
+        # if list of str, must be a list of valid str specifiers
+        #   in that case, X/y are passed through without conversion if on the list
+        #   if not on the list, converted to the first entry of the same scitype
+        #
+        # univariate-only controls whether internal X can be univariate/multivariate
+        # if True (only univariate), always applies vectorization over variables
+        "univariate-only": False,
+        # valid values: True = inner _fit, _transform receive only univariate serie
+        #   False = uni- and multivariate series are passed to inner methods
+        #
+        # requires_y = does y need to be passed in fit?
+        "requires_y": False,
+        # valid values: False (no), True = exception is raised if no y is seen in _fit
+        #   y can be passed or not in _transform for either value of requires_y
+        #
+        #
+        # capability tags: properties of the estimator
+        # --------------------------------------------
+        #
+        # fit_is_empty = is fit empty and can be skipped?
+        "fit_is_empty": True,
+        # valid values: True = _fit is considered empty and skipped, False = No
+        # CAUTION: default is "True", i.e., _fit will be skipped even if implemented
+        #
+        # capability:inverse_transform = is inverse_transform implemented?
+        "capability:inverse_transform": False,
+        # valid values: boolean True (yes), False (no)
+        # if True, _inverse_transform must be implemented
+        # if False, exception is raised if inverse_transform is called,
+        #   unless the skip-inverse-transform tag is set to True
+        #
+        # capability:unequal_length = can the transformer handle unequal length panels,
+        #   i.e., when passed unequal length instances in Panel or Hierarchical data
         "capability:unequal_length": True,
-        # can the transformer handle unequal length time series (if passed Panel)?
-        "capability:unequal_length:removes": False,
-        # is transform result always guaranteed to be equal length (and series)?
-        #   not relevant for transformers that return Primitives in transform-output
+        # valid values: boolean True (yes), False (no)
+        # if False, may raise exception when passed unequal length Panel/Hierarchical
+        #
+        # handles-missing-data = can the transformer handle missing data (np or pd.NA)?
         "handles-missing-data": False,  # can estimator handle missing data?
-        "capability:missing_values:removes": False,
-        # is transform result always guaranteed to contain no missing values?
+        # valid values: boolean True (yes), False (no)
+        # if False, may raise exception when passed time series with missing values
     }
 
     # todo: add any hyper-parameters and components to constructor
diff --git a/pyproject.toml b/pyproject.toml
index 6b9509db6..af7f779e0 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -59,7 +59,7 @@ all_extras = [
     "gluonts>=0.9.0",
     "matplotlib>=3.3.2",
     "mne",
-    "pmdarima>=1.8.0,!=1.8.1,<2.0.0",
+    "pmdarima>=1.8.0,!=1.8.1,<3.0.0",
     "prophet>=1.1",
     "pykalman>=0.9.5",
     "pyod>=0.8.0",
diff --git a/sktime/annotation/hmm.py b/sktime/annotation/hmm.py
index 141357201..4f7b77d13 100644
--- a/sktime/annotation/hmm.py
+++ b/sktime/annotation/hmm.py
@@ -6,6 +6,7 @@ Implements a basic Hidden Markov Model (HMM) as an annotation estimator.
 To read more about the algorithm, check out the `HMM wikipedia page
 <https://en.wikipedia.org/wiki/Hidden_Markov_model>`_.
 """
+import warnings
 from typing import Tuple
 
 import numpy as np
@@ -259,7 +260,7 @@ class HMM(BaseSeriesAnnotator):
             trans_prob[:, i] = np.max(paths, axis=0)
 
         if np.any(np.isinf(trans_prob[:, -1])):
-            raise ValueError("Change parameters, the distribution doesn't work")
+            warnings.warn("Change parameters, the distribution doesn't work")
 
         return trans_prob, trans_id
 
@@ -409,7 +410,7 @@ class HMM(BaseSeriesAnnotator):
         params : dict or list of dict
         """
         centers = [3.5, -5]
-        sd = [0.25 for _ in centers]
+        sd = [100 for _ in centers]
         emi_funcs = [
             (norm.pdf, {"loc": mean, "scale": sd[ind]})
             for ind, mean in enumerate(centers)
diff --git a/sktime/annotation/stray.py b/sktime/annotation/stray.py
new file mode 100644
index 000000000..5d5d83d14
--- /dev/null
+++ b/sktime/annotation/stray.py
@@ -0,0 +1,266 @@
+# -*- coding: utf-8 -*-
+"""Tests for STRAY (Search TRace AnomalY) outlier estimator."""
+
+import warnings
+from typing import Dict
+
+import numpy as np
+import numpy.typing as npt
+from sklearn.neighbors import NearestNeighbors
+
+from sktime.transformations.base import BaseTransformer
+
+__author__ = ["KatieBuc"]
+__all__ = ["STRAY"]
+
+
+class STRAY(BaseTransformer):
+    """STRAY: robust anomaly detection in data streams with concept drift.
+
+    This is based on STRAY (Search TRace AnomalY) _[1], which is a modification
+    of HDoutliers _[2]. HDoutliers is a powerful algorithm for the detection of
+    anomalous observations in a dataset, which has (among other advantages) the
+    ability to detect clusters of outliers in multi-dimensional data without
+    requiring a model of the typical behavior of the system. However, it suffers
+    from some limitations that affect its accuracy. STRAY is an extension of
+    HDoutliers that uses extreme value theory for the anomolous threshold
+    calculation, to deal with data streams that exhibit non-stationary behavior.
+
+    Parameters
+    ----------
+    alpha : float, optional (default=0.01)
+        Threshold for determining the cutoff for outliers. Observations are
+        considered outliers if they fall in the (1 - alpha) tail of
+        the distribution of the nearest-neighbor distances between exemplars.
+    k : int, optional (default=10)
+        Number of neighbours considered.
+    knn_algorithm : str {"auto", "ball_tree", "kd_tree", "brute"}, optional
+        (default="brute")
+        Algorithm used to compute the nearest neighbors, from
+        sklearn.neighbors.NearestNeighbors
+    p : float, optional (default=0.5)
+        Proportion of possible candidates for outliers. This defines the starting point
+        for the bottom up searching algorithm.
+    size_threshold : int, optional (default=50)
+        Sample size to calculate an emperical threshold.
+    outlier_tail : str {"min", "max"}, optional (default="max")
+        Direction of the outlier tail.
+
+    Attributes
+    ----------
+    score_ : pd.Series
+        Outlier score of each data point in X.
+    y_ : pd.Series
+        Outlier boolean flag for each data point in X.
+
+    References
+    ----------
+    .. [1] Talagala, Priyanga Dilini, Rob J. Hyndman, and Kate Smith-Miles.
+    "Anomaly detection in high-dimensional data." Journal of Computational
+    and Graphical Statistics 30.2 (2021): 360-374.
+    .. [2] Wilkinson, Leland. "Visualizing big data outliers through
+    distributed aggregation." IEEE transactions on visualization and
+    computer graphics 24.1 (2017): 256-266.
+
+    Examples
+    --------
+    >>> from sktime.annotation.stray import STRAY
+    >>> from sktime.datasets import load_airline
+    >>> from sklearn.preprocessing import MinMaxScaler
+    >>> import numpy as np
+    >>> X = load_airline().to_frame().to_numpy()
+    >>> scaler = MinMaxScaler()
+    >>> X = scaler.fit_transform(X)
+    >>> model = STRAY(k=3)
+    >>> y = model.fit_transform(X)
+    >>> y[:5]
+    array([False, False, False, False, False])
+    """
+
+    _tags = {
+        "handles-missing-data": True,
+        "X_inner_mtype": "np.ndarray",
+        "fit_is_empty": False,
+        "skip-inverse-transform": True,
+    }
+
+    def __init__(
+        self,
+        alpha: float = 0.01,
+        k: int = 10,
+        knn_algorithm: str = "brute",
+        p: float = 0.5,
+        size_threshold: int = 50,
+        outlier_tail: str = "max",
+    ):
+        self.alpha = alpha
+        self.k = k
+        self.knn_algorithm = knn_algorithm
+        self.p = p
+        self.size_threshold = size_threshold
+        self.outlier_tail = outlier_tail
+        super(STRAY, self).__init__()
+
+    def _find_threshold(self, outlier_score: npt.ArrayLike, n: int) -> npt.ArrayLike:
+        """Find Outlier Threshold.
+
+        Parameters
+        ----------
+        outlier_score : np.ArrayLike
+            The outlier scores determined by k nearast neighbours distance
+        n : int
+            The number of rows remaining in X when NA's are removed.
+
+        Returns
+        -------
+        array of indices of the observations determined to be outliers.
+        """
+        if self.outlier_tail == "min":
+            outlier_score = -outlier_score
+
+        order = np.argsort(outlier_score)
+        gaps = np.append(0, np.diff(outlier_score[order]))
+        n4 = int(max(min(self.size_threshold, np.floor(n / 4)), 2))
+
+        J = np.array(range(2, n4 + 1))
+        start = int(max(np.floor(n * (1 - self.p)), 1))
+
+        ghat = [
+            0.0 if i < start else sum((J / (n4 - 1)) * gaps[i - J + 1])
+            for i in range(n)
+        ]
+
+        log_alpha = np.log(1 / self.alpha)
+        bound = np.Inf
+
+        for i in range(start, n):
+            if gaps[i] > log_alpha * ghat[i]:
+                bound = outlier_score[order][i - 1]
+                break
+
+        return np.where(outlier_score > bound)[0]
+
+    def _find_outliers_kNN(self, X: npt.ArrayLike, n: int) -> Dict:
+        """Find outliers using kNN distance with maximum gap.
+
+        Parameters
+        ----------
+        X : np.ArrayLike
+            Data for anomaly detection (time series).
+        n : int
+            The number of rows remaining in X when NA's are removed.
+
+        Returns
+        -------
+        dict of index of outliers and the outlier scores
+        """
+        if len(X.shape) == 1:
+            nbrs = NearestNeighbors(n_neighbors=self.k + 1).fit(X.reshape(-1, 1))
+            distances, _ = nbrs.kneighbors(X.reshape(-1, 1))
+        else:
+            nbrs = NearestNeighbors(
+                n_neighbors=n if self.k >= n else self.k + 1,
+                algorithm=self.knn_algorithm,
+            ).fit(X)
+            distances, _ = nbrs.kneighbors(X)
+
+        if self.k == 1:
+            d = distances[:, 1]
+        else:
+            diff = np.apply_along_axis(np.diff, 1, distances)
+            d = distances[range(n), np.apply_along_axis(np.argmax, 1, diff) + 1]
+
+        out_index = self._find_threshold(d, n)
+        return {"idx_outliers": out_index, "out_scores": d}
+
+    def _find_outliers(self, X):
+        """Detect Anomalies in High Dimensional Data.
+
+        Parameters
+        ----------
+        X : np.ArrayLike
+            Data for anomaly detection (time series).
+
+        Returns
+        -------
+        dict of anomalies and their corresponding scores
+        """
+        r = np.shape(X)[0]
+        idx_dropna = np.array([i for i in range(r) if not np.isnan(X[i]).any()])
+        X_dropna = X[
+            idx_dropna,
+        ]
+
+        n = np.shape(X_dropna)[0]
+        outliers = self._find_outliers_kNN(X_dropna, n)
+
+        # adjusted back to length r, for missing data
+        slice_ = [True if i in outliers["idx_outliers"] else False for i in range(n)]
+        idx_outliers = idx_dropna[slice_]  # index values from 1:r
+        outlier_bool = np.array([1 if i in idx_outliers else 0 for i in range(r)])
+
+        list_scores = outliers["out_scores"].tolist()
+        outlier_scores = np.array(
+            [list_scores.pop(0) if i in idx_dropna else np.nan for i in range(r)]
+        )
+
+        return {
+            "outlier_scores": outlier_scores,
+            "outlier_bool": outlier_bool,
+        }
+
+    def _fit(self, X: npt.ArrayLike, y: npt.ArrayLike = None):
+        """Find outliers using STRAY (Search TRace AnomalY).
+
+        Parameters
+        ----------
+        X : np.ArrayLike
+            Data for anomaly detection (time series).
+        y : pd.Series, optional
+            Not used for this unsupervsed method.
+
+        Returns
+        -------
+        self :
+            Reference to self.
+        """
+        # remember X for transform
+        self._X = X
+
+        info_dict = self._find_outliers(X)
+        self.score_ = info_dict["outlier_scores"]
+        self.y_ = info_dict["outlier_bool"]
+
+        return self
+
+    def _transform(self, X: npt.ArrayLike, y: npt.ArrayLike = None) -> npt.ArrayLike:
+        """Return anomaly detection.
+
+        Parameters
+        ----------
+        X : np.ArrayLike
+            Data for anomaly detection (time series).
+
+        Returns
+        -------
+        y_ : np.ArrayLike
+            Anomaly detection, boolean.
+        """
+        # fit again if data is different to fit, but don't store anything
+        if not np.allclose(X, self._X, equal_nan=True):
+            new_obj = STRAY(
+                alpha=self.alpha,
+                k=self.k,
+                knn_algorithm=self.knn_algorithm,
+                p=self.p,
+                size_threshold=self.size_threshold,
+                outlier_tail=self.outlier_tail,
+            ).fit(X)
+            warnings.warn(
+                "Warning: Input data X differs from that given to fit(). "
+                "Refitting with new input data, not storing updated public class "
+                "attributes. For this, explicitly use fit(X) or fit_transform(X)."
+            )
+            return new_obj.y_.astype(bool)
+
+        return self.y_.astype(bool)
diff --git a/sktime/base/_base.py b/sktime/base/_base.py
index 09fb7eead..14f2ca9cd 100644
--- a/sktime/base/_base.py
+++ b/sktime/base/_base.py
@@ -816,7 +816,7 @@ class BaseEstimator(BaseObject):
         """
         if not self.is_fitted:
             raise NotFittedError(
-                f"parameter estimator of type {type(self).__name__} has not been "
+                f"estimator of type {type(self).__name__} has not been "
                 "fitted yet, please call fit on data before get_fitted_params"
             )
 
@@ -831,9 +831,11 @@ class BaseEstimator(BaseObject):
                 return x
 
         for c in c_dict.keys():
-            c_f_params = c_dict[c].get_fitted_params()
-            c_f_params = {f"{sh(c)}__{k}": c_f_params[k] for k in c_f_params.keys()}
-            fitted_params.update(c_f_params)
+            comp = c_dict[c]
+            if comp._is_fitted:
+                c_f_params = c_dict[c].get_fitted_params()
+                c_f_params = {f"{sh(c)}__{k}": c_f_params[k] for k in c_f_params.keys()}
+                fitted_params.update(c_f_params)
 
         fitted_params.update(self._get_fitted_params())
 
diff --git a/sktime/classification/_delegate.py b/sktime/classification/_delegate.py
index cef272bd0..e5873c526 100644
--- a/sktime/classification/_delegate.py
+++ b/sktime/classification/_delegate.py
@@ -122,3 +122,18 @@ class _DelegatedClassifier(BaseClassifier):
         """
         estimator = self._get_delegate()
         return estimator.predict_proba(X=X)
+
+    def _get_fitted_params(self):
+        """Get fitted parameters.
+
+        private _get_fitted_params, called from get_fitted_params
+
+        State required:
+            Requires state to be "fitted".
+
+        Returns
+        -------
+        fitted_params : dict
+        """
+        estimator = self._get_delegate()
+        return estimator.get_fitted_params()
diff --git a/sktime/datatypes/__init__.py b/sktime/datatypes/__init__.py
index 67b1c981e..8e5e22e48 100644
--- a/sktime/datatypes/__init__.py
+++ b/sktime/datatypes/__init__.py
@@ -13,6 +13,7 @@ from sktime.datatypes._check import (
 from sktime.datatypes._convert import convert, convert_to
 from sktime.datatypes._examples import get_examples
 from sktime.datatypes._registry import (
+    ALL_TIME_SERIES_MTYPES,
     MTYPE_LIST_HIERARCHICAL,
     MTYPE_LIST_PANEL,
     MTYPE_LIST_PROBA,
@@ -28,6 +29,7 @@ from sktime.datatypes._utilities import get_cutoff, update_data
 from sktime.datatypes._vectorize import VectorizedDF
 
 __all__ = [
+    "ALL_TIME_SERIES_MTYPES",
     "check_is_mtype",
     "check_is_scitype",
     "check_raise",
diff --git a/sktime/datatypes/_check.py b/sktime/datatypes/_check.py
index 2e2d9fb8d..d26d49958 100644
--- a/sktime/datatypes/_check.py
+++ b/sktime/datatypes/_check.py
@@ -32,7 +32,7 @@ from sktime.datatypes._alignment import check_dict_Alignment
 from sktime.datatypes._hierarchical import check_dict_Hierarchical
 from sktime.datatypes._panel import check_dict_Panel
 from sktime.datatypes._proba import check_dict_Proba
-from sktime.datatypes._registry import SCITYPE_LIST, mtype_to_scitype
+from sktime.datatypes._registry import AMBIGUOUS_MTYPES, SCITYPE_LIST, mtype_to_scitype
 from sktime.datatypes._series import check_dict_Series
 from sktime.datatypes._table import check_dict_Table
 
@@ -46,10 +46,6 @@ check_dict.update(check_dict_Table)
 check_dict.update(check_dict_Proba)
 
 
-# mtypes to exclude
-AMBIGUOUS_MTYPES = ["numpyflat", "alignment_loc"]
-
-
 def _check_scitype_valid(scitype: str = None):
     """Check validity of scitype."""
     valid_scitypes = list(set([x[1] for x in check_dict.keys()]))
diff --git a/sktime/datatypes/_registry.py b/sktime/datatypes/_registry.py
index bdf475873..910173d11 100644
--- a/sktime/datatypes/_registry.py
+++ b/sktime/datatypes/_registry.py
@@ -53,6 +53,15 @@ MTYPE_REGISTER += MTYPE_REGISTER_ALIGNMENT
 MTYPE_REGISTER += MTYPE_REGISTER_TABLE
 MTYPE_REGISTER += MTYPE_REGISTER_PROBA
 
+# mtypes to exclude in checking since they are ambiguous and rare
+AMBIGUOUS_MTYPES = ["numpyflat", "alignment_loc"]
+
+# all time series mtypes excluding ambiguous ones
+ALL_TIME_SERIES_MTYPES = (
+    list(MTYPE_LIST_PANEL) + list(MTYPE_LIST_SERIES) + list(MTYPE_LIST_HIERARCHICAL)
+)
+ALL_TIME_SERIES_MTYPES = list(set(ALL_TIME_SERIES_MTYPES).difference(AMBIGUOUS_MTYPES))
+
 
 __all__ = [
     "MTYPE_REGISTER",
diff --git a/sktime/dists_kernels/signature_kernel.py b/sktime/dists_kernels/signature_kernel.py
new file mode 100644
index 000000000..e4653a09a
--- /dev/null
+++ b/sktime/dists_kernels/signature_kernel.py
@@ -0,0 +1,1006 @@
+# -*- coding: utf-8 -*-
+"""Signature kernels from Kiraly et al, 2016."""
+
+__author__ = ["fkiraly"]
+
+import collections
+from functools import partial
+
+import numpy as np
+from scipy.sparse.linalg import svds
+from sklearn.base import BaseEstimator, TransformerMixin
+
+from sktime.dists_kernels._base import BasePairwiseTransformerPanel
+
+# cumsum varia
+# ------------
+
+
+def _coerce_to_list_or_tuple(i):
+    """Coerce integers to list of integers."""
+    if not isinstance(i, (list, tuple)):
+        return [i]
+    else:
+        return i
+
+
+def cumsum_rev_first(array):
+    """Reverse cumsum over 0th axis."""
+    out = np.zeros_like(array)
+    out[:-1] = np.cumsum(array[:0:-1], axis=0)[::-1]
+    return out
+
+
+def cumsum_rev(array):
+    """Reverse cumsum over both axes 0 and 1."""
+    out = np.zeros_like(array)
+    out[:-1, :-1] = np.cumsum(np.cumsum(array[:0:-1, :0:-1], 0), 1)[::-1, ::-1]
+    return out
+
+
+def cumsum_mult(array, dims):
+    """Cumsum over all axes in dims."""
+    dims = _coerce_to_list_or_tuple(dims)
+    for dimind in dims:
+        array = np.cumsum(array, axis=dimind)
+    return array
+
+
+def roll_mult(array, shift, dims):
+    """Roll over all axes in dims."""
+    dims = _coerce_to_list_or_tuple(dims)
+    for dimind in dims:
+        array = np.roll(array, shift, axis=dimind)
+    return array
+
+
+def makeinds(indlist):
+    """Make mesh from sequences."""
+    return np.ix_(*indlist)
+
+
+def cumsum_shift_mult(array, dims):
+    """Apply cumsum and shift to all axes in dims."""
+    array = cumsum_mult(array, dims)
+    array = roll_mult(array, 1, dims)
+
+    arrayshape = array.shape
+    indarr = []
+    for ind in range(len(arrayshape)):
+        indarr = indarr + [range(arrayshape[ind])]
+
+    for dimind in dims:
+        slicearr = indarr[:]
+        slicearr[dimind] = [0]
+        array[makeinds(slicearr)] = 0
+
+    return array
+
+
+# low rank reduction utilities
+# ----------------------------
+
+
+def rankreduce(array, rankbound):
+    """Project 2D array on top rankbound singular values."""
+    arraysvd = svds(array.astype("f"), k=rankbound)
+    return np.dot(arraysvd[0], np.diag(arraysvd[1]))
+
+
+def rankreduce_batch(arrays, rankbound):
+    """Apply rankreduce to axis 0 stack of arrays."""
+    resultarrays = np.zeros([arrays.shape[0], arrays.shape[1], rankbound])
+    for i in range(arrays.shape[0]):
+        resultarrays[i, :, :] = rankreduce(arrays[i, :, :], rankbound)
+    return resultarrays
+
+
+# kernels and distances
+# ---------------------
+
+
+def sqdist(X, Y):
+    """Row-wise squared distance between 2D array X and 2D array Y."""
+    M = np.shape(X)[0]
+    N = np.shape(Y)[0]
+    X_squared = np.tile((X * X).sum(-1), [N, 1]).T
+    Y_squared = np.tile((Y * Y).sum(-1), [M, 1])
+    X_times_Y = np.inner(X, Y)
+    return X_squared + Y_squared - 2 * X_times_Y
+
+
+def k_polynom(x, y, scale, deg):
+    """Polynomial kernel of degree deg, with scale coeff."""
+    return (1 + scale * np.inner(x, y)) ** deg
+
+
+def k_gauss(x, y, scale):
+    """Gaussian kernel with scale coeff."""
+    return np.exp(-(scale**2) * sqdist(x, y) / 2)
+
+
+def k_euclid(x, y, scale):
+    """Euclidean kernel with scale coeff."""
+    return scale * np.inner(x, y)
+
+
+def k_laplace(x, y, scale):
+    """Laplace kernel with scale coeff."""
+    return np.exp(-scale * np.sqrt(np.inner(x - y, x - y)))
+
+
+def k_tanh(x, y, off, scale):
+    """Tanh kernel with scale and offset."""
+    return np.tanh(off + scale * np.inner(x, y))
+
+
+def mirror(K):
+    """Mirrors an upper triangular kernel matrix, helper for Sqize_kernel."""
+    return K - np.diag(np.diag(K)) + np.transpose(K)
+
+
+# sequential kernel implementation
+# --------------------------------
+
+
+def sqize_kernel(K, L, theta=1.0, normalize=False):
+    """Compute the sequential kernel from a pairwise kernel matrix.
+
+    aka compute "sequentialization of the kernel K"
+
+    Parameters
+    ----------
+    K : 2D np.ndarray
+        the kernel matrix of increments, i.e.,
+        K[i,j] is the kernel between the i-th increment of path 1 (rows),
+        and the j-th increment of path 2 (columns)
+    L : an integer >= 1, representing the level of truncation
+    theta : float, optional, default=1.0
+        a positive scaling factor for the levels, i-th level is scaled by theta^i
+    normalize : bool, optional, default = False
+        whether the output kernel matrix is normalized
+        if True, sums and cumsums are divided by prod(K.shape)
+
+    Returns
+    -------
+    a real number, the sequential kernel between path 1 (rows of K) and path 2 (cols)
+    """
+    # L-1 runs through loop;
+    # returns R_ij=(1+\sum_i2>i,j2>j A_i2,j2(1+\sum A_iLjL)...)
+    if normalize:
+        normfac = np.prod(K.shape)
+        Id = np.ones(K.shape)
+        R = np.ones(K.shape)
+        for _ in range(L - 1):
+            R = (Id + theta * cumsum_rev(K * R) / normfac) / (1 + theta)
+        return (1 + theta * np.sum(K * R) / normfac) / (1 + theta)
+    else:
+        Id = np.ones(K.shape)
+        R = np.ones(K.shape)
+        for _ in range(L - 1):
+            R = Id + cumsum_rev(K * R)  # A*R is componentwise
+        # outermost bracket: since i1>=1 and not i1>1 we do it outside of loop
+        return 1 + np.sum(K * R)
+
+
+def sqize_kernel_ho(K, L, D=1, theta=1.0, normalize=False):
+    """Compute the higher-order sequential kernel from a pairwise kernel matrix.
+
+    aka compute "sequentialization of the kernel K", higher-order approximation
+
+    Parameters
+    ----------
+    K : 2D np.ndarray
+        the kernel matrix of increments, i.e.,
+        K[i,j] is the kernel between the i-th increment of path 1 (rows),
+        and the j-th increment of path 2 (columns)
+    L : an integer >= 1, representing the level of truncation
+    D : int, optional, default = 1
+        an integer >= 1, representing the order of approximation
+    theta : float, optional, default=1.0
+        a positive scaling factor for the levels, i-th level is scaled by theta^i
+    normalize : bool, optional, default = False
+        whether the output kernel matrix is normalized
+        if True, sums and cumsums are divided by prod(K.shape)
+
+    Returns
+    -------
+    a real number, the sequential kernel between path 1 (rows of K) and path 2 (cols)
+    """
+    A = np.zeros(np.concatenate(([L, D, D], K.shape)))
+    Id = np.ones(K.shape)
+
+    for ell in range(1, L):
+        Dprime = min(D, ell)
+        Acs = cumsum_shift_mult(np.sum(A[ell - 1, :, :, :, :], (0, 1)), (0, 1))
+        A[ell, 0, 0, :, :] = K * (Id + Acs)
+        for d1 in range(1, Dprime):
+            Acs1 = cumsum_shift_mult(np.sum(A[ell - 1, d1 - 1, :, :, :], 0), 1)
+            Acs2 = cumsum_shift_mult(np.sum(A[ell - 1, :, d1 - 1, :, :], 0), 0)
+            A[ell, d1, 0, :, :] = A[ell, d1, 0, :, :] + (1 / d1) * K * Acs1
+            A[ell, :, d1, :, :] = A[ell, 0, d1, :, :] + (1 / d1) * K * Acs2
+
+            for d2 in range(1, Dprime):
+                Acs12 = cumsum_shift_mult(
+                    np.sum(A[ell - 1, d1 - 1, d2 - 1, :, :], 0), 0
+                )
+                A[ell, d1, d2, :, :] = (
+                    A[ell, d1, d2, :, :] + (1 / (d1 * d2)) * K * Acs12
+                )
+    return 1 + np.sum(A[L - 1, :, :, :, :])
+
+
+# low-rank decomposition
+# ----------------------
+
+# LRdec object
+#  models matrix A = U x V.T
+#  U and V should be *arrays*, not *matrices*
+LRdec = collections.namedtuple("LRdec", ["U", "V"])
+
+
+def get_low_rank_matrix(K):
+    """Produce the matrix from the LRdec object.
+
+    Parameters
+    ----------
+    K : LRdec type object
+
+    Returns
+    -------
+    2D np.ndarray the matrix K.U x K.V.T modelled by the LRdec object
+    """
+    return np.inner(K.U, K.V)
+
+
+def add_low_rank(K, R):
+    """Efficient addition of two low-rank matrices.
+
+    efficient computation of sum of low-rank representations
+    using this and then get_low_rank_matrix is more efficient than an
+    explicit computation if the rank of the final matrix is not full
+
+    Parameters
+    ----------
+    K, R : LRdec type objects to add
+
+    Returns
+    -------
+    LRdec type object for sum of K and R
+    """
+    return LRdec(np.concatenate((K.U, R.U), axis=1), np.concatenate((K.V, R.V), axis=1))
+
+
+def add_low_rank_one(U, P):
+    """Efficient addition of low-rank matrices (symmetric)."""
+    return np.concatenate((U, P), axis=1)
+
+
+def mult_low_rank(K, theta):
+    """Efficient multiplication of two low-rank matrices.
+
+    efficient multiplication of sum of low-rank representations
+    using this and then get_low_rank_matrix is more efficient than an
+    explicit computation if the rank of the final matrix is not full
+
+    Parameters
+    ----------
+    K, R : LRdec type objects to add
+
+    Returns
+    -------
+    LRdec type object for product of K and R
+    """
+    return LRdec(theta * K.U, theta * K.V)
+
+
+def hadamard_low_rank(K, R):
+    """Efficient Hadamard product of two low-rank matrices.
+
+    efficient multiplication of sum of low-rank representations
+    using this and then get_low_rank_matrix is more efficient than an
+    explicit computation if the rank of the final matrix is not full
+
+    Parameters
+    ----------
+    K, R : LRdec type objects to add
+
+    Returns
+    -------
+    LRdec type object for Hadamard product of K and R
+    """
+    rankK = K.U.shape[1]
+    rankR = R.U.shape[1]
+    U = np.tile(K.U, rankR) * np.repeat(R.U, rankK, 1)
+    V = np.tile(K.V, rankR) * np.repeat(R.V, rankK, 1)
+    return LRdec(U, V)
+
+
+# multiplies U with every component (1st index) of P
+# def HadamardLowRankBatch(U, P):
+#    rankU = U.shape[1]
+#    N = P.shape[0]
+#    rankP = P.shape[2]
+#    return (np.repeat(np.repeat(np.array(U,ndmin = 3), rankP, 2),N,0)
+#           *np.repeat(P,rankU,2))
+
+
+def hadamard_low_rank_batch(U, P):
+    """Hadamard multiply U and P component-wise (1st)."""
+    rankU = U.shape[2]
+    rankP = P.shape[2]
+    return np.tile(U, rankP) * np.repeat(P, rankU, 2)
+
+
+def hadamard_low_rank_subsample(U, P, rho):
+    """Hadamard multiply U and P component-wise (1st), with Nyström type subsampling."""
+    rankU = U.shape[2]
+    rankP = P.shape[2]
+    permut = np.sort(np.random.permutation(range(rankU * rankP))[range(rho)])
+    return (np.tile(U, rankP) * np.repeat(P, rankU, 2))[:, :, permut]
+
+
+def cumsum_rev_low_rank(K):
+    """Compute cumsum for LRdec type collections.
+
+    equivalent of cumsum_rev for LRdec type objects
+
+    Parameters
+    ----------
+    K : LRdec type object to compute cumsum of
+
+    Return
+    ------
+    LRdec type object for cumsum_rev of K
+    """
+    return LRdec(cumsum_rev_first(K.U), cumsum_rev_first(K.V))
+
+
+def sum_low_rank(K):
+    """Compute cumsum for LRdec type collections.
+
+    equivalent of sum for LRdec type objects
+
+    Parameters
+    ----------
+    K : LRdec type object to compute sum of
+
+    Return
+    ------
+    LRdec type object for sum of K
+    """
+    return np.inner(sum(K.U), sum(K.V))
+
+
+# sequential kernel - low-rank version
+# ------------------------------------
+
+
+def sqize_kernel_low_rank(K, L, theta=1.0, normalize=False, rankbound=float("inf")):
+    """Compute the sequential kernel from kernel matrix, with low-rank approximation.
+
+    Parameters
+    ----------
+    K : 2D np.ndarray
+        the kernel matrix of increments, i.e.,
+        K[i,j] is the kernel between the i-th increment of path 1 (rows),
+        and the j-th increment of path 2 (columns)
+    L : an integer >= 1, representing the level of truncation
+    theta : float, optional, default=1.0
+        a positive scaling factor for the levels, i-th level is scaled by theta^i
+    normalize : bool, optional, default = False
+        whether the output kernel matrix is normalized
+        if True, sums and cumsums are divided by prod(K.shape)
+    rankbound : int, optional, default = infinity
+        a hard threshold for the rank of the level matrices
+
+    Returns
+    -------
+    a real number, the sequential kernel between path 1 (rows of K) and path 2 (cols)
+    """
+    # L-1 runs through loop;
+    # returns R_ij=(1+\sum_i2>i,j2>j A_i2,j2(1+\sum A_iLjL)...)
+    if normalize:
+        K = get_low_rank_matrix(K)
+        normfac = np.prod(K.shape)
+        Id = np.ones(K.shape)
+        R = np.ones(K.shape)
+        for _ in range(L - 1):
+            R = (Id + theta * cumsum_rev(K * R) / normfac) / (1 + theta)
+        return (1 + theta * np.sum(K * R) / normfac) / (1 + theta)
+    else:
+        Id = LRdec(np.ones([K.U.shape[0], 1]), np.ones([K.V.shape[0], 1]))
+        # Id = np.ones(K.shape)
+        R = Id
+        for _ in range(L - 1):
+            # todo: execute only if rank is lower than rankbound
+            #       reduce to rank
+            R = add_low_rank(
+                Id, mult_low_rank(cumsum_rev_low_rank(hadamard_low_rank(K, R)), theta)
+            )
+            # R = Id + cumsum_rev(K * R)
+        return 1 + theta * sum_low_rank(hadamard_low_rank(K, R))
+        # return 1 + np.sum(K * R)
+        # outermost bracket: since i1>=1 and not i1>1 we do it outside of loop
+
+
+def sqize_kernel_low_rank_fast(
+    K,
+    L,
+    theta=1.0,
+    normalize=False,
+    rankbound=float("inf"),
+):
+    """Compute the sequential kernel from kernel matrix, with low-rank approximation.
+
+    Vectorized across series, and using (faster) low-rank approximation across series.
+
+    Parameters
+    ----------
+    K : 3D np.ndarray
+        containing joint low-rank factors (symmetric singular factor)
+            1st index counts sequences
+            2nd index counts time
+            3rd index counts features
+        so K[m,:,:] is the mth factor,
+        and K[m,:,:] x K[m,:,:]^t is the kernel matrix of the mth factor
+    L : an integer >= 1, representing the level of truncation
+    theta : float, optional, default=1.0
+        a positive scaling factor for the levels, i-th level is scaled by theta^i
+    normalize : bool, optional, default = False
+        whether the output kernel matrix is normalized
+        if True, sums and cumsums are divided by prod(K.shape)
+    rankbound : int, optional, default = infinity
+        a hard threshold for the rank of the level matrices
+
+    Returns
+    -------
+    np.ndarray of shape (m, r), where r = min(rankbound, K.shape[2])
+        a matrix R such that R*R^t is the sequential kernel matrix
+        R*R^t[i,j] is sequential kernel (low-rank) between i-th and j-th sequence in K
+    """
+    if normalize:
+
+        Ksize = K.shape[0]
+        B = np.ones([Ksize, 1, 1])
+        R = np.ones([Ksize, 1])
+
+        for _ in range(L):
+            P = np.sqrt(theta) * hadamard_low_rank_batch(K, B) / Ksize
+            B = cumsum_shift_mult(P, [1])
+
+            if rankbound < B.shape[2]:
+                # B = rankreduce_batch(B, rankbound)
+                permut = np.sort(
+                    np.random.permutation(range(B.shape[2]))[range(rankbound)]
+                )
+                B = B[:, :, permut]
+
+            R = np.concatenate((R, np.sum(B, axis=1)), axis=1) / (np.sqrt(1 + theta))
+
+        return R
+
+    else:
+        Ksize = K.shape[0]
+        B = np.ones([Ksize, 1, 1])
+        R = np.ones([Ksize, 1])
+
+        for _ in range(L):
+            # todo: execute only if rank is lower than rankbound
+            #       reduce to rank
+            P = np.sqrt(theta) * hadamard_low_rank_batch(K, B)
+            B = cumsum_shift_mult(P, [1])
+
+            if rankbound < B.shape[2]:
+                # B = rankreduce_batch(B, rankbound)
+                permut = np.sort(
+                    np.random.permutation(range(B.shape[2]))[range(rankbound)]
+                )
+                B = B[:, :, permut]
+
+            R = np.concatenate((R, np.sum(B, axis=1)), axis=1)
+
+        return R
+
+
+# sequential kernel - wraps all versions
+# --------------------------------------
+
+
+def seq_kernel(
+    X,
+    kernelfun=None,
+    L=2,
+    D=1,
+    theta=1.0,
+    normalize=False,
+    lowrank=False,
+    rankbound=float("inf"),
+):
+    """Compute the sequential kernel between seqeuence/time series.
+
+    Provides interface for vanilla sequential kernel, low-rank, and higher-order.
+
+    Parameters
+    ----------
+    X : 3D np.ndarray of shape (N, d, _)
+        collection of sequences/time series
+        1st index = instance index
+        2nd index = variable/feature index
+        3rd index = time index
+    kernelfun : function (2D np.ndarray x 2D np.ndarray) -> 2D np.ndarray
+        pairwise kernel function, matrix sizes (n, d) x (m, d) -> (n x m)
+        optional, default = Euclidean (linear) kernel with scale parameter 1
+    L : int, optional, default = 2
+        an integer >= 1, representing the level of truncation
+    D : int, optional, default = 1
+        an integer >= 1, representing the order of approximation
+        can be set only if lowrank = False, otherwise ignored (always = 1)
+    theta : float, optional, default=1.0
+        a positive scaling factor for the levels, i-th level is scaled by theta^i
+    normalize : bool, optional, default = False
+        whether the output kernel matrix is normalized
+        if True, sums and cumsums are divided by prod(K.shape)
+    lowrank : bool, optional, default = False
+        whether to use low rank approximation in computing the kernel
+    rankbound : int, optional, default = infinity
+        a hard threshold for the rank of the level matrices
+        used only if lowrank = True
+
+    Returns
+    -------
+    np.ndarray of shape (N, N), sequential kernel matrix
+        [i,j]-th entry is sequential kernel between X[i] and X[j]
+    """
+    N = np.shape(X)[0]
+    KSeq = np.zeros((N, N))
+
+    if kernelfun is None:
+        kernelfun = partial(k_euclid, scale=1)
+
+    if not lowrank:
+        if D == 1:
+            for row1ind in range(N):
+                for row2ind in range(row1ind + 1):
+                    KSeq[row1ind, row2ind] = sqize_kernel(
+                        K=kernelfun(X[row1ind].T, X[row2ind].T),
+                        L=L,
+                        theta=theta,
+                        normalize=normalize,
+                    )
+        else:
+            for row1ind in range(N):
+                for row2ind in range(row1ind + 1):
+                    KSeq[row1ind, row2ind] = sqize_kernel_ho(
+                        K=kernelfun(X[row1ind].T, X[row2ind].T),
+                        L=L,
+                        D=D,
+                        theta=theta,
+                        normalize=normalize,
+                    )
+    else:
+        R = sqize_kernel_low_rank_fast(
+            K=X.transpose([0, 2, 1]),
+            L=L,
+            theta=theta,
+            normalize=normalize,
+            rankbound=rankbound,
+        )
+        KSeq = np.inner(R, R)
+        # todo: kernelfun gives back a LRdec object
+        #  for now, linear low-rank approximation is done
+        # KSeq[row1ind,row2ind] = qqize_kernel_low_rank(
+        #   kernelfun(X[row1ind].T,X[row2ind].T),L,theta,normalize = True)
+    return mirror(KSeq)
+
+
+def seq_kernel_XY(
+    X,
+    Y=None,
+    kernelfun=None,
+    L=2,
+    D=1,
+    theta=1.0,
+    normalize=False,
+    lowrank=False,
+    rankbound=float("inf"),
+):
+    """Compute the sequential kernel between two different collections of seqeuence.
+
+    Provides interface for vanilla sequential kernel, low-rank, and higher-order.
+
+    Parameters
+    ----------
+    X : 3D np.ndarray of shape (N, d, _)
+        collection of sequences/time series
+        1st index = instance index
+        2nd index = variable/feature index
+        3rd index = time index
+    Y : 3D np.ndarray of shape (M, d, _)
+        collection of sequences/time series
+        1st index = instance index
+        2nd index = variable/feature index
+        3rd index = time index
+    kernelfun : function (2D np.ndarray x 2D np.ndarray) -> 2D np.ndarray
+        pairwise kernel function, matrix sizes (n, d) x (m, d) -> (n x m)
+        optional, default = Euclidean (linear) kernel with scale parameter 1
+    L : int, optional, default = 2
+        an integer >= 1, representing the level of truncation
+    D : int, optional, default = 1
+        an integer >= 1, representing the order of approximation
+        can be set only if lowrank = False, otherwise ignored (always = 1)
+    theta : float, optional, default=1.0
+        a positive scaling factor for the levels, i-th level is scaled by theta^i
+    normalize : bool, optional, default = False
+        whether the output kernel matrix is normalized
+        if True, sums and cumsums are divided by prod(K.shape)
+    lowrank : bool, optional, default = False
+        whether to use low rank approximation in computing the kernel
+    rankbound : int, optional, default = infinity
+        a hard threshold for the rank of the level matrices
+        used only if lowrank = True
+
+    Returns
+    -------
+    np.ndarray of shape (N, M), sequential kernel matrix
+        [i,j]-th entry is sequential kernel between X[i] and Y[j]
+    """
+    # if no Y is passed, call seq_kernel
+    if Y is None:
+        return seq_kernel(
+            X=X,
+            kernelfun=kernelfun,
+            L=L,
+            D=D,
+            theta=theta,
+            normalize=normalize,
+            lowrank=lowrank,
+            rankbound=rankbound,
+        )
+
+    N = np.shape(X)[0]
+    M = np.shape(Y)[0]
+
+    KSeq = np.zeros((N, M))
+
+    if kernelfun is None:
+        kernelfun = partial(k_euclid, scale=1)
+
+    kwargs = {"L": L, "theta": theta, "normalize": normalize}
+
+    if not lowrank:
+        if D == 1:
+            for row1ind in range(N):
+                for row2ind in range(M):
+                    KSeq[row1ind, row2ind] = sqize_kernel(
+                        K=kernelfun(X[row1ind].T, Y[row2ind].T), **kwargs
+                    )
+        else:
+            for row1ind in range(N):
+                for row2ind in range(M):
+                    KSeq[row1ind, row2ind] = sqize_kernel_ho(
+                        K=kernelfun(X[row1ind].T, Y[row2ind].T), D=D, **kwargs
+                    )
+    else:
+        U = sqize_kernel_low_rank_fast(
+            K=X.transpose([0, 2, 1]), rankbound=rankbound, **kwargs
+        )
+        V = sqize_kernel_low_rank_fast(
+            K=Y.transpose([0, 2, 1]), rankbound=rankbound, **kwargs
+        )
+        KSeq = np.inner(U, V)
+        # KSeq = np.inner(sqize_kernel_low_rank_fast(X, L, theta, normalize),
+        # wqize_kernel_low_rank_fast(Y, L, theta, normalize))
+
+    return KSeq
+
+
+def data_tabulator(X):
+    """Tabulates sequence 3D np.ndarray into sklearn compatible 2D np.ndarray format."""
+    Xshape = np.shape(X)
+    return np.reshape(X, (Xshape[0], np.prod(Xshape[1:])))
+
+
+def time_series_reshaper(Xflat, numfeatures, subsample=1, differences=True):
+    """Convert 2D np.ndarray into a time series 3D np.ndarray.
+
+    Useful as part of sklearn pipeline for internal conversion to time series 3D format.
+
+    optionally, subsamples or differences time series
+
+    Parameters
+    ----------
+    Xflat : 2D np.ndarray (instances, flattened time series)
+    numfeatures : number of features/variables in the time series
+    subsample : int, optional, default = 1
+        time index step size to sub-sample
+    differences : bool, optional, default = True
+        whether to take first temporal differences (True) or not (False)
+
+    Returns
+    -------
+    Xflat, as a 3D array
+        regular subsampling is applied if subsample > 1
+        differencing is applied (after subsampling) if differences=True
+    """
+    flatXshape = np.shape(Xflat)
+    Xshape = (flatXshape[0], numfeatures, flatXshape[1] / numfeatures)
+    X = np.reshape(Xflat, Xshape)[:, :, ::subsample]
+
+    if differences:
+        return np.diff(X)
+    else:
+        return X
+
+
+# scikit-learn estimator (original 2016 paper)
+# --------------------------------------------
+# for historical reasons - old scikit-learn version of the estimator
+class SeqKernelizer(BaseEstimator, TransformerMixin):
+    """Compute the sequential kernel matrix row features on collection of series.
+
+    Original implementation of signature kernel in [1]_ and [2]_.
+    The algorithm used therein is sklearn GridSearchCV wrapping
+    the sklearn pipeline of SeqKernelizer and SVC.
+
+    Identical behaviour can be obtained in sktime as
+    a grid search wrapped DistanceFeatures(SignatureKernel()) * SVC() pipeline
+
+    Included for historical purposes only, as reference to original paper code,
+    and for reproduction of the original experiments in the JMLR publication.
+    Users and developers should use/modify SequentialKernel instead.
+
+    This sklearn estimator requires passing of integer "numfeatures" as parameter,
+    and will interpret rows of X as time series with `numfeatures` features/vars,
+    and X.shape[1]/numfeatures time stamps, reshaped in (vars, time stamps) order.
+
+    In transform, will transform a series to the row of the kernel matrix
+    between that series and all the series seen in fit,, via seq_kernel_XY.
+
+    Parameters
+    ----------
+    level : int, optional, default = 2
+        an integer >= 1, representing the level of truncation of the sequential kernel
+    degree : int, optional, default = 1
+        an integer >= 1, representing the order of approximation of sequential kernel
+        can be set only if lowrank = False, otherwise ignored (always = 1)
+    theta : float, optional, default=1.0
+        a positive scaling factor for the levels, i-th level is scaled by theta^i
+    kernel : str, one of "linear", "Gauss", "Laplace", "poly"
+        code for inner kernel in the sequential kernel, with kernel parameters
+        "linear" - Euclidean kernel with scale parameter
+        "Gauss" - Gaussian kernel with scale parameter
+        "Laplace" - Laplace kernel with scale parameter
+        "poly" - polynomial kernel with degree deg and scale parameter
+    scale : float, optional, default = 1.0
+        a positive scaling factor for the inner kernel
+    degree : int, optional, default = 1, used only for polynomial kernel (kernel="poly")
+        degree of the polynomial kernel (if used)
+    numfeatures : int, optional, default = 2
+        number of features/variables in the time series
+    subsample : int, optional, default = 1
+        time index step size to sub-sample
+    differences : bool, optional, default = True
+        whether to take first temporal differences (True) or not (False)
+    normalize : bool, optional, default = False
+        whether the output kernel matrix is normalized
+        if True, sums and cumsums are divided by prod(K.shape)
+    lowrank : bool, optional, default = False
+        whether to use low rank approximation in computing the kernel
+    rankbound : int, optional, default = infinity
+        a hard threshold for the rank of the level matrices
+        used only if lowrank = True
+
+    References
+    ----------
+    .. [1] F. Kiraly, H. Oberhauser. 2016. "Kernels for sequentially ordered data.",
+        arXiv: 1601.08169.
+    .. [2] F. Kiraly, H. Oberhauser. 2019. "Kernels for sequentially ordered data.",
+        Journal of Machine Learning Research.
+    """
+
+    def __init__(
+        self,
+        level=2,
+        degree=1,
+        theta=1,
+        kernel="linear",
+        scale=1,
+        deg=2,
+        numfeatures=2,
+        subsample=100,
+        differences=True,
+        normalize=False,
+        lowrank=False,
+        rankbound=float("inf"),
+    ):
+        self.level = level
+        self.degree = degree
+        self.theta = theta
+        self.kernel = kernel
+        self.scale = scale
+        self.deg = deg
+        self.numfeatures = numfeatures
+        self.subsample = subsample
+        self.differences = differences
+        self.normalize = normalize
+        self.lowrank = lowrank
+        self.rankbound = rankbound
+
+        self._reshape_kwargs = {
+            "numfeatures": numfeatures,
+            "subsample": subsample,
+            "differences": differences,
+        }
+        self._kern_kwargs = {
+            "level": level,
+            "degree": degree,
+            "theta": theta,
+            "normalize": normalize,
+            "lowrank": lowrank,
+            "rankbound": rankbound,
+        }
+
+    def fit(self, X, y=None):
+        """Fit = reshape the series X."""
+        self._X = time_series_reshaper(X, **self._reshape_kwargs)
+        return self
+
+    def transform(self, X):
+        """Transform the data to kernel matrix rows."""
+        X = time_series_reshaper(X, **self._reshape_kwargs)
+
+        def kernselect(kername):
+            switcher = {
+                "linear": partial(k_euclid, scale=self.scale),
+                "Gauss": partial(k_gauss, scale=self.scale),
+                "Laplace": partial(k_laplace, scale=self.scale),
+                "poly": partial(k_polynom, scale=self.scale, deg=self.deg),
+            }
+            return switcher.get(kername, "nothing")
+
+        KSeq = seq_kernel_XY(X, self._X, kernselect(self.kernel), **self._kern_kwargs)
+
+        return KSeq
+
+
+# sktime interface - pairwise transformer
+# ---------------------------------------
+
+
+class SignatureKernel(BasePairwiseTransformerPanel):
+    """Compute the sequential kernel matrix row features on collection of series.
+
+    Implements the signature kernel of Kiraly et al, see [1]_ and [2]_,
+    including higher-order and low-rank approximation variants described therein.
+
+    Parameters
+    ----------
+    kernel : sktime pairwise (tabular) transformer, callable, or None
+        inner (tabular) kernel used in the signature sequence kernel
+        if callable: function (2D np.ndarray x 2D np.ndarray) -> 2D np.ndarray
+        pairwise kernel function, matrix sizes (n, d) x (m, d) -> (n x m)
+        optional, default = None = Euclidean (linear) kernel with scale parameter 1
+    level : int, optional, default = 2
+        an integer >= 1, representing the level of truncation of the sequential kernel
+    degree : int, optional, default = 1
+        an integer >= 1, representing the order of approximation of sequential kernel
+        can be set only if lowrank = False, otherwise ignored (always = 1)
+    theta : float, optional, default=1.0
+        a positive scaling factor for the levels, i-th level is scaled by theta^i
+    normalize : bool, optional, default = False
+        whether the output kernel matrix is normalized
+        if True, sums and cumsums are divided by prod(K.shape)
+    lowrank : bool, optional, default = False
+        whether to use low rank approximation in computing the kernel
+    rankbound : int, optional, default = infinity
+        a hard threshold for the rank of the level matrices
+        used only if lowrank = True
+
+    References
+    ----------
+    .. [1] F. Kiraly, H. Oberhauser. 2016. "Kernels for sequentially ordered data.",
+        arXiv: 1601.08169.
+    .. [2] F. Kiraly, H. Oberhauser. 2019. "Kernels for sequentially ordered data.",
+        Journal of Machine Learning Research.
+    """
+
+    _tags = {"X_inner_mtype": "numpy3D"}
+
+    def __init__(
+        self,
+        kernel=None,
+        level=2,
+        degree=1,
+        theta=1,
+        normalize=False,
+        lowrank=False,
+        rankbound=float("inf"),
+    ):
+        self.kernel = kernel
+        self.level = level
+        self.degree = degree
+        self.theta = theta
+        self.normalize = normalize
+        self.lowrank = lowrank
+        self.rankbound = rankbound
+
+        if kernel is None:
+            self._kernel = partial(k_euclid, scale=1)
+        else:
+            self._kernel = kernel
+
+        self._kern_kwargs = {
+            "L": level,
+            "D": degree,
+            "theta": theta,
+            "normalize": normalize,
+            "lowrank": lowrank,
+            "rankbound": rankbound,
+        }
+
+        super(SignatureKernel, self).__init__()
+
+    def _transform(self, X, X2=None):
+        """Compute distance/kernel matrix.
+
+        private _transform containing core logic, called from public transform
+
+        Behaviour: returns pairwise distance/kernel matrix
+            between samples in X and X2
+                if X2 is not passed, is equal to X
+                if X/X2 is a pd.DataFrame and contains non-numeric columns,
+                    these are removed before computation
+
+        Parameters
+        ----------
+        X: 3D np.array of shape [num_instances, num_vars, num_time_points]
+        X2: 3D np.array of shape [num_instances, num_vars, num_time_points], optional
+            default X2 = X
+
+        Returns
+        -------
+        distmat: np.array of shape [n, m]
+            (i,j)-th entry contains distance/kernel between X[i] and X2[j]
+        """
+        kwargs = self._kern_kwargs
+        kernel = self._kernel
+
+        if X2 is None:
+            return seq_kernel(X, kernelfun=kernel, **kwargs)
+        else:
+            return seq_kernel_XY(X, X2, kernelfun=kernel, **kwargs)
+
+    @classmethod
+    def get_test_params(cls, parameter_set="default"):
+        """Return testing parameter settings for the estimator.
+
+        Parameters
+        ----------
+        parameter_set : str, default="default"
+            Name of the set of test parameters to return, for use in tests. If no
+            special parameters are defined for a value, will return `"default"` set.
+            There are currently no reserved values for distance/kernel transformers.
+
+        Returns
+        -------
+        params : dict or list of dict, default = {}
+            Parameters to create testing instances of the class
+            Each dict are parameters to construct an "interesting" test instance, i.e.,
+            `MyClass(**params)` or `MyClass(**params[i])` creates a valid test instance.
+            `create_test_instance` uses the first (or only) dictionary in `params`
+        """
+        param1 = {}
+
+        # test higher level and normalization
+        param2 = {"level": 3, "normalize": True}
+
+        # test higher-order function
+        param3 = {"degree": 2}
+
+        # test low-rank approximation
+        param4 = {"lowrank": True}
+
+        paramlist = [param1, param2, param3, param4]
+
+        return paramlist
diff --git a/sktime/forecasting/arima.py b/sktime/forecasting/arima.py
index d457bdbc4..e8c292296 100644
--- a/sktime/forecasting/arima.py
+++ b/sktime/forecasting/arima.py
@@ -376,6 +376,8 @@ class AutoARIMA(_PmdArimaAdapter):
 
         super(AutoARIMA, self).__init__()
 
+        self._sp = sp if sp else 1
+
     def _instantiate_model(self):
         # import inside method to avoid hard dependency
         from pmdarima.arima import AutoARIMA as _AutoARIMA  # type: ignore
@@ -396,7 +398,7 @@ class AutoARIMA(_PmdArimaAdapter):
             max_D=self.max_D,
             max_Q=self.max_Q,
             max_order=self.max_order,
-            m=self.sp,
+            m=self._sp,
             seasonal=self.seasonal,
             stationary=self.stationary,
             information_criterion=self.information_criterion,
diff --git a/sktime/forecasting/base/_base.py b/sktime/forecasting/base/_base.py
index b001290d7..7f24a0673 100644
--- a/sktime/forecasting/base/_base.py
+++ b/sktime/forecasting/base/_base.py
@@ -1137,14 +1137,40 @@ class BaseForecaster(BaseEstimator):
     def get_fitted_params(self):
         """Get fitted parameters.
 
+        Overrides BaseEstimator default in case of vectorization.
+
         State required:
             Requires state to be "fitted".
 
         Returns
         -------
-        fitted_params : dict
+        fitted_params : dict of fitted parameters, keys are str names of parameters
+            parameters of components are indexed as [componentname]__[paramname]
         """
-        raise NotImplementedError("abstract method")
+        # if self is not vectorized, run the default get_fitted_params
+        if not getattr(self, "_is_vectorized", False):
+            return super(BaseForecaster, self).get_fitted_params()
+
+        # otherwise, we delegate to the instances' get_fitted_params
+        # instances' parameters are returned at dataframe-slice-like keys
+        fitted_params = {}
+
+        # forecasters contains a pd.DataFrame with the individual forecasters
+        forecasters = self.forecasters_
+
+        # return forecasters in the "forecasters" param
+        fitted_params["forecasters"] = forecasters
+
+        # populate fitted_params with forecasters and their parameters
+        for ix, col in zip(forecasters.index, forecasters.columns):
+            fcst = forecasters.loc[ix, col]
+            fcst_key = f"forecasters.loc[{ix},{col}]"
+            fitted_params[fcst_key] = fcst
+            fcst_params = fcst.get_fitted_params()
+            for key, val in fcst_params.items():
+                fitted_params[f"{fcst_key}__{key}"] = val
+
+        return fitted_params
 
     def _check_X_y(self, X=None, y=None):
         """Check and coerce X/y for fit/predict/update functions.
diff --git a/sktime/forecasting/base/_delegate.py b/sktime/forecasting/base/_delegate.py
index 18fa8f270..96431aab8 100644
--- a/sktime/forecasting/base/_delegate.py
+++ b/sktime/forecasting/base/_delegate.py
@@ -285,3 +285,18 @@ class _DelegatedForecaster(BaseForecaster):
         """
         estimator = self._get_delegate()
         return estimator.predict_proba(fh=fh, X=X, marginal=marginal)
+
+    def _get_fitted_params(self):
+        """Get fitted parameters.
+
+        private _get_fitted_params, called from get_fitted_params
+
+        State required:
+            Requires state to be "fitted".
+
+        Returns
+        -------
+        fitted_params : dict
+        """
+        estimator = self._get_delegate()
+        return estimator.get_fitted_params()
diff --git a/sktime/forecasting/base/adapters/_pmdarima.py b/sktime/forecasting/base/adapters/_pmdarima.py
index bd10d6bcd..def5b8363 100644
--- a/sktime/forecasting/base/adapters/_pmdarima.py
+++ b/sktime/forecasting/base/adapters/_pmdarima.py
@@ -171,10 +171,12 @@ class _PmdArimaAdapter(BaseForecaster):
                 pred_int.loc[fh_abs] = result[1][fh_idx, :]
                 pred_ints.append(pred_int)
             # unpack results
-            y_pred.loc[fh_abs] = result[0][fh_idx]
+            result = pd.Series(result[0]).iloc[fh_idx]
+            y_pred.loc[fh_abs] = result
             return y_pred, pred_ints
         else:
-            y_pred.loc[fh_abs] = result[fh_idx]
+            result = pd.Series(result).iloc[fh_idx]
+            y_pred.loc[fh_abs] = result
             return y_pred
 
     def _predict_fixed_cutoff(
@@ -220,7 +222,9 @@ class _PmdArimaAdapter(BaseForecaster):
                 pred_ints.append(pred_int)
             return result[0], pred_ints
         else:
-            return pd.Series(result[fh_idx], index=fh_abs)
+            result = pd.Series(result).iloc[fh_idx]
+            result.index = fh_abs
+            return result
 
     def _predict_interval(self, fh, X=None, coverage=0.90):
         """Compute/return prediction quantiles for a forecast.
diff --git a/sktime/forecasting/compose/_multiplexer.py b/sktime/forecasting/compose/_multiplexer.py
index f4dc47452..0a4bde71e 100644
--- a/sktime/forecasting/compose/_multiplexer.py
+++ b/sktime/forecasting/compose/_multiplexer.py
@@ -4,6 +4,7 @@
 """Implements forecaster for selecting among different model classes."""
 
 from sktime.base import _HeterogenousMetaEstimator
+from sktime.datatypes import ALL_TIME_SERIES_MTYPES
 from sktime.forecasting.base._base import BaseForecaster
 from sktime.forecasting.base._delegate import _DelegatedForecaster
 
@@ -83,7 +84,8 @@ class MultiplexForecaster(_DelegatedForecaster, _HeterogenousMetaEstimator):
         "requires-fh-in-fit": False,
         "handles-missing-data": False,
         "scitype:y": "both",
-        "y_inner_mtype": ["pd.DataFrame", "pd.Series"],
+        "y_inner_mtype": ALL_TIME_SERIES_MTYPES,
+        "X_inner_mtype": ALL_TIME_SERIES_MTYPES,
         "fit_is_empty": False,
     }
 
@@ -108,8 +110,12 @@ class MultiplexForecaster(_DelegatedForecaster, _HeterogenousMetaEstimator):
             clone_ests=False,
         )
         self._set_forecaster()
+
         self.clone_tags(self.forecaster_)
         self.set_tags(**{"fit_is_empty": False})
+        # this ensures that we convert in the inner estimator, not in the multiplexer
+        self.set_tags(**{"y_inner_mtype": ALL_TIME_SERIES_MTYPES})
+        self.set_tags(**{"X_inner_mtype": ALL_TIME_SERIES_MTYPES})
 
     @property
     def _forecasters(self):
diff --git a/sktime/forecasting/compose/_pipeline.py b/sktime/forecasting/compose/_pipeline.py
index 38f6a90f3..302286ad9 100644
--- a/sktime/forecasting/compose/_pipeline.py
+++ b/sktime/forecasting/compose/_pipeline.py
@@ -1129,8 +1129,8 @@ class ForecastX(BaseForecaster):
     """
 
     _tags = {
-        "X_inner_mtype": "pd.DataFrame",
-        "y_inner_mtype": "pd.DataFrame",
+        "X_inner_mtype": SUPPORTED_MTYPES,
+        "y_inner_mtype": SUPPORTED_MTYPES,
         "X-y-must-have-same-index": False,
         "fit_is_empty": False,
         "ignores-exogeneous-X": False,
diff --git a/sktime/forecasting/conformal.py b/sktime/forecasting/conformal.py
index 7b99b6e6e..ca8809a0e 100644
--- a/sktime/forecasting/conformal.py
+++ b/sktime/forecasting/conformal.py
@@ -8,6 +8,7 @@ Code based partially on NaiveVariance by ilyasmoutawwakil.
 __all__ = ["ConformalIntervals"]
 __author__ = ["fkiraly", "bethrice44"]
 
+from math import floor
 from warnings import warn
 
 import numpy as np
@@ -53,8 +54,12 @@ class ConformalIntervals(BaseForecaster):
             Caveat: this does not give frequentist but conformal predictive intervals
         "conformal": as in Stankeviciute et al, but with H=1,
             i.e., no Bonferroni correction under number of indices in the horizon
-    initial_window : int, optional, default=1
-        number of minimum initial indices to use for fitting when computing residuals
+    initial_window : float, int or None, optional (default=max(10, 0.1*len(y)))
+        Defines the size of the initial training window
+        If float, should be between 0.0 and 1.0 and represent the proportion
+        of the dataset to include for the initial window for the train split.
+        If int, represents the relative number of train samples in the initial window.
+        If None, the value is set to the larger of 0.1*len(y) and 10
     sample_frac : float, optional, default=None
         value in range (0,1) corresponding to fraction of y index to calculate
         residuals matrix values for (for speeding up calculation)
@@ -273,6 +278,39 @@ class ConformalIntervals(BaseForecaster):
 
         return pred_int
 
+    def _parse_initial_window(self, y, initial_window=None):
+
+        n_samples = len(y)
+
+        if initial_window is None:
+            initial_window = max(10, 0.1 * n_samples)
+
+        initial_window_type = np.asarray(initial_window).dtype.kind
+
+        if (
+            initial_window_type == "i"
+            and (initial_window >= n_samples or initial_window <= 0)
+            or initial_window_type == "f"
+            and (initial_window <= 0 or initial_window >= 1)
+        ):
+            raise ValueError(
+                "initial_window={0} should be either positive and smaller"
+                " than the number of samples {1} or a float in the "
+                "(0, 1) range".format(initial_window, n_samples)
+            )
+
+        if initial_window is not None and initial_window_type not in ("i", "f"):
+            raise ValueError(
+                "Invalid value for initial_window: {}".format(initial_window)
+            )
+
+        if initial_window_type == "f":
+            n_initial_window = int(floor(initial_window * n_samples))
+        elif initial_window_type == "i":
+            n_initial_window = int(initial_window)
+
+        return n_initial_window
+
     def _compute_sliding_residuals(self, y, X, forecaster, initial_window, sample_frac):
         """Compute sliding residuals used in uncertainty estimates.
 
@@ -284,10 +322,15 @@ class ConformalIntervals(BaseForecaster):
             sktime compatible exogeneous time series to use in forecasts
         forecaster : sktime compatible forecaster
             forecaster to use in computing the sliding residuals
-        initial_window : int
-            minimum length of initial window to use in fitting
+        initial_window : float, int or None, optional (default=max(10, 0.1*len(y)))
+            Defines the size of the initial training window
+            If float, should be between 0.0 and 1.0 and represent the proportion
+            of the dataset to include for the initial window for the train split.
+            If int, represents the relative number of train samples in the
+            initial window.
+            If None, the value is set to the larger of 0.1*len(y) and 10
         sample_frac : float
-            for speeding up computing of residuals matrix.
+            For speeding up computing of residuals matrix.
             sample value in range (0, 1) to obtain a fraction of y indices to
             compute residuals matrix for
         Returns
@@ -300,7 +343,9 @@ class ConformalIntervals(BaseForecaster):
         """
         y = convert_to(y, "pd.Series")
 
-        y_index = y.index[initial_window:]
+        n_initial_window = self._parse_initial_window(y, initial_window=initial_window)
+
+        y_index = y.iloc[n_initial_window:].index
 
         residuals_matrix = pd.DataFrame(columns=y_index, index=y_index, dtype="float")
 
diff --git a/sktime/forecasting/stream/_update.py b/sktime/forecasting/stream/_update.py
index 7ebec92ea..7a885b280 100644
--- a/sktime/forecasting/stream/_update.py
+++ b/sktime/forecasting/stream/_update.py
@@ -6,9 +6,17 @@ __author__ = ["fkiraly"]
 
 import pandas as pd
 
+from sktime.datatypes import ALL_TIME_SERIES_MTYPES
 from sktime.datatypes._utilities import get_window
 from sktime.forecasting.base._delegate import _DelegatedForecaster
 
+# prepare tags to clone - exceptions are TAGS_TO_KEEP
+TAGS_TO_KEEP = ["fit_is_empty", "X_inner_mtype", "y_inner_mtype"]
+# fit must be executed to fit the wrapped estimator and remember the cutoff
+# mtype tags are set so X/y is passed through, conversions happen in wrapped estimator
+TAGS_TO_CLONE = _DelegatedForecaster().get_tags().keys()
+TAGS_TO_CLONE = list(set(TAGS_TO_CLONE).difference(TAGS_TO_KEEP))
+
 
 class UpdateRefitsEvery(_DelegatedForecaster):
     """Refits periodically when update is called.
@@ -39,7 +47,12 @@ class UpdateRefitsEvery(_DelegatedForecaster):
 
     _delegate_name = "forecaster_"
 
-    _tags = {"fit_is_empty": False, "requires-fh-in-fit": False}
+    _tags = {
+        "fit_is_empty": False,
+        "requires-fh-in-fit": False,
+        "y_inner_mtype": ALL_TIME_SERIES_MTYPES,
+        "X_inner_mtype": ALL_TIME_SERIES_MTYPES,
+    }
 
     def __init__(
         self, forecaster, refit_interval=0, refit_window_size=None, refit_window_lag=0
@@ -53,10 +66,7 @@ class UpdateRefitsEvery(_DelegatedForecaster):
 
         super(UpdateRefitsEvery, self).__init__()
 
-        self.clone_tags(forecaster)
-
-        # fit must be executed to fit the wrapped estimator and remember the cutoff
-        self.set_tags(fit_is_empty=False)
+        self.clone_tags(forecaster, TAGS_TO_CLONE)
 
     def _fit(self, y, X=None, fh=None):
         """Fit forecaster to training data.
@@ -219,7 +229,12 @@ class UpdateEvery(_DelegatedForecaster):
 
     _delegate_name = "forecaster_"
 
-    _tags = {"fit_is_empty": False, "requires-fh-in-fit": False}
+    _tags = {
+        "fit_is_empty": False,
+        "requires-fh-in-fit": False,
+        "y_inner_mtype": ALL_TIME_SERIES_MTYPES,
+        "X_inner_mtype": ALL_TIME_SERIES_MTYPES,
+    }
 
     def __init__(self, forecaster, update_interval=None):
         self.forecaster = forecaster
@@ -229,10 +244,7 @@ class UpdateEvery(_DelegatedForecaster):
 
         super(UpdateEvery, self).__init__()
 
-        self.clone_tags(forecaster)
-
-        # fit must be executed to fit the wrapped estimator and remember the cutoff
-        self.set_tags(fit_is_empty=False)
+        self.clone_tags(forecaster, TAGS_TO_KEEP)
 
     def _fit(self, y, X=None, fh=None):
         """Fit forecaster to training data.
@@ -379,7 +391,12 @@ class DontUpdate(_DelegatedForecaster):
 
     _delegate_name = "forecaster_"
 
-    _tags = {"fit_is_empty": False, "requires-fh-in-fit": False}
+    _tags = {
+        "fit_is_empty": False,
+        "requires-fh-in-fit": False,
+        "y_inner_mtype": ALL_TIME_SERIES_MTYPES,
+        "X_inner_mtype": ALL_TIME_SERIES_MTYPES,
+    }
 
     def __init__(self, forecaster):
         self.forecaster = forecaster
@@ -387,7 +404,7 @@ class DontUpdate(_DelegatedForecaster):
 
         super(DontUpdate, self).__init__()
 
-        self.clone_tags(forecaster)
+        self.clone_tags(forecaster, TAGS_TO_CLONE)
 
     def _update(self, y, X=None, update_params=True):
         """Update time series to incremental training data.
diff --git a/sktime/param_est/seasonality.py b/sktime/param_est/seasonality.py
index 0e16b4d5d..4cb4f523e 100644
--- a/sktime/param_est/seasonality.py
+++ b/sktime/param_est/seasonality.py
@@ -50,6 +50,13 @@ class SeasonalityACF(BaseParamFitter):
             estimate the autocovariance. When using "conservative",
             n is set to the number of non-missing observations.
 
+    Attributes
+    ----------
+    sp_ : int, seasonality period at lowest p-level, if any sub-threshold, else 1
+        if `candidate_sp` is passed, will be in `candidate_sp` or 1
+    sp_significant_ : list of int, seasonality periods with sub-threshold p-levels
+        ordered increasingly by p-level. Empty list, not [1], if none are sub-threshold
+
     Examples
     --------
     >>> from sktime.datasets import load_airline
@@ -162,8 +169,8 @@ class SeasonalityACF(BaseParamFitter):
             self.sp_ = sp_significant[0]
             self.sp_significant_ = sp_significant
         else:
-            self.sp_ = None
-            self.sp_significant_ = None
+            self.sp_ = 1
+            self.sp_significant_ = []
 
         return self
 
@@ -240,6 +247,13 @@ class SeasonalityACFqstat(BaseParamFitter):
             estimate the autocovariance. When using "conservative",
             n is set to the number of non-missing observations.
 
+    Attributes
+    ----------
+    sp_ : int, seasonality period at lowest p-level, if any sub-threshold, else 1
+        if `candidate_sp` is passed, will be in `candidate_sp` or 1
+    sp_significant_ : list of int, seasonality periods with sub-threshold p-levels
+        ordered increasingly by p-level. Empty list, not [1], if none are sub-threshold
+
     Examples
     --------
     >>> from sktime.datasets import load_airline
@@ -358,8 +372,8 @@ class SeasonalityACFqstat(BaseParamFitter):
             self.sp_ = sp_significant[0]
             self.sp_significant_ = sp_significant
         else:
-            self.sp_ = None
-            self.sp_significant_ = None
+            self.sp_ = 1
+            self.sp_significant_ = []
 
         return self
 
diff --git a/sktime/regression/_delegate.py b/sktime/regression/_delegate.py
index 655fbeedc..1260ab08c 100644
--- a/sktime/regression/_delegate.py
+++ b/sktime/regression/_delegate.py
@@ -91,3 +91,18 @@ class _DelegatedRegressor(BaseRegressor):
         """
         estimator = self._get_delegate()
         return estimator.predict(X=X)
+
+    def _get_fitted_params(self):
+        """Get fitted parameters.
+
+        private _get_fitted_params, called from get_fitted_params
+
+        State required:
+            Requires state to be "fitted".
+
+        Returns
+        -------
+        fitted_params : dict
+        """
+        estimator = self._get_delegate()
+        return estimator.get_fitted_params()
diff --git a/sktime/transformations/_delegate.py b/sktime/transformations/_delegate.py
index c8055ba1c..1507b5ed8 100644
--- a/sktime/transformations/_delegate.py
+++ b/sktime/transformations/_delegate.py
@@ -166,3 +166,18 @@ class _DelegatedTransformer(BaseTransformer):
         estimator = self._get_delegate()
         estimator.update(y=y, X=X, update_params=update_params)
         return self
+
+    def _get_fitted_params(self):
+        """Get fitted parameters.
+
+        private _get_fitted_params, called from get_fitted_params
+
+        State required:
+            Requires state to be "fitted".
+
+        Returns
+        -------
+        fitted_params : dict
+        """
+        estimator = self._get_delegate()
+        return estimator.get_fitted_params()
diff --git a/sktime/transformations/base.py b/sktime/transformations/base.py
index 7e706ee5b..53b583fe4 100644
--- a/sktime/transformations/base.py
+++ b/sktime/transformations/base.py
@@ -407,6 +407,7 @@ class BaseTransformer(BaseEstimator):
         # checks and conversions complete, pass to inner fit
         #####################################################
         vectorization_needed = isinstance(X_inner, VectorizedDF)
+        self._is_vectorized = vectorization_needed
         # we call the ordinary _fit if no looping/vectorization needed
         if not vectorization_needed:
             self._fit(X=X_inner, y=y_inner)
@@ -678,6 +679,44 @@ class BaseTransformer(BaseEstimator):
 
         return self
 
+    def get_fitted_params(self):
+        """Get fitted parameters.
+
+        Overrides BaseEstimator default in case of vectorization.
+
+        State required:
+            Requires state to be "fitted".
+
+        Returns
+        -------
+        fitted_params : dict of fitted parameters, keys are str names of parameters
+            parameters of components are indexed as [componentname]__[paramname]
+        """
+        # if self is not vectorized, run the default get_fitted_params
+        if not getattr(self, "_is_vectorized", False):
+            return super(BaseTransformer, self).get_fitted_params()
+
+        # otherwise, we delegate to the instances' get_fitted_params
+        # instances' parameters are returned at dataframe-slice-like keys
+        fitted_params = {}
+
+        # transformers contains a pd.DataFrame with the individual transformers
+        transformers = self.transformers_
+
+        # return forecasters in the "forecasters" param
+        fitted_params["transformers"] = transformers
+
+        # populate fitted_params with ftransformers and their parameters
+        for ix, col in zip(transformers.index, transformers.columns):
+            fcst = transformers.loc[ix, col]
+            fcst_key = f"transformers.loc[{ix},{col}]"
+            fitted_params[fcst_key] = fcst
+            fcst_params = fcst.get_fitted_params()
+            for key, val in fcst_params.items():
+                fitted_params[f"{fcst_key}__{key}"] = val
+
+        return fitted_params
+
     def _check_X_y(self, X=None, y=None, return_metadata=False):
         """Check and coerce X/y for fit/transform functions.
 
diff --git a/sktime/transformations/compose.py b/sktime/transformations/compose.py
index cca2d9ea5..1a3107056 100644
--- a/sktime/transformations/compose.py
+++ b/sktime/transformations/compose.py
@@ -9,6 +9,7 @@ from sklearn import clone
 from sklearn.utils.metaestimators import if_delegate_has_method
 
 from sktime.base import _HeterogenousMetaEstimator
+from sktime.datatypes import ALL_TIME_SERIES_MTYPES
 from sktime.transformations._delegate import _DelegatedTransformer
 from sktime.transformations.base import BaseTransformer
 from sktime.utils.multiindex import flatten_multiindex
@@ -906,6 +907,7 @@ class MultiplexTransformer(_DelegatedTransformer, _HeterogenousMetaEstimator):
     _tags = {
         "fit_is_empty": False,
         "univariate-only": False,
+        "X_inner_mtype": ALL_TIME_SERIES_MTYPES,
     }
 
     _delegate_name = "transformer_"
@@ -928,6 +930,8 @@ class MultiplexTransformer(_DelegatedTransformer, _HeterogenousMetaEstimator):
         self._set_transformer()
         self.clone_tags(self.transformer_)
         self.set_tags(**{"fit_is_empty": False})
+        # this ensures that we convert in the inner estimator, not in the multiplexer
+        self.set_tags(**{"X_inner_mtype": ALL_TIME_SERIES_MTYPES})
 
     @property
     def _transformers(self):
diff --git a/sktime/transformations/panel/compose_distance.py b/sktime/transformations/panel/compose_distance.py
index 663e96cf2..e2d39fee5 100644
--- a/sktime/transformations/panel/compose_distance.py
+++ b/sktime/transformations/panel/compose_distance.py
@@ -64,6 +64,9 @@ class DistanceFeatures(BaseTransformer):
         "capability:unequal_length:removes": False,
         "handles-missing-data": True,
         "capability:missing_values:removes": False,
+        # we leave remember_data as False, since updating self._X in update
+        # would increase the number of columns in the transform return
+        "remember_data": False,
     }
 
     def __init__(self, distance=None, distance_mtype=None, flatten_hierarchy=False):
diff --git a/sktime/transformations/series/fourier.py b/sktime/transformations/series/fourier.py
new file mode 100644
index 000000000..8bfa17d62
--- /dev/null
+++ b/sktime/transformations/series/fourier.py
@@ -0,0 +1,212 @@
+# -*- coding: utf-8 -*-
+# copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
+"""Fourier features for time series with long/complex seasonality."""
+
+__author__ = ["ltsaprounis"]
+
+import warnings
+from copy import deepcopy
+from typing import List, Union
+
+import numpy as np
+import pandas as pd
+
+from sktime.transformations.base import BaseTransformer
+
+
+class FourierFeatures(BaseTransformer):
+    r"""Fourier Features for time series seasonality.
+
+    Fourier Series terms can be used as explanatory variables for the cases of multiple
+    seasonal periods and or complex / long seasonal periods [1]_, [2]_. For every
+    seasonal period, :math:`sp` and fourier term :math:`k` pair there are 2 fourier
+    terms sin_sp_k and cos_sp_k:
+        - sin_sp_k = :math:`sin(\frac{2 \pi k t}{sp})`
+        - cos_sp_k = :math:`cos(\frac{2 \pi k t}{sp})`
+
+    Where :math:`t` is the number of time steps elapsed from the beginning of the time
+    series.
+
+    The output of the trasform is a pandas DataFrame that includes the fourier terms as
+    additional columns with the naming convention stated above (sin_sp_k and cos_sp_k).
+    For instance for sp_list = [12, 3] and fourier_terms_list = [2, 1] the transformed
+    series will have the additional columns:
+    "cos_12_1", "sin_12_1", "cos_12_2", "sin_12_2", "cos_3_1", "sin_3_1"
+
+    The implementation is based on the fourier function from the R forecast package [3]_
+
+    Parameters
+    ----------
+    sp_list : List[Union[int, float]]
+        list of seasonal periods
+    fourier_terms_list : List[int]
+        list of number of fourier terms (K) for each seasonal period.
+        Each K matches to the sp (seasonal period) of the sp_list.
+        For example, if sp_list = [7, 365] and fourier_terms_list = [3, 9], the seasonal
+        frequency of 7 will have 3 fourier terms and the seasonal frequency of 365
+        will have 9 fourier terms.
+
+    References
+    ----------
+    .. [1] Hyndsight - Forecasting with long seasonal periods:
+        https://robjhyndman.com/hyndsight/longseasonality/
+    .. [2] Hyndman, R.J., & Athanasopoulos, G. (2021) Forecasting: principles and
+        practice, 3rd edition, OTexts: Melbourne, Australia. OTexts.com/fpp3.
+        Accessed on August 14th 2022.
+    .. [3] https://pkg.robjhyndman.com/forecast/reference/fourier.html
+
+    Examples
+    --------
+    >>> from sktime.transformations.series.fourier import FourierFeatures
+    >>> from sktime.datasets import load_airline
+    >>> y = load_airline()
+    >>> transformer = FourierFeatures(sp_list=[12], fourier_terms_list=[4])
+    >>> y_hat = transformer.fit_transform(y)
+    """
+
+    _tags = {
+        "scitype:transform-input": "Series",
+        # what is the scitype of X: Series, or Panel
+        "scitype:transform-output": "Series",
+        # what scitype is returned: Primitives, Series, Panel
+        "scitype:transform-labels": "None",
+        # what is the scitype of y: None (not needed), Primitives, Series, Panel
+        "scitype:instancewise": True,  # is this an instance-wise transform?
+        "capability:inverse_transform": False,  # can the transformer inverse transform?
+        "univariate-only": False,  # can the transformer handle multivariate X?
+        "X_inner_mtype": "pd.DataFrame",  # which mtypes do _fit/_predict support for X?
+        # this can be a Panel mtype even if transform-input is Series, vectorized
+        "y_inner_mtype": "None",  # which mtypes do _fit/_predict support for y?
+        "requires_y": False,  # does y need to be passed in fit?
+        "enforce_index_type": [pd.PeriodIndex],  # index type that needs to be enforced
+        # in X/y
+        "fit_is_empty": False,  # is fit empty and can be skipped? Yes = True
+        "X-y-must-have-same-index": False,  # can estimator handle different X/y index?
+        "transform-returns-same-time-index": True,
+        # does transform return have the same time index as input X
+        "skip-inverse-transform": True,  # is inverse-transform skipped when called?
+        "capability:unequal_length": False,
+        # can the transformer handle unequal length time series (if passed Panel)?
+        "capability:unequal_length:removes": False,
+        # is transform result always guaranteed to be equal length (and series)?
+        #   not relevant for transformers that return Primitives in transform-output
+        "handles-missing-data": False,  # can estimator handle missing data?
+        # todo: rename to capability:missing_values
+        "capability:missing_values:removes": False,
+        # is transform result always guaranteed to contain no missing values?
+        "python_version": None,  # PEP 440 python version specifier to limit versions
+    }
+
+    def __init__(self, sp_list: List[Union[int, float]], fourier_terms_list: List[int]):
+        self.sp_list = sp_list
+        self.fourier_terms_list = fourier_terms_list
+
+        if len(self.sp_list) != len(self.fourier_terms_list):
+            raise ValueError(
+                "In FourierFeatures the length of the sp_list needs to be equal "
+                "to the length of fourier_terms_list."
+            )
+
+        if np.any(np.array(self.sp_list) / np.array(self.fourier_terms_list) < 1):
+            raise ValueError(
+                "In FourierFeatures the number of each element of fourier_terms_list"
+                "needs to be lower from the corresponding element of the sp_list"
+            )
+
+        super(FourierFeatures, self).__init__()
+
+    def _fit(self, X, y=None):
+        """Fit transformer to X and y.
+
+        private _fit containing the core logic, called from fit
+
+        Parameters
+        ----------
+        X : Series or Panel of mtype X_inner_mtype
+            if X_inner_mtype is list, _fit must support all types in it
+            Data to fit transform to
+        y : Series or Panel of mtype y_inner_mtype, default=None
+            Additional data, e.g., labels for transformation
+
+        Returns
+        -------
+        self: reference to self
+        """
+        # Create the sp, k pairs
+        # Don't add pairs where the coefficient k/sp already exists
+        self.sp_k_pairs_list_ = []
+        coefficient_list = []
+        for i, sp in enumerate(self.sp_list):
+            for k in range(1, self.fourier_terms_list[i] + 1):
+                coef = k / sp
+                if coef not in coefficient_list:
+                    coefficient_list.append(coef)
+                    self.sp_k_pairs_list_.append((sp, k))
+                else:
+                    warnings.warn(
+                        f"The terms sin_{sp}_{k} and cos_{sp}_{k} from FourierFeatures "
+                        "will be skipped because the resulting coefficient already "
+                        "exists from other seasonal period, fourier term pairs."
+                    )
+
+        # store the integer form of the minimum date in the prediod index
+        # this is used to make sure that time t is calculated with reference to
+        # the data passed on fit
+        self.min_t_ = np.min(X.index.astype(int))
+
+        return self
+
+    def _transform(self, X, y=None):
+        """Transform X and return a transformed version.
+
+        private _transform containing core logic, called from transform
+
+        Parameters
+        ----------
+        X : Series or Panel of mtype X_inner_mtype
+            if X_inner_mtype is list, _transform must support all types in it
+            Data to be transformed
+        y : Series or Panel of mtype y_inner_mtype, default=None
+            Additional data, e.g., labels for transformation
+
+        Returns
+        -------
+        transformed version of X
+        """
+        X_transformed = deepcopy(X)
+        # get the integer form of the PeriodIndex
+        int_index = X_transformed.index.astype(int) - self.min_t_
+
+        for sp_k in self.sp_k_pairs_list_:
+            sp = sp_k[0]
+            k = sp_k[1]
+
+            X_transformed[f"sin_{sp}_{k}"] = np.sin(int_index * 2 * k * np.pi / sp)
+            X_transformed[f"cos_{sp}_{k}"] = np.cos(int_index * 2 * k * np.pi / sp)
+
+        return X_transformed
+
+    @classmethod
+    def get_test_params(cls, parameter_set="default"):
+        """Return testing parameter settings for the estimator.
+
+        Parameters
+        ----------
+        parameter_set : str, default="default"
+            Name of the set of test parameters to return, for use in tests. If no
+            special parameters are defined for a value, will return `"default"` set.
+            There are currently no reserved values for transformers.
+
+        Returns
+        -------
+        params : dict or list of dict, default = {}
+            Parameters to create testing instances of the class
+            Each dict are parameters to construct an "interesting" test instance, i.e.,
+            `MyClass(**params)` or `MyClass(**params[i])` creates a valid test instance.
+            `create_test_instance` uses the first (or only) dictionary in `params`
+        """
+        params = [
+            {"sp_list": [12], "fourier_terms_list": [4]},
+            {"sp_list": [12, 6.2], "fourier_terms_list": [3, 4]},
+        ]
+        return params


Test patch:
diff --git a/sktime/annotation/tests/test_stray.py b/sktime/annotation/tests/test_stray.py
new file mode 100644
index 000000000..bcb6127ce
--- /dev/null
+++ b/sktime/annotation/tests/test_stray.py
@@ -0,0 +1,289 @@
+# -*- coding: utf-8 -*-
+"""Tests for STRAY outlier estimator."""
+
+__author__ = ["KatieBuc"]
+
+import numpy as np
+from sklearn.preprocessing import MinMaxScaler, RobustScaler
+
+from sktime.annotation.stray import STRAY
+
+
+def test_default_1D():
+    """Test with default parameters and 1D input array."""
+    X = np.array(
+        [
+            -7.207066,
+            -5.722571,
+            -4.915559,
+            -8.345698,
+            -5.570875,
+            -5.493944,
+            -6.574740,
+            -6.546632,
+            -6.564452,
+            -6.890038,
+            -6.477193,
+            -6.998386,
+            -6.776254,
+            -5.935541,
+            -5.040506,
+            0.000000,
+            5.889715,
+            5.488990,
+            5.088805,
+            5.162828,
+            8.415835,
+            6.134088,
+            5.509314,
+            5.559452,
+            6.459589,
+            5.306280,
+            4.551795,
+            6.574756,
+            4.976344,
+            5.984862,
+            5.064051,
+            7.102298,
+        ]
+    )
+
+    y_expected = np.array(
+        [
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            True,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+        ]
+    )
+
+    scaler = MinMaxScaler()
+    X = scaler.fit_transform(X.reshape(-1, 1))
+    model = STRAY()
+    y_actual = model.fit_transform(X)
+    assert np.allclose(y_actual, y_expected)
+
+
+def test_default_2D():
+    """Test with default parameters and 2D input array."""
+    X = np.array(
+        [
+            [-1.207, -0.776, -0.694],
+            [0.277, 0.064, -1.448],
+            [1.084, 0.959, 0.575],
+            [-2.346, -0.11, -1.024],
+            [0.429, -0.511, -0.015],
+            [0.506, -0.911, -0.936],
+            [-0.575, -0.837, 1.102],
+            [-0.547, 2.416, -0.476],
+            [-0.564, 0.134, -0.709],
+            [-0.89, -0.491, -0.501],
+            [-0.477, -0.441, -1.629],
+            [-0.998, 0.46, -1.168],
+            [10.0, 12.0, 10.0],
+            [3.0, 7.0, 10.0],
+        ]
+    )
+
+    y_expected = np.array(
+        [
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            True,
+            True,
+        ]
+    )
+
+    scaler = MinMaxScaler()
+    X = scaler.fit_transform(X)
+    model = STRAY()
+    y_actual = model.fit_transform(X)
+    assert np.allclose(y_actual, y_expected)
+
+
+def test_1D_score_with_na():
+    """Test score with 1D input array with missing values."""
+    X = np.array(
+        [
+            np.nan,
+            -5.72257076,
+            -4.91555882,
+            -8.3456977,
+            -5.57087531,
+            0.0,
+            6.50605589,
+            5.42526004,
+            5.45336814,
+            5.435548,
+            5.10996217,
+        ]
+    )
+
+    y_scores_expected = np.array(
+        [
+            np.nan,
+            0.17662069,
+            0.23095851,
+            0.17662069,
+            0.18683466,
+            0.33097498,
+            0.07087969,
+            0.02122967,
+            0.02312225,
+            0.02192238,
+            0.02122967,
+        ]
+    )
+
+    scaler = MinMaxScaler()
+    X = scaler.fit_transform(X.reshape(-1, 1))
+    model = STRAY(k=3)
+    fitted_model = model.fit(X)
+    y_scores_actual = fitted_model.score_
+    assert np.allclose(y_scores_actual, y_scores_expected, equal_nan=True)
+
+
+def test_1D_bool_with_na():
+    """Test anomaly detection with 1D input array with missing values."""
+    X = np.array(
+        [
+            np.nan,
+            -5.72257076,
+            -4.91555882,
+            -8.3456977,
+            -5.57087531,
+            0.0,
+            6.50605589,
+            5.42526004,
+            5.45336814,
+            5.435548,
+            5.10996217,
+        ]
+    )
+
+    y_expected = np.array([0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0])
+
+    scaler = MinMaxScaler()
+    X = scaler.fit_transform(X.reshape(-1, 1))
+    model = STRAY(k=3)
+    fitted_model = model.fit(X)
+    y_actual = fitted_model.y_
+    assert np.allclose(y_actual, y_expected)
+
+
+def test_2D_score_with_na():
+    """Test score with 2D input array with missing values."""
+    X = np.array(
+        [
+            [-1.20706575, -0.57473996],
+            [0.27742924, -0.54663186],
+            [1.08444118, np.nan],
+            [-2.3456977, -0.89003783],
+            [0.42912469, -0.4771927],
+            [0.50605589, -0.99838644],
+        ]
+    )
+
+    y_scores_expected = np.array(
+        [0.43612712, 0.43612712, np.nan, 0.69004259, 0.51240831, 0.85697804]
+    )
+
+    scaler = MinMaxScaler()
+    X = scaler.fit_transform(X)
+    model = STRAY(k=2, size_threshold=4)
+    fitted_model = model.fit(X)
+    y_scores_actual = fitted_model.score_
+    assert np.allclose(y_scores_actual, y_scores_expected, equal_nan=True)
+
+
+def test_2D_bool_with_na():
+    """Test anomaly detection with 2D input array with missing values."""
+    X = np.array(
+        [
+            [-1.20706575, -0.57473996],
+            [0.27742924, -0.54663186],
+            [1.08444118, np.nan],
+            [-2.3456977, -0.89003783],
+            [0.42912469, -0.4771927],
+            [0.50605589, -0.99838644],
+        ]
+    )
+
+    y_expected = np.array([0, 0, 0, 1, 1, 1])
+
+    scaler = MinMaxScaler()
+    X = scaler.fit_transform(X)
+    model = STRAY(k=2, size_threshold=4)
+    fitted_model = model.fit(X)
+    y_actual = fitted_model.y_
+    assert np.allclose(y_actual, y_expected)
+
+
+def test_2D_score_with_standardize():
+    """Test score with 2D input array and median/IQR normalization."""
+    X = np.array(
+        [
+            [-1.20706575, -0.57473996],
+            [0.27742924, -0.54663186],
+            [1.08444118, -0.5644520],
+            [-2.3456977, -0.89003783],
+            [0.42912469, -0.4771927],
+            [0.50605589, -0.99838644],
+        ]
+    )
+
+    y_scores_expected = np.array(
+        [
+            1.1274565,
+            0.6139288,
+            0.5982989,
+            1.4866554,
+            0.5982989,
+            1.7245212,
+        ]
+    )
+
+    scaler = RobustScaler()
+    X = scaler.fit_transform(X)
+    model = STRAY(k=2, size_threshold=4)
+    fitted_model = model.fit(X)
+    y_scores_actual = fitted_model.score_
+    assert np.allclose(y_scores_actual, y_scores_expected)
diff --git a/sktime/transformations/series/tests/test_fourier.py b/sktime/transformations/series/tests/test_fourier.py
new file mode 100644
index 000000000..cbca0520b
--- /dev/null
+++ b/sktime/transformations/series/tests/test_fourier.py
@@ -0,0 +1,74 @@
+# -*- coding: utf-8 -*-
+"""Tests for the FourierFeatures transformer."""
+
+import numpy as np
+import pandas as pd
+import pytest
+from pandas.util.testing import assert_frame_equal
+
+from sktime.datasets import load_airline
+from sktime.transformations.series.fourier import FourierFeatures
+
+Y = load_airline()
+
+
+def test_fourier_list_length_missmatch():
+    """Tests exception raised when sp_list & fourier_terms_list lengths don't match."""
+    with pytest.raises(ValueError) as ex:
+        FourierFeatures(sp_list=[365, 52], fourier_terms_list=[1])
+        assert ex.value == (
+            "In FourierFeatures the length of the sp_list needs to be equal "
+            "to the length of fourier_terms_list."
+        )
+
+
+def test_fourier_k_larger_than_sp():
+    """Tests exception raised when fourier_terms_list elements larger than sp_list."""
+    with pytest.raises(ValueError) as ex:
+        FourierFeatures(sp_list=[2], fourier_terms_list=[3])
+        assert ex.value == (
+            "In FourierFeatures the number of each element of fourier_terms_list"
+            "needs to be lower from the corresponding element of the sp_list"
+        )
+
+
+def test_fourier_redundant_terms_dropped():
+    """Tests redundant sp and k pairs are dropped if their equivalents exist."""
+    transformer = FourierFeatures(sp_list=[12, 9, 3], fourier_terms_list=[4, 3, 1])
+    transformer.fit(Y)
+    assert transformer.sp_k_pairs_list_ == [
+        (12, 1),
+        (12, 2),
+        (12, 3),
+        (12, 4),
+        (9, 1),
+        (9, 2),
+    ]
+
+
+def test_fit_transform_outputs():
+    """Tests that we get the expected outputs."""
+    y = Y.iloc[:3]
+    y_transformed = FourierFeatures(sp_list=[12], fourier_terms_list=[2]).fit_transform(
+        y
+    )
+    expected = (
+        y.to_frame()
+        .assign(sin_12_1=[np.sin(2 * np.pi * i / 12) for i in range(3)])
+        .assign(cos_12_1=[np.cos(2 * np.pi * i / 12) for i in range(3)])
+        .assign(sin_12_2=[np.sin(4 * np.pi * i / 12) for i in range(3)])
+        .assign(cos_12_2=[np.cos(4 * np.pi * i / 12) for i in range(3)])
+    )
+    assert_frame_equal(y_transformed, expected)
+
+
+def test_fit_transform_behaviour():
+    """Tests that the transform method evaluates time steps passed based on X in fit."""
+    transformer = FourierFeatures(sp_list=[12], fourier_terms_list=[2])
+    # fit_transform on one part of the dataset
+    y_tr_1 = transformer.fit_transform(Y.iloc[:100])
+    # transform the rest
+    y_tr_2 = transformer.transform(Y.iloc[100:])
+    # fit_transform the entire dataset
+    y_tr_complete = transformer.fit_transform(Y)
+    assert_frame_equal(pd.concat([y_tr_1, y_tr_2]), y_tr_complete)
diff --git a/sktime/utils/_testing/scenarios_transformers.py b/sktime/utils/_testing/scenarios_transformers.py
index a6835f8a9..1adef6cb3 100644
--- a/sktime/utils/_testing/scenarios_transformers.py
+++ b/sktime/utils/_testing/scenarios_transformers.py
@@ -195,8 +195,8 @@ class TransformerFitTransformSeriesUnivariate(TransformerTestScenario):
     }
 
     args = {
-        "fit": {"X": _make_series(n_timepoints=10, random_state=RAND_SEED)},
-        "transform": {"X": _make_series(n_timepoints=10, random_state=RAND_SEED)},
+        "fit": {"X": _make_series(n_timepoints=11, random_state=RAND_SEED)},
+        "transform": {"X": _make_series(n_timepoints=11, random_state=RAND_SEED2)},
         # "inverse_transform": {"X": _make_series(n_timepoints=10)},
     }
     default_method_sequence = ["fit", "transform"]