Repository: alan-turing-institute/sktime
Pull request #2972 (2022-07-10 21:15:17)

[ENH] add check for forecast to have correct columns
This adds a test that forecast `predict` output has the same column names as the training data, if both have column names.

Surprisingly, this was not yet tested.

Relies on PRs with fixes of bugs that would make this fail:
* `Prophet.predict` https://github.com/alan-turing-institute/sktime/pull/2973

aiwalter: Should this also je testted for transformers in `transform` and `inverse_transform`?
fkiraly: > Should this also je testted for transformers in `transform` and `inverse_transform`?

Not all transformers are 1:1 in terms of columns, or column names.

The counterpart for transformers are these:
* uniqueness (weaker but catches the most common problem) https://github.com/alan-turing-institute/sktime/pull/2969
* automated naming https://github.com/alan-turing-institute/sktime/pull/2958

aiwalter: But we could test that after `transform` _and_ `inverse_transform` we  have the same columns?
fkiraly: > But we could test that after `transform` _and_ `inverse_transform` we  have the same columns?

Yes, good idea.
Naively, I thought this was already tested, but indeed, `test_transform_inverse_transform_equivalent` does not seem to check for that.
I fear there might be pandora's box opening if we add the check to the test ... perhaps we may like to do it in a separate PR...
fkiraly: should we merge this and #2973, @aiwalter?
We can approximate step by step the desired end state of stronger output checks across estimator types.

Patch:
diff --git a/docs/source/api_reference/forecasting.rst b/docs/source/api_reference/forecasting.rst
index b7a7a153e..f04d76e56 100644
--- a/docs/source/api_reference/forecasting.rst
+++ b/docs/source/api_reference/forecasting.rst
@@ -201,6 +201,14 @@ All "ARIMA" models below include SARIMAX capability.
 
     VAR
 
+.. currentmodule:: sktime.forecasting.varmax
+
+.. autosummary::
+    :toctree: auto_generated/
+    :template: class.rst
+
+    VARMAX
+
 Structural time series models
 -----------------------------
 
diff --git a/examples/01_forecasting.ipynb b/examples/01_forecasting.ipynb
index 51e6bc5b0..062192828 100644
--- a/examples/01_forecasting.ipynb
+++ b/examples/01_forecasting.ipynb
@@ -1128,7 +1128,7 @@
     "#  e.g., the non-symmetric verion\n",
     "# it also allows for inspection of metric properties\n",
     "#  e.g., are higher values better (answer: no)?\n",
-    "mape.greater_is_better"
+    "mape.get_tag(\"lower_is_better\")"
    ]
   },
   {
diff --git a/setup.cfg b/setup.cfg
index 2ea7039d1..117fd3fa0 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -20,7 +20,7 @@ addopts =
     --showlocals
     --pyargs
     --verbose
-    -n 2
+    -n auto
 filterwarnings =
     ignore::UserWarning
     ignore:numpy.dtype size changed
diff --git a/sktime/base/_base.py b/sktime/base/_base.py
index 9dea5fa86..a6abfe77f 100644
--- a/sktime/base/_base.py
+++ b/sktime/base/_base.py
@@ -630,14 +630,6 @@ class TagAliaserMixin:
     ensure toremove this class as a parent of BaseObject or BaseEstimator.
     """
 
-    # dictionary of aliases
-    # key = old tag; value = new tag, aliased by old tag
-    alias_dict = {"capability:early_prediction": ""}
-
-    # dictionary of removal version
-    # key = old tag; value = version in which tag will be removed, as string
-    deprecate_dict = {"capability:early_prediction": "0.13.0"}
-
     def __init__(self):
         super(TagAliaserMixin, self).__init__()
 
@@ -795,7 +787,7 @@ class TagAliaserMixin:
                 warnings.warn(msg, category=DeprecationWarning)
 
 
-class BaseEstimator(TagAliaserMixin, BaseObject):
+class BaseEstimator(BaseObject):
     """Base class for defining estimators in sktime.
 
     Extends sktime's BaseObject to include basic functionality for fittable estimators.
diff --git a/sktime/classification/feature_based/_fresh_prince.py b/sktime/classification/feature_based/_fresh_prince.py
index b8cd614a9..6f39d80d1 100644
--- a/sktime/classification/feature_based/_fresh_prince.py
+++ b/sktime/classification/feature_based/_fresh_prince.py
@@ -66,6 +66,7 @@ class FreshPRINCE(BaseClassifier):
         "capability:train_estimate": True,
         "classifier_type": "feature",
         "python_version": "<3.10",
+        "python_dependencies": "tsfresh",
     }
 
     def __init__(
diff --git a/sktime/classification/feature_based/_tsfresh_classifier.py b/sktime/classification/feature_based/_tsfresh_classifier.py
index 94e267674..6e98bc349 100644
--- a/sktime/classification/feature_based/_tsfresh_classifier.py
+++ b/sktime/classification/feature_based/_tsfresh_classifier.py
@@ -69,6 +69,7 @@ class TSFreshClassifier(BaseClassifier):
         "capability:multithreading": True,
         "classifier_type": "feature",
         "python_version": "<3.10",
+        "python_dependencies": "tsfresh",
     }
 
     def __init__(
diff --git a/sktime/datatypes/_panel/_convert.py b/sktime/datatypes/_panel/_convert.py
index 976f11491..7b8df36d2 100644
--- a/sktime/datatypes/_panel/_convert.py
+++ b/sktime/datatypes/_panel/_convert.py
@@ -683,8 +683,8 @@ def from_multi_index_to_nested(
     multi_ind_dataframe : pd.DataFrame
         Input multi-indexed pandas DataFrame
 
-    instance_index_name : str
-        The name of multi-index level corresponding to the DataFrame's instances
+    instance_index_name : int or str, default=0 (first level = 0-th index)
+        Index or name of multi-index level corresponding to the DataFrame's instances
 
     cells_as_numpy : bool, default = False
         If True, then nested cells contain NumPy array
@@ -764,54 +764,21 @@ def from_nested_to_multi_index(X, instance_index=None, time_index=None):
         The multi-indexed pandas DataFrame
 
     """
-    if time_index is None:
-        time_index_name = "timepoints"
-    else:
-        time_index_name = time_index
+    # this contains the right values, but does not have the right index
+    #   need convert_dtypes or dtypes will always be object
+    # explode by column to ensure we deal with unequal length series properly
+    X_mi = pd.DataFrame()
 
-    # n_columns = X.shape[1]
-    nested_col_mask = [*are_columns_nested(X)]
+    for c in X.columns:
+        X_col = X[[c]].explode(c)
+        X_col = X_col.infer_objects()
 
-    if instance_index is None:
-        instance_idxs = X.index.get_level_values(-1).unique()
-        # n_instances = instance_idxs.shape[0]
-        instance_index_name = "instance"
-
-    else:
-        if instance_index in X.index.names:
-            instance_idxs = X.index.get_level_values(instance_index).unique()
-        else:
-            instance_idxs = X.index.get_level_values(-1).unique()
-        # n_instances = instance_idxs.shape[0]
-        instance_index_name = instance_index
-
-    instances = []
-    for instance_idx in instance_idxs:
-        iidx = instance_idx
-        series = [i[1] for i in X.loc[iidx, :].iteritems()]
-        colnames = [i[0] for i in X.loc[iidx, :].iteritems()]
-        for x in series:
-            if hasattr(x, "name"):
-                x.name = None
-
-        instance = [pd.DataFrame(s, columns=[c]) for s, c in zip(series, colnames)]
-        instance = pd.concat(instance, axis=1)
-        # For primitive (non-nested column) assume the same
-        # primitive value applies to every timepoint of the instance
-        for col_idx, is_nested in enumerate(nested_col_mask):
-            if not is_nested:
-                instance.iloc[:, col_idx] = instance.iloc[:, col_idx].ffill()
-
-        # Correctly assign multi-index
-        multi_index = pd.MultiIndex.from_product(
-            [[instance_idx], instance.index],
-            names=[instance_index_name, time_index_name],
-        )
-        instance.index = multi_index
-        instances.append(instance)
+        # create the right MultiIndex and assign to X_mi
+        idx_df = X[[c]].applymap(lambda x: x.index).explode(c)
+        idx_df = idx_df.set_index(c, append=True)
+        X_col.index = idx_df.index.set_names([instance_index, time_index])
 
-    X_mi = pd.concat(instances)
-    X_mi.columns = X.columns
+        X_mi[[c]] = X_col
 
     return X_mi
 
diff --git a/sktime/datatypes/_series/_convert.py b/sktime/datatypes/_series/_convert.py
index 61881e937..fea20e113 100644
--- a/sktime/datatypes/_series/_convert.py
+++ b/sktime/datatypes/_series/_convert.py
@@ -100,7 +100,7 @@ def convert_MvS_to_np_as_Series(obj: pd.DataFrame, store=None) -> np.ndarray:
         store["columns"] = obj.columns
         store["index"] = obj.index
 
-    return obj.to_numpy()
+    return obj.to_numpy(dtype="float")
 
 
 convert_dict[("pd.DataFrame", "np.ndarray", "Series")] = convert_MvS_to_np_as_Series
@@ -114,7 +114,7 @@ def convert_UvS_to_np_as_Series(obj: pd.Series, store=None) -> np.ndarray:
     if isinstance(store, dict):
         store["index"] = obj.index
 
-    return pd.DataFrame(obj).to_numpy()
+    return pd.DataFrame(obj).to_numpy(dtype="float")
 
 
 convert_dict[("pd.Series", "np.ndarray", "Series")] = convert_UvS_to_np_as_Series
diff --git a/sktime/datatypes/_vectorize.py b/sktime/datatypes/_vectorize.py
index 5a169ec47..9ad6d6a23 100644
--- a/sktime/datatypes/_vectorize.py
+++ b/sktime/datatypes/_vectorize.py
@@ -10,6 +10,7 @@ import pandas as pd
 
 from sktime.datatypes._check import check_is_scitype, mtype
 from sktime.datatypes._convert import convert_to
+from sktime.utils.multiindex import flatten_multiindex
 
 
 class VectorizedDF:
@@ -105,37 +106,34 @@ class VectorizedDF:
         self.X_multiindex = self._init_conversion(X)
         self.iter_indices = self._init_iter_indices()
 
-    def _init_conversion(self, X):
-        """Convert X to a pandas multiindex format."""
-        is_scitype = self.is_scitype
+    def _coerce_to_df(self, obj, scitype=None, store=None, store_behaviour=None):
+        """Coerce obj to a pandas multiindex format."""
+        pandas_dict = {
+            "Series": "pd.DataFrame",
+            "Panel": "pd-multiindex",
+            "Hierarchical": "pd_multiindex_hier",
+            None: ["pd.DataFrame", "pd-multiindex", "pd_multiindex_hier"],
+        }
 
-        if is_scitype == "Series":
-            return convert_to(
-                X,
-                to_type="pd.DataFrame",
-                as_scitype="Series",
-                store=self.converter_store,
-            )
-        elif is_scitype == "Panel":
-            return convert_to(
-                X,
-                to_type="pd-multiindex",
-                as_scitype="Panel",
-                store=self.converter_store,
-            )
-        elif is_scitype == "Hierarchical":
-            return convert_to(
-                X,
-                to_type="pd_multiindex_hier",
-                as_scitype="Hierarchical",
-                store=self.converter_store,
-            )
-        else:
+        if scitype not in pandas_dict.keys():
             raise RuntimeError(
-                f"unexpected value found for attribute self.is_scitype: {is_scitype}"
-                'must be "Panel" or "Hierarchical"'
+                f"unexpected value found for attribute scitype: {scitype}"
+                f"must be one of {pandas_dict.keys()}"
             )
 
+        return convert_to(
+            obj,
+            to_type=pandas_dict[scitype],
+            as_scitype=scitype,
+            store=store,
+            store_behaviour=store_behaviour,
+        )
+
+    def _init_conversion(self, X):
+        """Convert X to a pandas multiindex format."""
+        is_scitype = self.is_scitype
+        return self._coerce_to_df(X, is_scitype, store=self.converter_store)
+
     def _init_iter_indices(self):
         """Initialize indices that are iterated over in vectorization."""
         iterate_as = self.iterate_as
@@ -162,6 +160,11 @@ class VectorizedDF:
 
         return row_ix, col_ix
 
+    @property
+    def index(self):
+        """Defaults to pandas index of X converted to pandas type."""
+        return self.X_multiindex.index
+
     def get_iter_indices(self):
         """Get indices that are iterated over in vectorization.
 
@@ -248,7 +251,13 @@ class VectorizedDF:
         """Shorthand to retrieve self (iterator) as list."""
         return list(self)
 
-    def reconstruct(self, df_list, convert_back=False, overwrite_index=True):
+    def reconstruct(
+        self,
+        df_list,
+        convert_back=False,
+        overwrite_index=True,
+        col_multiindex="none",
+    ):
         """Reconstruct original format from iterable of vectorization instances.
 
         Parameters
@@ -256,41 +265,90 @@ class VectorizedDF:
         df_list : iterable of objects of same type and sequence as __getitem__ returns.
             can be self, but will in general be another object to be useful.
             Example: [some_operation(df) for df in self] that leaves types the same
-        convert_back : bool, default = False
+        convert_back : bool, optional, default = False
             whether to convert output back to mtype of X in __init__
             if False, the return will be a pandas.DataFrame with Index or multiIndex
             if True, the return is converted to the mtype of X in __init__
-        overwrite_index : bool, default = True
+        overwrite_index : bool, optional, default = True
             if True, the resulting return will have index overwritten by that of X
                 only if applies, i.e., overwrite is possible and X had an index
             if False, no index overwrite will happen
+        col_multiindex : str, one of "none", "flat", "multiindex", default = "none"
+            whether column multiindex is introduced, and if yes, what kind of,
+            in case of column vectorization being applied
+            "none" - df_list are simply concatenated, unless:
+                If at least one var is transformed to two or more, "flat" is enforced.
+                If this would cause non-unique column names, "flat" is enforced.
+            "multiindex" - additional level is introduced, by names of X passed to init,
+                if there would be more than one column underneath at least one level
+                this is added as an additional pandas MultiIndex level
+            "flat" - additional level is introduced, by names of X passed to init
+                if there would be more than one column underneath at least one level
+                like "multiindex", but new level is flattened to "var__colname" strings
+                no new multiindex level is introduced
 
         Returns
         -------
         X_reconstructed_orig_format : row-concatenation of df-list,
             with keys and additional level from self.get_iter_indices
             if convert_back=False, always a pd.DataFrame in an sktime MultiIndex format
-                (pd-multiindex mtype for Panel, or pd_multiindex_hier for Hierarchical)
+                (pd-muliindex mtype for Panel, or pd_multiindex_hier for Hierarchical)
             if convert_back=True, will have same format and mtype as X input to __init__
         """
+
+        def coerce_to_df(x):
+            if not isinstance(x, pd.DataFrame):
+                return self._coerce_to_df(x)
+            else:
+                return x
+
+        df_list = [coerce_to_df(x) for x in df_list]
+
+        # condition where "flat" behaviour is enforced if "none":
+        def _force_flat(df_list):
+            # transformer produces at least one multivariate output
+            # and there is more than one data frame to concatenate columns of
+            force_flat = len(df_list) > 1 and any(len(x.columns) > 1 for x in df_list)
+            # or, there would be duplicate columns, after the transformation
+            all_col_idx = [ix for df in df_list for ix in df.columns]
+            force_flat = force_flat or len(set(all_col_idx)) != len(all_col_idx)
+            return force_flat
+
         row_ix, col_ix = self.get_iter_indices()
+        force_flat = False
         if row_ix is None and col_ix is None:
             X_mi_reconstructed = self.X_multiindex
         elif col_ix is None:
             X_mi_reconstructed = pd.concat(df_list, keys=row_ix, axis=0)
         elif row_ix is None:
-            X_mi_reconstructed = pd.concat(df_list, axis=1)
+            force_flat = _force_flat(df_list)
+            if col_multiindex in ["flat", "multiindex"] or force_flat:
+                col_keys = self.X_multiindex.columns
+            else:
+                col_keys = None
+            X_mi_reconstructed = pd.concat(df_list, axis=1, keys=col_keys)
         else:  # both col_ix and row_ix are not None
             col_concats = []
             row_n = len(row_ix)
             col_n = len(col_ix)
             for i in range(row_n):
                 ith_col_block = df_list[i * col_n : (i + 1) * col_n]
-                col_concats += [pd.concat(ith_col_block, axis=1)]
+                force_flat = force_flat or _force_flat(ith_col_block)
+                if col_multiindex in ["flat", "multiindex"] or force_flat:
+                    col_keys = self.X_multiindex.columns
+                else:
+                    col_keys = None
+                col_concats += [pd.concat(ith_col_block, axis=1, keys=col_keys)]
+
             X_mi_reconstructed = pd.concat(col_concats, keys=row_ix, axis=0)
 
         X_mi_index = X_mi_reconstructed.index
         X_orig_row_index = self.X_multiindex.index
+
+        flatten = col_multiindex == "flat" or (col_multiindex == "none" and force_flat)
+        if flatten and isinstance(X_mi_reconstructed.columns, pd.MultiIndex):
+            X_mi_reconstructed.columns = flatten_multiindex(X_mi_reconstructed.columns)
+
         if overwrite_index and len(X_mi_index.names) == len(X_orig_row_index.names):
             X_mi_reconstructed.index = X_mi_index.set_names(X_orig_row_index.names)
 
diff --git a/sktime/forecasting/base/_fh.py b/sktime/forecasting/base/_fh.py
index f871dfdd0..8d9bc6712 100644
--- a/sktime/forecasting/base/_fh.py
+++ b/sktime/forecasting/base/_fh.py
@@ -643,7 +643,7 @@ def _to_relative(fh: ForecastingHorizon, cutoff=None) -> ForecastingHorizon:
         # Out: Index([<0 * Hours>, <4 * Hours>, <8 * Hours>], dtype = 'object')
         # [v - periods[0] for v in periods]
         # Out: Index([<0 * Hours>, <2 * Hours>, <4 * Hours>], dtype='object')
-        # TODO: v0.12.1 OR 0.13.0: Check if this comment below can be removed,
+        # TODO: 0.14.0: Check if this comment below can be removed,
         # so check if pandas has released the fix to PyPI:
         # This bug was reported: https://github.com/pandas-dev/pandas/issues/45999
         # and fixed: https://github.com/pandas-dev/pandas/pull/46006
diff --git a/sktime/forecasting/model_selection/_tune.py b/sktime/forecasting/model_selection/_tune.py
index 412efa115..c13e9aa9d 100644
--- a/sktime/forecasting/model_selection/_tune.py
+++ b/sktime/forecasting/model_selection/_tune.py
@@ -12,15 +12,13 @@ import numpy as np
 import pandas as pd
 from joblib import Parallel, delayed
 from sklearn.model_selection import ParameterGrid, ParameterSampler, check_cv
-from sklearn.utils.metaestimators import if_delegate_has_method
 
-from sktime.exceptions import NotFittedError
-from sktime.forecasting.base import BaseForecaster
+from sktime.forecasting.base._delegate import _DelegatedForecaster
 from sktime.forecasting.model_evaluation import evaluate
 from sktime.utils.validation.forecasting import check_scoring
 
 
-class BaseGridSearch(BaseForecaster):
+class BaseGridSearch(_DelegatedForecaster):
 
     _tags = {
         "scitype:y": "both",
@@ -68,108 +66,11 @@ class BaseGridSearch(BaseForecaster):
         ]
         self.clone_tags(forecaster, tags_to_clone)
 
-    @if_delegate_has_method(delegate=("best_forecaster_", "forecaster"))
-    def _update(self, y, X=None, update_params=False):
-        """Call _update on the forecaster with the best found parameters."""
-        self.check_is_fitted("update")
-        self.best_forecaster_._update(y, X, update_params=update_params)
-        return self
-
-    @if_delegate_has_method(delegate=("best_forecaster_", "forecaster"))
-    def _predict(self, fh=None, X=None):
-        """Call _predict on the forecaster with the best found parameters."""
-        self.check_is_fitted("predict")
-        return self.best_forecaster_._predict(fh=fh, X=X)
-
-    @if_delegate_has_method(delegate=("best_forecaster_", "forecaster"))
-    def _predict_interval(self, fh=None, X=None, coverage=None):
-        """Call _predict_interval on the forecaster with the best found parameters."""
-        self.check_is_fitted("predict")
-        return self.best_forecaster_._predict_interval(fh=fh, X=X, coverage=coverage)
-
-    @if_delegate_has_method(delegate=("best_forecaster_", "forecaster"))
-    def _predict_quantiles(self, fh=None, X=None, alpha=None):
-        """Call _predict_quantiles on the forecaster with the best found parameters."""
-        self.check_is_fitted("predict")
-        return self.best_forecaster_._predict_quantiles(fh=fh, X=X, alpha=alpha)
-
-    @if_delegate_has_method(delegate=("best_forecaster_", "forecaster"))
-    def _predict_var(self, fh=None, X=None, cov=False):
-        """Call _predict_var on the forecaster with the best found parameters."""
-        self.check_is_fitted("predict")
-        return self.best_forecaster_._predict_var(fh=fh, X=X, cov=cov)
-
-    @if_delegate_has_method(delegate=("best_forecaster_", "forecaster"))
-    def _predict_proba(self, fh=None, X=None, marginal=True):
-        """Call _predict_proba on the forecaster with the best found parameters."""
-        self.check_is_fitted("predict")
-        return self.best_forecaster_._predict_proba(fh=fh, X=X, marginal=marginal)
-
-    @if_delegate_has_method(delegate=("best_forecaster_", "forecaster"))
-    def transform(self, y, X=None):
-        """Call transform on the forecaster with the best found parameters."""
-        self.check_is_fitted("transform")
-        return self.best_forecaster_.transform(y, X)
-
-    @if_delegate_has_method(delegate=("best_forecaster_", "forecaster"))
-    def get_fitted_params(self):
-        """Get fitted parameters.
-
-        Returns
-        -------
-        fitted_params : dict
-        """
-        self.check_is_fitted("get_fitted_params")
-        return self.best_forecaster_.get_fitted_params()
-
-    @if_delegate_has_method(delegate=("best_forecaster_", "forecaster"))
-    def inverse_transform(self, y, X=None):
-        """Call inverse_transform on the forecaster with the best found params.
-
-        Only available if the underlying forecaster implements
-        ``inverse_transform`` and ``refit=True``.
-
-        Parameters
-        ----------
-        y : indexable, length n_samples
-            Must fulfill the input assumptions of the
-            underlying forecaster.
-        """
-        self.check_is_fitted("inverse_transform")
-        return self.best_forecaster_.inverse_transform(y, X)
+    _delegate_name = "best_forecaster_"
 
     def _run_search(self, evaluate_candidates):
         raise NotImplementedError("abstract method")
 
-    def check_is_fitted(self, method_name=None):
-        """Check if `fit` has been called.
-
-        Parameters
-        ----------
-        method_name : str
-            Name of the calling method.
-
-        Raises
-        ------
-        NotFittedError
-            If forecaster has not been fitted yet.
-        """
-        super(BaseGridSearch, self).check_is_fitted()
-
-        # We additionally check if the tuned forecaster has been fitted.
-        if method_name is not None:
-            if not self.refit:
-                raise NotFittedError(
-                    "This %s instance was initialized "
-                    "with refit=False. %s is "
-                    "available only after refitting on the "
-                    "best parameters. You can refit an forecaster "
-                    "manually using the ``best_params_`` "
-                    "attribute" % (type(self).__name__, method_name)
-                )
-            else:
-                self.best_forecaster_.check_is_fitted()
-
     def _fit(self, y, X=None, fh=None):
         """Fit to training data.
 
@@ -257,7 +158,7 @@ class BaseGridSearch(BaseForecaster):
 
         # Rank results, according to whether greater is better for the given scoring.
         results[f"rank_{scoring_name}"] = results.loc[:, f"mean_{scoring_name}"].rank(
-            ascending=not scoring.greater_is_better
+            ascending=scoring.get_tag("lower_is_better")
         )
 
         self.cv_results_ = results
@@ -274,7 +175,7 @@ class BaseGridSearch(BaseForecaster):
 
         # Sort values according to rank
         results = results.sort_values(
-            by=f"rank_{scoring_name}", ascending=not scoring.greater_is_better
+            by=f"rank_{scoring_name}", ascending=scoring.get_tag("lower_is_better")
         )
         # Select n best forecaster
         self.n_best_forecasters_ = []
@@ -499,6 +400,7 @@ class ForecastingGridSearchCV(BaseGridSearch):
         -------
         params : dict or list of dict
         """
+        from sktime.forecasting.exp_smoothing import ExponentialSmoothing
         from sktime.forecasting.model_selection._split import SingleWindowSplitter
         from sktime.forecasting.naive import NaiveForecaster
         from sktime.performance_metrics.forecasting import MeanAbsolutePercentageError
@@ -509,7 +411,13 @@ class ForecastingGridSearchCV(BaseGridSearch):
             "param_grid": {"window_length": [2, 5]},
             "scoring": MeanAbsolutePercentageError(symmetric=True),
         }
-        return params
+        params2 = {
+            "forecaster": ExponentialSmoothing(),
+            "cv": SingleWindowSplitter(fh=1),
+            "param_grid": {"initialization_method": ["estimated", "heuristic"]},
+            "scoring": MeanAbsolutePercentageError(symmetric=True),
+        }
+        return [params, params2]
 
 
 class ForecastingRandomizedSearchCV(BaseGridSearch):
diff --git a/sktime/forecasting/naive.py b/sktime/forecasting/naive.py
index 2beceff16..000d9984b 100644
--- a/sktime/forecasting/naive.py
+++ b/sktime/forecasting/naive.py
@@ -25,22 +25,79 @@ from sktime.datatypes._convert import convert, convert_to
 from sktime.datatypes._utilities import get_slice
 from sktime.forecasting.base._base import DEFAULT_ALPHA, BaseForecaster
 from sktime.forecasting.base._sktime import _BaseWindowForecaster
-from sktime.forecasting.compose import ColumnEnsembleForecaster
 from sktime.utils.validation import check_window_length
 from sktime.utils.validation.forecasting import check_sp
 
 
-class _NaiveForecaster(_BaseWindowForecaster):
-    """Univariate NaiveForecaster."""
+class NaiveForecaster(_BaseWindowForecaster):
+    """Forecast based on naive assumptions about past trends continuing.
+
+    NaiveForecaster is a forecaster that makes forecasts using simple
+    strategies. Two out of three strategies are robust against NaNs. The
+    NaiveForecaster can also be used for multivariate data and it then
+    applies internally the ColumnEnsembleForecaster, so each column
+    is forecasted with the same strategy.
+
+    Internally, this forecaster does the following:
+    - obtains the so-called "last window", a 1D array that denotes the
+      most recent time window that the forecaster is allowed to use
+    - reshapes the last window into a 2D array according to the given
+      seasonal periodicity (prepended with NaN values to make it fit);
+    - make a prediction for each column, using the given strategy:
+      - "last": last non-NaN row
+      - "mean": np.nanmean over rows
+    - tile the predictions using the seasonal periodicity
+
+    Parameters
+    ----------
+    strategy : {"last", "mean", "drift"}, default="last"
+        Strategy used to make forecasts:
+
+        * "last":   (robust against NaN values)
+                    forecast the last value in the
+                    training series when sp is 1.
+                    When sp is not 1,
+                    last value of each season
+                    in the last window will be
+                    forecasted for each season.
+        * "mean":   (robust against NaN values)
+                    forecast the mean of last window
+                    of training series when sp is 1.
+                    When sp is not 1, mean of all values
+                    in a season from last window will be
+                    forecasted for each season.
+        * "drift":  (not robust against NaN values)
+                    forecast by fitting a line between the
+                    first and last point of the window and
+                    extrapolating it into the future.
+
+    sp : int, default=1
+        Seasonal periodicity to use in the seasonal forecasting.
+
+    window_length : int or None, default=None
+        Window length to use in the `mean` strategy. If None, entire training
+            series will be used.
+
+    Examples
+    --------
+    >>> from sktime.datasets import load_airline
+    >>> from sktime.forecasting.naive import NaiveForecaster
+    >>> y = load_airline()
+    >>> forecaster = NaiveForecaster(strategy="drift")
+    >>> forecaster.fit(y)
+    NaiveForecaster(...)
+    >>> y_pred = forecaster.predict(fh=[1,2,3])
+    """
 
     _tags = {
+        "y_inner_mtype": "pd.Series",
         "requires-fh-in-fit": False,
-        "handles-missing-data": True,  # todo: switch to True if GH1367 is fixed
+        "handles-missing-data": True,
         "scitype:y": "univariate",
     }
 
     def __init__(self, strategy="last", window_length=None, sp=1):
-        super(_NaiveForecaster, self).__init__()
+        super(NaiveForecaster, self).__init__()
         self.strategy = strategy
         self.sp = sp
         self.window_length = window_length
@@ -247,111 +304,6 @@ class _NaiveForecaster(_BaseWindowForecaster):
         fh_idx = fh.to_indexer(self.cutoff)
         return y_pred[fh_idx]
 
-
-class NaiveForecaster(BaseForecaster):
-    """Forecast based on naive assumptions about past trends continuing.
-
-    NaiveForecaster is a forecaster that makes forecasts using simple
-    strategies. Two out of three strategies are robust against NaNs. The
-    NaiveForecaster can also be used for multivariate data and it then
-    applies internally the ColumnEnsembleForecaster, so each column
-    is forecasted with the same strategy.
-
-    Internally, this forecaster does the following:
-    - obtains the so-called "last window", a 1D array that denotes the
-      most recent time window that the forecaster is allowed to use
-    - reshapes the last window into a 2D array according to the given
-      seasonal periodicity (prepended with NaN values to make it fit);
-    - make a prediction for each column, using the given strategy:
-      - "last": last non-NaN row
-      - "mean": np.nanmean over rows
-    - tile the predictions using the seasonal periodicity
-
-    Parameters
-    ----------
-    strategy : {"last", "mean", "drift"}, default="last"
-        Strategy used to make forecasts:
-
-        * "last":   (robust against NaN values)
-                    forecast the last value in the
-                    training series when sp is 1.
-                    When sp is not 1,
-                    last value of each season
-                    in the last window will be
-                    forecasted for each season.
-        * "mean":   (robust against NaN values)
-                    forecast the mean of last window
-                    of training series when sp is 1.
-                    When sp is not 1, mean of all values
-                    in a season from last window will be
-                    forecasted for each season.
-        * "drift":  (not robust against NaN values)
-                    forecast by fitting a line between the
-                    first and last point of the window and
-                    extrapolating it into the future.
-
-    sp : int, default=1
-        Seasonal periodicity to use in the seasonal forecasting.
-
-    window_length : int or None, default=None
-        Window length to use in the `mean` strategy. If None, entire training
-            series will be used.
-
-    Examples
-    --------
-    >>> from sktime.datasets import load_airline
-    >>> from sktime.forecasting.naive import NaiveForecaster
-    >>> y = load_airline()
-    >>> forecaster = NaiveForecaster(strategy="drift")
-    >>> forecaster.fit(y)
-    NaiveForecaster(...)
-    >>> y_pred = forecaster.predict(fh=[1,2,3])
-    """
-
-    _tags = {
-        "y_inner_mtype": ["pd.DataFrame", "pd.Series"],
-        "scitype:y": "both",
-        "requires-fh-in-fit": False,
-        "handles-missing-data": True,  # todo: switch to True if GH1367 is fixed
-    }
-
-    def __init__(self, strategy="last", window_length=None, sp=1):
-        self.strategy = strategy
-        self.sp = sp
-        self.window_length = window_length
-        super(NaiveForecaster, self).__init__()
-
-    def _fit(self, y, X=None, fh=None):
-        """Fit to training data.
-
-        Parameters
-        ----------
-        y : pd.Series, pd.DataFrame
-            Target time series to which to fit the forecaster.
-        fh : int, list or np.array, default=None
-            The forecasters horizon with the steps ahead to to predict.
-        X : pd.DataFrame, default=None
-            Exogenous variables are ignored.
-
-        Returns
-        -------
-        self : returns an instance of self.
-        """
-        self._forecaster = ColumnEnsembleForecaster(
-            _NaiveForecaster(
-                strategy=self.strategy, sp=self.sp, window_length=self.window_length
-            )
-        )
-
-        # 1st part of preserving name/columns, otherwise they get lost on occasion
-        # provisionally hard-fixing the mess of 100 nested function calls
-        # until someone finds it in their heart to simplify the NaiveForecaster
-        if isinstance(y, pd.Series):
-            self.cols_ = y.name
-        else:
-            self.cols_ = y.columns
-        self._forecaster.fit(y=y, X=X, fh=fh)
-
     def _predict(self, fh=None, X=None):
         """Forecast time series at future horizon.
 
@@ -362,41 +314,23 @@ class NaiveForecaster(BaseForecaster):
         X : pd.DataFrame, optional (default=None)
             Exogenous time series
         """
-        y_pred = self._forecaster.predict(fh=fh, X=X)
+        y_pred = super(NaiveForecaster, self)._predict(fh=fh, X=X)
+
+        # test_predict_time_index_in_sample_full[ForecastingPipeline-0-int-int-True]
+        #   causes a pd.DataFrame to appear as y_pred, which upsets the next lines
+        #   reasons are unclear, this is coming from the _BaseWindowForecaster
+        # todo: investigate this
+        if isinstance(y_pred, pd.DataFrame):
+            y_pred = y_pred.iloc[:, 0]
 
         # check for in-sample prediction, if first time point needs to be imputed
         if self._y.index[0] in y_pred.index:
-            # fill NaN with next row values
-            y_pred.loc[self._y.index[0]] = y_pred.loc[self._y.index[1]]
-
-        # 2nd part of preserving name/columns, otherwise they get lost on occasion
-        # provisionally hard-fixing the mess of 100 nested function calls
-        # until someone finds it in their heart to simplify the NaiveForecaster
-        if isinstance(y_pred, pd.Series):
-            y_pred.name = self.cols_
-        else:  # is pd.DataFrame
-            y_pred.columns = self.cols_
+            if y_pred.loc[[self._y.index[0]]].hasnans:
+                # fill NaN with observed values
+                y_pred.loc[self._y.index[0]] = self._y[self._y.index[1]]
 
         return y_pred
 
-    def _update(self, y, X=None, update_params=True):
-        """Update cutoff value and, optionally, fitted parameters.
-
-        Parameters
-        ----------
-        y : pd.Series, pd.DataFrame, or np.array
-            Target time series to which to fit the forecaster.
-        X : pd.DataFrame, optional (default=None)
-            Exogeneous data
-        update_params : bool, optional (default=True)
-            whether model parameters should be updated
-
-        Returns
-        -------
-        self : reference to self
-        """
-        return self._forecaster.update(y=y, X=X, update_params=update_params)
-
     @classmethod
     def get_test_params(cls, parameter_set="default"):
         """Return testing parameter settings for the estimator.
@@ -417,9 +351,9 @@ class NaiveForecaster(BaseForecaster):
         """
         params_list = [
             {},
-            {"sp": 2},
-            {"strategy": "mean"},
+            {"strategy": "mean", "sp": 2},
             {"strategy": "drift"},
+            {"strategy": "last"},
             {"strategy": "mean", "window_length": 5},
         ]
 
diff --git a/sktime/forecasting/varmax.py b/sktime/forecasting/varmax.py
new file mode 100644
index 000000000..56b30b668
--- /dev/null
+++ b/sktime/forecasting/varmax.py
@@ -0,0 +1,405 @@
+# -*- coding: utf-8 -*-
+"""Vector Autoregressive Moving Average with eXogenous regressors model (VARMAX)."""
+__all__ = ["VARMAX"]
+__author__ = ["KatieBuc"]
+
+import warnings
+
+import pandas as pd
+from statsmodels.tsa.statespace.varmax import VARMAX as _VARMAX
+
+from sktime.forecasting.base.adapters import _StatsModelsAdapter
+
+
+class VARMAX(_StatsModelsAdapter):
+    r"""Wrapper for statsmodels VARMAX model.
+
+    Vector Autoregressive Moving Average with eXogenous regressors model (VARMAX)
+
+    Parameters
+    ----------
+    order : iterable
+        The (p,q) order of the model for the number of AR and MA parameters to
+        use.
+    trend : str{'n','c','t','ct'} or iterable, optional
+        Parameter controlling the deterministic trend polynomial :math:`A(t)`.
+        Can be specified as a string where 'c' indicates a constant (i.e. a
+        degree zero component of the trend polynomial), 't' indicates a
+        linear trend with time, and 'ct' is both. Can also be specified as an
+        iterable defining the non-zero polynomial exponents to include, in
+        increasing order. For example, `[1,1,0,1]` denotes
+        :math:`a + bt + ct^3`. Default is a constant trend component.
+    error_cov_type : {'diagonal', 'unstructured'}, optional
+        The structure of the covariance matrix of the error term, where
+        "unstructured" puts no restrictions on the matrix and "diagonal"
+        requires it to be a diagonal matrix (uncorrelated errors). Default is
+        "unstructured".
+    measurement_error : bool, optional
+        Whether or not to assume the endogenous observations `endog` were
+        measured with error. Default is False.
+    enforce_stationarity : bool, optional
+        Whether or not to transform the AR parameters to enforce stationarity
+        in the autoregressive component of the model. Default is True.
+    enforce_invertibility : bool, optional
+        Whether or not to transform the MA parameters to enforce invertibility
+        in the moving average component of the model. Default is True.
+    trend_offset : int, optional
+        The offset at which to start time trend values. Default is 1, so that
+        if `trend='t'` the trend is equal to 1, 2, ..., n_obs. Typically is only
+        set when the model created by extending a previous dataset.
+    start_params : array_like, optional
+        Initial guess of the solution for the loglikelihood maximization. If
+        None, the default is given by Model.start_params.
+    transformed : bool, optional
+        Whether or not start_params is already transformed. Default is True.
+    includes_fixed : bool, optional
+        If parameters were previously fixed with the fix_params method, this
+        argument describes whether or not start_params also includes the
+        fixed parameters, in addition to the free parameters. Default is False.
+    cov_type : str, optional
+        The `cov_type` keyword governs the method for calculating the
+        covariance matrix of parameter estimates. Can be one of:
+         - 'opg' for the outer product of gradient estimator
+         - 'oim' for the observed information matrix estimator, calculated
+            using the method of Harvey (1989)
+         - 'approx' for the observed information matrix estimator,
+            calculated using a numerical approximation of the Hessian matrix.
+         - 'robust' for an approximate (quasi-maximum likelihood) covariance
+            matrix that may be valid even in the presence of some
+            misspecifications. Intermediate calculations use the 'oim'
+            method.
+         - 'robust_approx' is the same as 'robust' except that the
+            intermediate calculations use the 'approx' method.
+         - 'none' for no covariance matrix calculation.
+        Default is 'opg' unless memory conservation is used to avoid computing the
+        loglikelihood values for each observation, in which case the default is
+        'approx'.
+    cov_kwds : dict or None, optional
+        A dictionary of arguments affecting covariance matrix computation.
+        opg, oim, approx, robust, robust_approx
+         - 'approx_complex_step' : bool, optional - If True, numerical
+            approximations are computed using complex-step methods. If False,
+            numerical approximations are computed using finite difference
+            methods. Default is True.
+         - 'approx_centered' : bool, optional - If True, numerical
+            approximations computed using finite difference methods use a
+            centered approximation. Default is False.
+    method : str, optional
+        The `method` determines which solver from `scipy.optimize`
+        is used, and it can be chosen from among the following strings:
+         - 'newton' for Newton-Raphson
+         - 'nm' for Nelder-Mead
+         - 'bfgs' for Broyden-Fletcher-Goldfarb-Shanno (BFGS)
+         - 'lbfgs' for limited-memory BFGS with optional box constraints
+         - 'powell' for modified Powell's method
+         - 'cg' for conjugate gradient
+         - 'ncg' for Newton-conjugate gradient
+         - 'basinhopping' for global basin-hopping solver
+        The explicit arguments in `fit` are passed to the solver,
+        with the exception of the basin-hopping solver. Each
+        solver has several optional arguments that are not the same across
+        solvers. See the notes section below (or scipy.optimize) for the
+        available arguments and for the list of explicit arguments that the
+        basin-hopping solver supports.
+    maxiter : int, optional
+        The maximum number of iterations to perform.
+    full_output : bool, optional
+        Set to True to have all available output in the Results object's
+        mle_retvals attribute. The output is dependent on the solver.
+        See LikelihoodModelResults notes section for more information.
+    disp : bool, optional
+        Set to True to print convergence messages.
+    callback : callable callback(xk), optional
+        Called after each iteration, as callback(xk), where xk is the
+        current parameter vector.
+    return_params : bool, optional
+        Whether or not to return only the array of maximizing parameters.
+        Default is False.
+    optim_score : {'harvey', 'approx'} or None, optional
+        The method by which the score vector is calculated. 'harvey' uses
+        the method from Harvey (1989), 'approx' uses either finite
+        difference or complex step differentiation depending upon the
+        value of `optim_complex_step`, and None uses the built-in gradient
+        approximation of the optimizer. Default is None. This keyword is
+        only relevant if the optimization method uses the score.
+    optim_complex_step : bool, optional
+        Whether or not to use complex step differentiation when
+        approximating the score; if False, finite difference approximation
+        is used. Default is True. This keyword is only relevant if
+        `optim_score` is set to 'harvey' or 'approx'.
+    optim_hessian : {'opg','oim','approx'}, optional
+        The method by which the Hessian is numerically approximated. 'opg'
+        uses outer product of gradients, 'oim' uses the information
+        matrix formula from Harvey (1989), and 'approx' uses numerical
+        approximation. This keyword is only relevant if the
+        optimization method uses the Hessian matrix.
+    low_memory : bool, optional
+        If set to True, techniques are applied to substantially reduce
+        memory usage. If used, some features of the results object will
+        not be available (including smoothed results and in-sample
+        prediction), although out-of-sample forecasting is possible.
+        Default is False.
+    dynamic : bool, int, str, or datetime, optional
+        Integer offset relative to `start` at which to begin dynamic
+        prediction. Can also be an absolute date string to parse or a
+        datetime type (these are not interpreted as offsets).
+        Prior to this observation, true endogenous values will be used for
+        prediction; starting with this observation and continuing through
+        the end of prediction, forecasted endogenous values will be used
+        instead.
+    information_set : str, optional
+        The information set to condition each prediction on. Default is
+        "predicted", which computes predictions of period t values
+        conditional on observed data through period t-1; these are
+        one-step-ahead predictions, and correspond with the typical
+        `fittedvalues` results attribute. Alternatives are "filtered",
+        which computes predictions of period t values conditional on
+        observed data through period t, and "smoothed", which computes
+        predictions of period t values conditional on the entire dataset
+        (including also future observations t+1, t+2, ...).
+    signal_only : bool, optional
+        Whether to compute predictions of only the "signal" component of
+        the observation equation. Default is False. For example, the
+        observation equation of a time-invariant model is
+        :math:`y_t = d + Z \alpha_t + \varepsilon_t`, and the "signal"
+        component is then :math:`Z \alpha_t`. If this argument is set to
+        True, then predictions of the "signal" :math:`Z \alpha_t` will be
+        returned. Otherwise, the default is for predictions of :math:`y_t`
+        to be returned.
+    suppress_warnings : bool, optional
+        Many warnings might be thrown inside of statsmodels. If
+        ``suppress_warnings`` is True, all of these warnings will be squelched.
+        Default is False.
+
+    Notes
+    -----
+    Generically, the VARMAX model is specified (see for example chapter 18 of
+    [1]_):
+    .. math::
+        y_t = A(t) + A_1 y_{t-1} + \dots + A_p y_{t-p} + B x_t + \epsilon_t +
+        M_1 \epsilon_{t-1} + \dots M_q \epsilon_{t-q}
+    where :math:`\epsilon_t \sim N(0, \Omega)`, and where :math:`y_t` is a
+    `k_endog x 1` vector. Additionally, this model allows considering the case
+    where the variables are measured with error.
+    Note that in the full VARMA(p,q) case there is a fundamental identification
+    problem in that the coefficient matrices :math:`\{A_i, M_j\}` are not
+    generally unique, meaning that for a given time series process there may
+    be multiple sets of matrices that equivalently represent it. See Chapter 12
+    of [1]_ for more information. Although this class can be used to estimate
+    VARMA(p,q) models, a warning is issued to remind users that no steps have
+    been taken to ensure identification in this case.
+
+    References
+    ----------
+    .. [1] Lütkepohl, Helmut. 2007.
+       New Introduction to Multiple Time Series Analysis.
+       Berlin: Springer.
+
+    Examples
+    --------
+    >>> from sktime.forecasting.varmax import VARMAX
+    >>> from sktime.datasets import load_macroeconomic
+    >>> from sktime.forecasting.model_selection import temporal_train_test_split
+    >>> y = load_macroeconomic()
+    >>> forecaster = VARMAX(suppress_warnings=True)
+    >>> forecaster.fit(y[['realgdp', 'unemp']])
+    VARMAX(...)
+    >>> y_pred = forecaster.predict(fh=[1,4,12])
+    """
+
+    _tags = {
+        "scitype:y": "multivariate",
+        "ignores-exogeneous-X": False,
+        "handles-missing-data": False,
+        "y_inner_mtype": "pd.DataFrame",
+        "X_inner_mtype": "pd.DataFrame",
+        "requires-fh-in-fit": False,
+        "X-y-must-have-same-index": True,
+        "enforce_index_type": None,
+        "capability:pred_int": False,
+    }
+
+    def __init__(
+        self,
+        order=(1, 0),
+        trend="c",
+        error_cov_type="unstructured",
+        measurement_error=False,
+        enforce_stationarity=True,
+        enforce_invertibility=True,
+        trend_offset=1,
+        start_params=None,
+        transformed=True,
+        includes_fixed=False,
+        cov_type=None,
+        cov_kwds=None,
+        method="lbfgs",
+        maxiter=50,
+        full_output=1,
+        disp=False,
+        callback=None,
+        return_params=False,
+        optim_score=None,
+        optim_complex_step=None,
+        optim_hessian=None,
+        flags=None,
+        low_memory=False,
+        dynamic=False,
+        information_set="predicted",
+        signal_only=False,
+        suppress_warnings=False,
+    ):
+        # Model parameters
+        self.order = order
+        self.trend = trend
+        self.error_cov_type = error_cov_type
+        self.measurement_error = measurement_error
+        self.enforce_stationarity = enforce_stationarity
+        self.enforce_invertibility = enforce_invertibility
+        self.trend_offset = trend_offset
+        self.start_params = start_params
+        self.transformed = transformed
+        self.includes_fixed = includes_fixed
+        self.cov_type = cov_type
+        self.cov_kwds = cov_kwds
+        self.method = method
+        self.maxiter = maxiter
+        self.full_output = full_output
+        self.disp = disp
+        self.callback = callback
+        self.return_params = return_params
+        self.optim_score = optim_score
+        self.optim_complex_step = optim_complex_step
+        self.optim_hessian = optim_hessian
+        self.flags = flags
+        self.low_memory = low_memory
+        self.dynamic = dynamic
+        self.information_set = information_set
+        self.signal_only = signal_only
+        self.suppress_warnings = suppress_warnings
+
+        super(VARMAX, self).__init__()
+
+    def _fit_forecaster(self, y, X=None):
+        """Fit forecaster to training data.
+
+        Writes to self:
+            Sets fitted model attributes ending in "_".
+
+        Parameters
+        ----------
+        y : array_like
+            The observed time-series process :math:`y`, shaped n_obs x k_endog.
+        X : array_like, optional (default=None)
+            Array of exogenous regressors, shaped n_obs x k.
+
+        Returns
+        -------
+        self : reference to self
+        """
+        if self.suppress_warnings:
+            warnings.filterwarnings("ignore")
+
+        self._forecaster = _VARMAX(
+            endog=y,
+            exog=X,
+            order=self.order,
+            trend=self.trend,
+            error_cov_type=self.error_cov_type,
+            measurement_error=self.measurement_error,
+            enforce_stationarity=self.enforce_stationarity,
+            enforce_invertibility=self.enforce_invertibility,
+            trend_offset=self.trend_offset,
+        )
+        self._fitted_forecaster = self._forecaster.fit(
+            start_params=self.start_params,
+            transformed=self.transformed,
+            includes_fixed=self.includes_fixed,
+            cov_type=self.cov_type,
+            cov_kwds=self.cov_kwds,
+            method=self.method,
+            maxiter=self.maxiter,
+            full_output=self.full_output,
+            disp=self.disp,
+            callback=self.callback,
+            return_params=self.return_params,
+            optim_score=self.optim_score,
+            optim_complex_step=self.optim_complex_step,
+            optim_hessian=self.optim_hessian,
+            flags=self.flags,
+            low_memory=self.low_memory,
+        )
+        return self
+
+    # defining `_predict`, instead of inheriting from `_StatsModelsAdapter`,
+    # for two reasons:
+    # 1. to pass in `dynamic`, `information_set` and `signal_only`
+    # 2. to deal with statsmodel integer indexing issue
+    def _predict(self, fh, X=None):
+        """
+        Wrap Statmodel's VARMAX forecast method.
+
+        Parameters
+        ----------
+        fh : ForecastingHorizon
+            The forecasters horizon with the steps ahead to to predict.
+            Default is one-step ahead forecast,
+            i.e. np.array([1])
+        X : pd.DataFrame, optional (default=None)
+            Exogenous variables.
+
+        Returns
+        -------
+        y_pred : np.ndarray
+            Returns series of predicted values.
+        """
+        start, end = fh.to_absolute_int(self._y.index[0], self.cutoff)[[0, -1]]
+
+        y_pred = self._fitted_forecaster.predict(
+            start=start,
+            end=end,
+            dynamic=self.dynamic,
+            information_set=self.information_set,
+            signal_only=self.signal_only,
+            exog=X,
+        )
+
+        # statsmodel returns zero-based index when index is of type int with the
+        # following warning
+        # ValueWarning: No supported index is available. Prediction results will be
+        # given with an integer index beginning at `start`...
+        # but only when out-of-sample forecasting, i.e. when forecasting horizon is
+        # greater than zero
+        if (type(self._y.index) == pd.core.indexes.numeric.Int64Index) & (
+            any(fh.to_relative(self.cutoff) > 0)
+        ):
+            y_pred.index = y_pred.index + self._y.index[0]
+
+        return y_pred.loc[fh.to_absolute(self.cutoff).to_pandas()]
+
+    @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 forecasters.
+
+        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 = [
+            {"order": (1, 0)},
+            {"order": (0, 1)},
+            {"order": (1, 1)},
+        ]
+
+        return params
diff --git a/sktime/performance_metrics/base/_base.py b/sktime/performance_metrics/base/_base.py
index 594e01c88..401fffb8f 100644
--- a/sktime/performance_metrics/base/_base.py
+++ b/sktime/performance_metrics/base/_base.py
@@ -5,8 +5,6 @@
 __author__ = ["rnkuhns", "fkiraly"]
 __all__ = ["BaseMetric"]
 
-from warnings import warn
-
 from sktime.base import BaseObject
 
 
@@ -16,22 +14,9 @@ class BaseMetric(BaseObject):
     Extends sktime BaseObject.
     """
 
-    # todo: 0.13.0, remove the func/name args here and to all the metrics
-    def __init__(
-        self,
-        func=None,
-        name=None,
-    ):
-
-        if func is not None or name is not None:
-            warn(
-                "func and name parameters of BaseMetric are deprecated from 0.12.0"
-                "and will be removed in 0.13.0",
-                DeprecationWarning,
-            )
+    def __init__(self):
 
-        self.func = func
-        self.name = name
+        super(BaseMetric, self).__init__()
 
     def __call__(self, y_true, y_pred, **kwargs):
         """Calculate metric value using underlying metric function.
diff --git a/sktime/performance_metrics/forecasting/_classes.py b/sktime/performance_metrics/forecasting/_classes.py
index f4f814dba..3eea4e338 100644
--- a/sktime/performance_metrics/forecasting/_classes.py
+++ b/sktime/performance_metrics/forecasting/_classes.py
@@ -65,21 +65,6 @@ __all__ = [
 ]
 
 
-class _BaseForecastingErrorMetric(BaseMetric):
-    """Base class for defining forecasting error metrics in sktime.
-
-    Extends sktime's BaseMetric to the forecasting interface. Forecasting error
-    metrics measure the error (loss) between forecasts and true values. Lower
-    values are better.
-    """
-
-    greater_is_better = False
-
-    def __init__(self, func, name=None, multioutput="uniform_average"):
-        self.multioutput = multioutput
-        super().__init__(func, name=name)
-
-
 class BaseForecastingErrorMetric(BaseMetric):
     """Base class for defining forecasting error metrics in sktime.
 
@@ -113,30 +98,9 @@ class BaseForecastingErrorMetric(BaseMetric):
         "inner_implements_multilevel": False,
     }
 
-    # todo: 0.13.0, remove the greater_is_better property
-    @property
-    def greater_is_better(self):
-        """Whether greater is better for the metric."""
-        warn(
-            "The greater_is_better attribute is deprecated from 0.12.0 "
-            "and will be removed in 0.13.0. Use the lower_is_better tag instead, "
-            'e.g., my_metric.get_tag("lower_is_better")',
-            DeprecationWarning,
-        )
-        return not self.get_tag("lower_is_better", False)
-
-    def __init__(
-        self,
-        func=None,
-        name=None,
-        multioutput="uniform_average",
-        multilevel="uniform_average",
-    ):
+    def __init__(self, multioutput="uniform_average", multilevel="uniform_average"):
         self.multioutput = multioutput
         self.multilevel = multilevel
-        # todo: in 0.13.0, when func/name are removed, remove func/name args
-        #   temporarily, this is how we avoid overwrite with None
-        # super().__init__(func=func, name=name)
         self.name = type(self).__name__
 
     def __call__(self, y_true, y_pred, **kwargs):
@@ -441,15 +405,25 @@ class BaseForecastingErrorMetricFunc(BaseForecastingErrorMetric):
         return res
 
 
-class _BaseForecastingScoreMetric(_BaseForecastingErrorMetric):
-    """Base class for defining forecasting score metrics in sktime.
+class _DynamicForecastingErrorMetric(BaseForecastingErrorMetricFunc):
+    """Class for defining forecasting error metrics from a function dynamically."""
 
-    Extends sktime's BaseMetric to the forecasting interface. Forecasting score
-    metrics measure the agreement between forecasts and true values. Higher
-    values are better.
-    """
+    def __init__(
+        self,
+        func,
+        name=None,
+        multioutput="uniform_average",
+        multilevel="uniform_average",
+        lower_is_better=True,
+    ):
+        self.multioutput = multioutput
+        self.multilevel = multilevel
+        self.func = func
+        self.name = name
+        super().__init__()
 
-    _tags = {"lower_is_better": False}
+        if not lower_is_better:
+            self.set_tags(**{"lower_is_better": False})
 
 
 class _ScaledMetricTags:
@@ -463,9 +437,13 @@ class _ScaledMetricTags:
 
 
 def make_forecasting_scorer(
-    func, name=None, greater_is_better=False, multioutput="uniform_average"
+    func,
+    name=None,
+    greater_is_better=False,
+    multioutput="uniform_average",
+    multilevel="uniform_average",
 ):
-    """Create a metric class from metric functions.
+    """Create a metric class from a metric functions.
 
     Parameters
     ----------
@@ -483,16 +461,24 @@ def make_forecasting_scorer(
         If array-like, values used as weights to average the errors.
         If 'raw_values', returns a full set of errors in case of multioutput input.
         If 'uniform_average', errors of all outputs are averaged with uniform weight.
+    multilevel : {'raw_values', 'uniform_average'}
+        Defines how to aggregate metric for hierarchical data (with levels).
+        If 'uniform_average' (default), errors are mean-averaged across levels.
+        If 'raw_values', does not average errors across levels, hierarchy is retained.
 
     Returns
     -------
     scorer:
         Metric class that can be used as forecasting scorer.
     """
-    if greater_is_better:
-        return _BaseForecastingScoreMetric(func, name=name, multioutput=multioutput)
-    else:
-        return _BaseForecastingErrorMetric(func, name=name, multioutput=multioutput)
+    lower_is_better = not greater_is_better
+    return _DynamicForecastingErrorMetric(
+        func,
+        name=name,
+        multioutput=multioutput,
+        multilevel=multilevel,
+        lower_is_better=lower_is_better,
+    )
 
 
 class MeanAbsoluteScaledError(_ScaledMetricTags, BaseForecastingErrorMetricFunc):
diff --git a/sktime/performance_metrics/forecasting/probabilistic/_classes.py b/sktime/performance_metrics/forecasting/probabilistic/_classes.py
index d935874ab..210a1566d 100644
--- a/sktime/performance_metrics/forecasting/probabilistic/_classes.py
+++ b/sktime/performance_metrics/forecasting/probabilistic/_classes.py
@@ -8,12 +8,12 @@ import pandas as pd
 from sklearn.utils import check_array, check_consistent_length
 
 from sktime.datatypes import check_is_scitype, convert
-from sktime.performance_metrics.forecasting._classes import _BaseForecastingErrorMetric
+from sktime.performance_metrics.forecasting._classes import BaseForecastingErrorMetric
 
 # TODO: Rework tests now
 
 
-class _BaseProbaForecastingErrorMetric(_BaseForecastingErrorMetric):
+class _BaseProbaForecastingErrorMetric(BaseForecastingErrorMetric):
     """Base class for probabilistic forecasting error metrics in sktime.
 
     Extends sktime's BaseMetric to the forecasting interface. Forecasting error
@@ -39,16 +39,10 @@ class _BaseProbaForecastingErrorMetric(_BaseForecastingErrorMetric):
         "lower_is_better": True,
     }
 
-    def __init__(
-        self,
-        func=None,
-        name=None,
-        multioutput="uniform_average",
-        score_average=True,
-    ):
+    def __init__(self, multioutput="uniform_average", score_average=True):
         self.multioutput = multioutput
         self.score_average = score_average
-        super().__init__(func, name=name, multioutput=multioutput)
+        super().__init__(multioutput=multioutput)
 
     def __call__(self, y_true, y_pred, **kwargs):
         """Calculate metric value using underlying metric function.
@@ -398,19 +392,11 @@ class PinballLoss(_BaseProbaForecastingErrorMetric):
         "lower_is_better": True,
     }
 
-    def __init__(
-        self,
-        multioutput="uniform_average",
-        score_average=True,
-        alpha=None,
-    ):
-        name = "PinballLoss"
+    def __init__(self, multioutput="uniform_average", score_average=True, alpha=None):
         self.score_average = score_average
         self.alpha = self._check_alpha(alpha)
         self.metric_args = {"alpha": alpha}
-        super().__init__(
-            name=name, multioutput=multioutput, score_average=score_average
-        )
+        super().__init__(multioutput=multioutput, score_average=score_average)
 
     def _evaluate_by_index(self, y_true, y_pred, multioutput, **kwargs):
         """Logic for finding the metric evaluated at each index.
@@ -489,12 +475,9 @@ class EmpiricalCoverage(_BaseProbaForecastingErrorMetric):
     }
 
     def __init__(self, multioutput="uniform_average", score_average=True):
-        name = "EmpiricalCoverage"
         self.score_average = score_average
         self.multioutput = multioutput
-        super().__init__(
-            name=name, score_average=score_average, multioutput=multioutput
-        )
+        super().__init__(score_average=score_average, multioutput=multioutput)
 
     def _evaluate_by_index(self, y_true, y_pred, multioutput, **kwargs):
         """Logic for finding the metric evaluated at each index.
@@ -559,12 +542,9 @@ class ConstraintViolation(_BaseProbaForecastingErrorMetric):
     }
 
     def __init__(self, multioutput="uniform_average", score_average=True):
-        name = "ConstraintViolation"
         self.score_average = score_average
         self.multioutput = multioutput
-        super().__init__(
-            name=name, score_average=score_average, multioutput=multioutput
-        )
+        super().__init__(score_average=score_average, multioutput=multioutput)
 
     def _evaluate_by_index(self, y_true, y_pred, multioutput, **kwargs):
         """Logic for finding the metric evaluated at each index.
diff --git a/sktime/transformations/bootstrap/_mbb.py b/sktime/transformations/bootstrap/_mbb.py
index ec5e9313b..9ee8b91d5 100644
--- a/sktime/transformations/bootstrap/_mbb.py
+++ b/sktime/transformations/bootstrap/_mbb.py
@@ -174,7 +174,7 @@ class STLBootstrapTransformer(BaseTransformer):
         # todo: what is the scitype of y: None (not needed), Primitives, Series, Panel
         "scitype:transform-labels": "None",
         "scitype:instancewise": True,  # is this an instance-wise transform?
-        "X_inner_mtype": "pd.Series",  # which mtypes do _fit/_predict support for X?
+        "X_inner_mtype": "pd.DataFrame",  # which mtypes do _fit/_predict support for X?
         # X_inner_mtype can be Panel mtype even if transform-input is Series, vectorized
         "y_inner_mtype": "None",  # which mtypes do _fit/_predict support for y?
         "capability:inverse_transform": False,
@@ -293,6 +293,9 @@ class STLBootstrapTransformer(BaseTransformer):
         -------
         transformed version of X
         """
+        Xcol = X.columns
+        X = X[X.columns[0]]
+
         if len(X) <= self.block_length_:
             raise ValueError(
                 "STLBootstrapTransformer requires that block_length is"
@@ -366,7 +369,10 @@ class STLBootstrapTransformer(BaseTransformer):
                 )
             )
 
-        return pd.concat(df_list)
+        Xt = pd.concat(df_list)
+        Xt.columns = Xcol
+
+        return Xt
 
     @classmethod
     def get_test_params(cls, parameter_set="default"):
@@ -479,7 +485,7 @@ class MovingBlockBootstrapTransformer(BaseTransformer):
         # todo: what is the scitype of y: None (not needed), Primitives, Series, Panel
         "scitype:transform-labels": "None",
         "scitype:instancewise": True,  # is this an instance-wise transform?
-        "X_inner_mtype": "pd.Series",  # which mtypes do _fit/_predict support for X?
+        "X_inner_mtype": "pd.DataFrame",  # which mtypes do _fit/_predict support for X?
         # X_inner_mtype can be Panel mtype even if transform-input is Series, vectorized
         "y_inner_mtype": "None",  # which mtypes do _fit/_predict support for y?
         "capability:inverse_transform": False,
@@ -525,6 +531,9 @@ class MovingBlockBootstrapTransformer(BaseTransformer):
         -------
         transformed version of X
         """
+        Xcol = X.columns
+        X = X[X.columns[0]]
+
         if len(X) <= self.block_length:
             raise ValueError(
                 "MovingBlockBootstrapTransformer requires that block_length is"
@@ -573,7 +582,10 @@ class MovingBlockBootstrapTransformer(BaseTransformer):
                 )
             )
 
-        return pd.concat(df_list)
+        Xt = pd.concat(df_list)
+        Xt.columns = Xcol
+
+        return Xt
 
     @classmethod
     def get_test_params(cls, parameter_set="default"):
diff --git a/sktime/transformations/panel/segment.py b/sktime/transformations/panel/segment.py
index 5931ebd1d..cf6c9e928 100644
--- a/sktime/transformations/panel/segment.py
+++ b/sktime/transformations/panel/segment.py
@@ -125,9 +125,10 @@ class IntervalSegmenter(BaseTransformer):
         new_column_names = []
         for interval in self.intervals_:
             start, end = interval[0], interval[-1]
-            interval = X[:, start:end]
-            intervals.append(interval)
-            new_column_names.append(f"{column_names}_{start}_{end}")
+            if f"{column_names}_{start}_{end}" not in new_column_names:
+                interval = X[:, start:end]
+                intervals.append(interval)
+                new_column_names.append(f"{column_names}_{start}_{end}")
 
         # Return nested pandas DataFrame.
         Xt = pd.DataFrame(_concat_nested_arrays(intervals))
@@ -383,7 +384,8 @@ class SlidingWindowSegmenter(BaseTransformer):
 
     Parameters
     ----------
-        window_length : int, length of interval.
+        window_length : int, optional, default=5.
+            length of sliding window interval
 
     Returns
     -------
diff --git a/sktime/transformations/series/difference.py b/sktime/transformations/series/difference.py
index 49da24d2a..f13ebba39 100644
--- a/sktime/transformations/series/difference.py
+++ b/sktime/transformations/series/difference.py
@@ -2,11 +2,10 @@
 # -*- coding: utf-8 -*-
 # copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
 """Class to iteratively apply differences to a time series."""
-__author__ = ["RNKuhns"]
+__author__ = ["RNKuhns", "fkiraly"]
 __all__ = ["Differencer"]
 
 from typing import Union
-from warnings import warn
 
 import numpy as np
 import pandas as pd
@@ -61,9 +60,6 @@ def _inverse_diff(Z, lag):
     return Z
 
 
-# todo: deprecation in 0.13.0
-#   remove the drop_na *argument* and its handling
-#   change the default behaviour form "drop_na" to "fill_zero"
 class Differencer(BaseTransformer):
     """Apply iterative differences to a timeseries.
 
@@ -87,13 +83,7 @@ class Differencer(BaseTransformer):
         The lags used to difference the data.
         If a single `int` value is
 
-    drop_na : bool, default = True
-        deprecated from 0.12.0, to be removed in 0.13.0
-        Whether the differencer should drop the initial observations that
-        contain missing values as a result of the differencing operation(s).
-
-    na_handling : str, default = "drop_na"
-        default will change to "fill_zero" from 0.13.0
+    na_handling : str, default = "fill_zero"
         How to handle the NaNs that appear at the start of the series from differencing
         Example: there are only 3 differences in a series of length 4,
             differencing [a, b, c, d] gives [?, b-a, c-b, d-c]
@@ -128,27 +118,9 @@ class Differencer(BaseTransformer):
 
     VALID_NA_HANDLING_STR = ["drop_na", "keep_na", "fill_zero"]
 
-    def __init__(self, lags=1, drop_na=None, na_handling="drop_na"):
+    def __init__(self, lags=1, na_handling="fill_zero"):
         self.lags = lags
-        self.drop_na = drop_na
         self.na_handling = self._check_na_handling(na_handling)
-        # note: internally, we will use self._na_handling
-        #   because we must never change input param saves for sklearn compatibility
-        #   and because we need to "translate" the old "drop_na" arg
-        #   this could, in certain cases, overwrite the self. parameter
-        #   and cause sklearn compatibility issues due to the point mentioned
-
-        translate_arg = {True: "drop_na", False: "keep_na"}
-
-        if drop_na is not None:
-            warn(
-                f"the drop_na parameter is deprecated and will be removed in 0.13.0, "
-                f'use na_handling="{translate_arg[drop_na]}" instead',
-                DeprecationWarning,
-            )
-            self._na_handling = translate_arg[drop_na]
-        else:
-            self._na_handling = na_handling
 
         self._Z = None
         self._lags = None
@@ -159,7 +131,7 @@ class Differencer(BaseTransformer):
 
         # if the na_handling is "fill_zero" or "keep_na"
         #   then the returned indices are same to the passed indices
-        if self._na_handling in ["fill_zero", "keep_na"]:
+        if self.na_handling in ["fill_zero", "keep_na"]:
             self.set_tags(**{"transform-returns-same-time-index": True})
 
     def _check_na_handling(self, na_handling):
@@ -255,7 +227,7 @@ class Differencer(BaseTransformer):
         """
         Xt = _diff_transform(X, self._lags)
 
-        na_handling = self._na_handling
+        na_handling = self.na_handling
         if na_handling == "drop_na":
             Xt = Xt.iloc[self._cumulative_lags[-1] :]
         elif na_handling == "fill_zero":
diff --git a/sktime/transformations/series/impute.py b/sktime/transformations/series/impute.py
index fb32f61f4..90f4b4ea1 100644
--- a/sktime/transformations/series/impute.py
+++ b/sktime/transformations/series/impute.py
@@ -7,8 +7,6 @@ __author__ = ["aiwalter"]
 __all__ = ["Imputer"]
 
 
-from warnings import warn
-
 import numpy as np
 from sklearn.utils import check_random_state
 
@@ -170,16 +168,6 @@ class Imputer(BaseTransformer):
         """
         X = X.copy()
 
-        # TODO v0.13.0: Remove this if statement and warning
-        if self.method in ["drift", "mean", "median", "random"]:
-            warn(
-                """Imputer methods "drift", "mean", "median", "random" have been
-                moved to the fit() method, so usage via transform() is
-                deprecated. To still fit on the transform data only, please use the
-                new FitInTransform transformer.
-                """,
-                DeprecationWarning,
-            )
         # replace missing_values with np.nan
         if self.missing_values:
             X = X.replace(to_replace=self.missing_values, value=np.nan)
@@ -190,7 +178,7 @@ class Imputer(BaseTransformer):
         if self.method == "random":
             for col in X.columns:
                 X[col] = X[col].apply(
-                    lambda i: self._get_random(col) if np.isnan(i) else i
+                    lambda i: self._get_random(col) if np.isnan(i) else i  # noqa: B023
                 )
         elif self.method == "constant":
             X = X.fillna(value=self.value)
diff --git a/sktime/transformations/series/summarize.py b/sktime/transformations/series/summarize.py
index 6c25f0ffc..9ac75e86c 100644
--- a/sktime/transformations/series/summarize.py
+++ b/sktime/transformations/series/summarize.py
@@ -6,8 +6,6 @@
 __author__ = ["mloning", "RNKuhns", "danbartl", "grzegorzrut"]
 __all__ = ["SummaryTransformer", "WindowSummarizer"]
 
-import warnings
-
 import pandas as pd
 from joblib import Parallel, delayed
 
@@ -152,19 +150,6 @@ class WindowSummarizer(BaseTransformer):
     >>> transformer = WindowSummarizer(**kwargs)
     >>> y_transformed = transformer.fit_transform(y)
 
-        Example where we transform on a different, later test set:
-    >>> y = load_airline()
-    >>> y_train, y_test = temporal_train_test_split(y)
-    >>> kwargs = {
-    ...     "lag_config": {
-    ...         "lag": ["lag", [[1, 0]]],
-    ...         "mean": ["mean", [[3, 0], [12, 0]]],
-    ...         "std": ["std", [[4, 0]]],
-    ...     }
-    ... }
-    >>> transformer = WindowSummarizer(**kwargs)
-    >>> y_test_transformed = transformer.fit(y_train).transform(y_test)
-
         Example with transforming multiple columns of exogeneous features
     >>> y, X = load_longley()
     >>> y_train, y_test, X_train, X_test = temporal_train_test_split(y, X)
@@ -218,15 +203,12 @@ class WindowSummarizer(BaseTransformer):
 
     def __init__(
         self,
-        lag_config=None,
         lag_feature=None,
         n_jobs=-1,
         target_cols=None,
         truncate=None,
     ):
 
-        # self._converter_store_X = dict()
-        self.lag_config = lag_config
         self.lag_feature = lag_feature
         self.n_jobs = n_jobs
         self.target_cols = target_cols
@@ -275,47 +257,32 @@ class WindowSummarizer(BaseTransformer):
             self._target_cols = self.target_cols
 
         # Convert lag config dictionary to pandas dataframe
-        if self.lag_config is not None:
-            func_dict = pd.DataFrame(self.lag_config).T.reset_index()
-            func_dict.rename(
-                columns={"index": "name", 0: "summarizer", 1: "window"},
-                inplace=True,
-            )
-            func_dict = func_dict.explode("window")
-            func_dict["window"] = func_dict["window"].apply(lambda x: [x[1] + 1, x[0]])
-            func_dict.drop("name", inplace=True, axis=1)
-            warnings.warn(
-                "Specifying lag features via lag_config is deprecated since 0.12.0,"
-                + " and will be removed in 0.13.0. Please use the lag_feature notation"
-                + " (see the documentation for the new notation)."
-            )
+        if self.lag_feature is None:
+            func_dict = pd.DataFrame(
+                {
+                    "lag": [1],
+                }
+            ).T.reset_index()
         else:
-            if self.lag_feature is None:
-                func_dict = pd.DataFrame(
-                    {
-                        "lag": [1],
-                    }
-                ).T.reset_index()
-            else:
-                func_dict = pd.DataFrame.from_dict(
-                    self.lag_feature, orient="index"
-                ).reset_index()
+            func_dict = pd.DataFrame.from_dict(
+                self.lag_feature, orient="index"
+            ).reset_index()
 
-            func_dict = pd.melt(
-                func_dict, id_vars="index", value_name="window", ignore_index=False
-            )
-            func_dict.sort_index(inplace=True)
-            func_dict.drop("variable", axis=1, inplace=True)
-            func_dict.rename(
-                columns={"index": "summarizer"},
-                inplace=True,
-            )
-            func_dict = func_dict.dropna(axis=0, how="any")
-            # Identify lags (since they can follow special notation)
-            lags = func_dict["summarizer"] == "lag"
-            # Convert lags to default list notation with window_length 1
-            boost_lag = func_dict.loc[lags, "window"].apply(lambda x: [int(x), 1])
-            func_dict.loc[lags, "window"] = boost_lag
+        func_dict = pd.melt(
+            func_dict, id_vars="index", value_name="window", ignore_index=False
+        )
+        func_dict.sort_index(inplace=True)
+        func_dict.drop("variable", axis=1, inplace=True)
+        func_dict.rename(
+            columns={"index": "summarizer"},
+            inplace=True,
+        )
+        func_dict = func_dict.dropna(axis=0, how="any")
+        # Identify lags (since they can follow special notation)
+        lags = func_dict["summarizer"] == "lag"
+        # Convert lags to default list notation with window_length 1
+        boost_lag = func_dict.loc[lags, "window"].apply(lambda x: [int(x), 1])
+        func_dict.loc[lags, "window"] = boost_lag
         self.truncate_start = func_dict["window"].apply(lambda x: x[0] + x[1]).max()
         self._func_dict = func_dict
 


Test patch:
diff --git a/sktime/classification/early_classification/tests/test_all_early_classifiers.py b/sktime/classification/early_classification/tests/test_all_early_classifiers.py
index 4b8e160bb..9b4a9029d 100644
--- a/sktime/classification/early_classification/tests/test_all_early_classifiers.py
+++ b/sktime/classification/early_classification/tests/test_all_early_classifiers.py
@@ -5,7 +5,9 @@ __author__ = ["mloning", "TonyBagnall", "fkiraly", "MatthewMiddlehurst"]
 
 import numpy as np
 
-from sktime.classification.tests.test_all_classifiers import TestAllClassifiers
+from sktime.classification.tests.test_all_classifiers import (
+    TestAllClassifiers as ClassifierTests,
+)
 from sktime.tests.test_all_estimators import BaseFixtureGenerator, QuickTester
 
 
@@ -35,7 +37,7 @@ class TestAllEarlyClassifiers(EarlyClassifierFixtureGenerator, QuickTester):
 
     def test_multivariate_input_exception(self, estimator_instance):
         """Test univariate early classifiers raise exception on multivariate X."""
-        test = TestAllClassifiers.test_multivariate_input_exception
+        test = ClassifierTests.test_multivariate_input_exception
         test(self, estimator_instance)
 
     def test_classifier_output(self, estimator_instance, scenario):
@@ -47,7 +49,7 @@ class TestAllEarlyClassifiers(EarlyClassifierFixtureGenerator, QuickTester):
         n_classes = scenario.get_tag("n_classes")
         X_new = scenario.args["predict"]["X"]
         y_train = scenario.args["fit"]["y"]
-        y_pred, decisions, state_info = scenario.run(
+        y_pred, decisions = scenario.run(
             estimator_instance, method_sequence=["fit", "predict"]
         )
 
@@ -63,7 +65,7 @@ class TestAllEarlyClassifiers(EarlyClassifierFixtureGenerator, QuickTester):
         assert isinstance(estimator_instance.get_state_info(), np.ndarray)
 
         # check predict proba (all classifiers have predict_proba by default)
-        y_proba, decisions, state_info = scenario.run(
+        y_proba, decisions = scenario.run(
             estimator_instance, method_sequence=["predict_proba"]
         )
         assert isinstance(y_proba, np.ndarray)
diff --git a/sktime/forecasting/tests/test_all_forecasters.py b/sktime/forecasting/tests/test_all_forecasters.py
index 75bc9c757..44b10bd52 100644
--- a/sktime/forecasting/tests/test_all_forecasters.py
+++ b/sktime/forecasting/tests/test_all_forecasters.py
@@ -30,6 +30,7 @@ from sktime.forecasting.tests._config import (
 from sktime.performance_metrics.forecasting import mean_absolute_percentage_error
 from sktime.tests.test_all_estimators import BaseFixtureGenerator, QuickTester
 from sktime.utils._testing.forecasting import (
+    _assert_correct_columns,
     _assert_correct_pred_time_index,
     _get_expected_index_for_update_predict,
     _get_n_columns,
@@ -224,6 +225,7 @@ class TestAllForecasters(ForecasterFixtureGenerator, QuickTester):
             estimator_instance.fit(y_train, fh=fh)
             y_pred = estimator_instance.predict()
             _assert_correct_pred_time_index(y_pred.index, cutoff, fh=fh_int)
+            _assert_correct_columns(y_pred, y_train)
         except NotImplementedError:
             pass
 
@@ -295,6 +297,7 @@ class TestAllForecasters(ForecasterFixtureGenerator, QuickTester):
             y_pred = estimator_instance.predict(X=X_test)
             cutoff = get_cutoff(y_train, return_index=True)
             _assert_correct_pred_time_index(y_pred.index, cutoff, fh)
+            _assert_correct_columns(y_pred, y_train)
         except NotImplementedError:
             pass
 
@@ -515,6 +518,7 @@ class TestAllForecasters(ForecasterFixtureGenerator, QuickTester):
         )
         cutoff = get_cutoff(y_train, return_index=True)
         _assert_correct_pred_time_index(y_pred.index, cutoff, fh_int_oos)
+        _assert_correct_columns(y_pred, y_train)
 
     @pytest.mark.parametrize(
         "fh_int_oos", TEST_OOS_FHS, ids=[f"fh={fh}" for fh in TEST_OOS_FHS]
diff --git a/sktime/forecasting/tests/test_varmax.py b/sktime/forecasting/tests/test_varmax.py
new file mode 100644
index 000000000..1c39e829f
--- /dev/null
+++ b/sktime/forecasting/tests/test_varmax.py
@@ -0,0 +1,66 @@
+# -*- coding: utf-8 -*-
+"""Tests the VARMAX model."""
+__author__ = ["KatieBuc"]
+
+import numpy as np
+import pandas as pd
+from numpy.testing import assert_allclose
+from statsmodels.tsa.api import VARMAX as _VARMAX
+
+from sktime.forecasting.base import ForecastingHorizon
+from sktime.forecasting.model_selection import temporal_train_test_split
+from sktime.forecasting.varmax import VARMAX
+
+np.random.seed(13455)
+index = pd.date_range(start="2020-01", end="2021-12", freq="M")
+df = pd.DataFrame(
+    np.random.randint(0, 100, size=(23, 3)),
+    columns=list("ABC"),
+    index=pd.PeriodIndex(index),
+)
+
+
+def test_VARMAX_against_statsmodels():
+    """Compares Sktime's and Statsmodel's VARMAX.
+
+    with default variables.
+    """
+    train, _ = temporal_train_test_split(df.astype("float64"))
+    y = train[["A", "B"]]
+
+    sktime_model = VARMAX()
+    fh = ForecastingHorizon([1, 3, 4, 5, 7, 9])
+    sktime_model.fit(y)
+    y_pred = sktime_model.predict(fh=fh)
+
+    stats = _VARMAX(y)
+    stats_fit = stats.fit()
+    start, end = len(train) + fh[0] - 1, len(train) + fh[-1] - 1
+    y_pred_stats = stats_fit.predict(start=start, end=end)
+    y_pred_stats = y_pred_stats.loc[fh.to_absolute(train.index[-1]).to_pandas()]
+
+    assert_allclose(y_pred, y_pred_stats)
+
+
+def test_VARMAX_against_statsmodels_with_exog():
+    """Compares Sktime's and Statsmodel's VARMAX.
+
+    with exogenous input.
+    """
+    train, test = temporal_train_test_split(df.astype("float64"))
+    y_train, X_train = train[["A", "B"]], train[["C"]]
+    _, X_test = test[["A", "B"]], test[["C"]]
+    fh = ForecastingHorizon([1, 2, 3, 4, 5, 6])
+    assert len(fh) == len(X_test)
+
+    sktime_model = VARMAX()
+    sktime_model.fit(y_train, X=X_train)
+    y_pred = sktime_model.predict(fh=fh, X=X_test)
+
+    stats = _VARMAX(y_train, exog=X_train)
+    stats_fit = stats.fit()
+    start, end = len(train) + fh[0] - 1, len(train) + fh[-1] - 1
+    y_pred_stats = stats_fit.predict(start=start, end=end, exog=X_test)
+    y_pred_stats = y_pred_stats.loc[fh.to_absolute(train.index[-1]).to_pandas()]
+
+    assert_allclose(y_pred, y_pred_stats)
diff --git a/sktime/tests/_config.py b/sktime/tests/_config.py
index f0a143ead..c7a408c69 100644
--- a/sktime/tests/_config.py
+++ b/sktime/tests/_config.py
@@ -81,6 +81,15 @@ EXCLUDED_TESTS = {
     "SeriesToSeriesRowTransformer": ["test_methods_do_not_change_state"],
     # ColumnTransformer still needs to be refactored, see #2537
     "ColumnTransformer": ["test_methods_do_not_change_state"],
+    # Early classifiers intentionally retain information from pervious predict calls
+    #   for #1.
+    # #2 amd #3 are due to predict/predict_proba returning two items and that breaking
+    #   assert_array_equal
+    "TEASER": [
+        "test_methods_do_not_change_state",
+        "test_fit_idempotent",
+        "test_persistence_via_pickle",
+    ],
 }
 
 # We here configure estimators for basic unit testing, including setting of
diff --git a/sktime/transformations/panel/tests/test_segment.py b/sktime/transformations/panel/tests/test_segment.py
index 5f2e1311c..f7f105952 100644
--- a/sktime/transformations/panel/tests/test_segment.py
+++ b/sktime/transformations/panel/tests/test_segment.py
@@ -1,21 +1,24 @@
 # -*- coding: utf-8 -*-
+"""Tests for RandomIntervalSegmenter."""
 import numpy as np
 import pandas as pd
 import pytest
 
-from sktime.transformations.panel.segment import RandomIntervalSegmenter
-from sktime.transformations.panel.segment import _rand_intervals_fixed_n
-from sktime.transformations.panel.segment import _rand_intervals_rand_n
+from sktime.transformations.panel.segment import (
+    RandomIntervalSegmenter,
+    _rand_intervals_fixed_n,
+    _rand_intervals_rand_n,
+)
 from sktime.utils._testing.panel import _make_nested_from_array
 
 N_ITER = 10
 
 
-# Test output format and dimensions.
 @pytest.mark.parametrize("n_instances", [1, 3])
 @pytest.mark.parametrize("n_timepoints", [10, 20])
 @pytest.mark.parametrize("n_intervals", [0.1, 1.0, 1, 3, 10, "sqrt", "random", "log"])
 def test_output_format_dim(n_timepoints, n_instances, n_intervals):
+    """Test output format and dimensions."""
     X = _make_nested_from_array(
         np.ones(n_timepoints), n_instances=n_instances, n_columns=1
     )
@@ -30,29 +33,29 @@ def test_output_format_dim(n_timepoints, n_instances, n_intervals):
     # Check number of generated intervals/columns.
     if n_intervals != "random":
         if np.issubdtype(type(n_intervals), np.floating):
-            assert Xt.shape[1] == np.maximum(1, int(n_timepoints * n_intervals))
+            assert Xt.shape[1] <= np.maximum(1, int(n_timepoints * n_intervals))
         elif np.issubdtype(type(n_intervals), np.integer):
-            assert Xt.shape[1] == n_intervals
+            assert Xt.shape[1] <= n_intervals
         elif n_intervals == "sqrt":
-            assert Xt.shape[1] == np.maximum(1, int(np.sqrt(n_timepoints)))
+            assert Xt.shape[1] <= np.maximum(1, int(np.sqrt(n_timepoints)))
         elif n_intervals == "log":
-            assert Xt.shape[1] == np.maximum(1, int(np.log(n_timepoints)))
+            assert Xt.shape[1] <= np.maximum(1, int(np.log(n_timepoints)))
 
 
-# Check that exception is raised for bad input args.
 @pytest.mark.parametrize("bad_interval", [0, -0, "str", 1.2, -1.2, -1])
 def test_bad_input_args(bad_interval):
+    """Check that exception is raised for bad input args."""
     X = _make_nested_from_array(np.ones(10), n_instances=10, n_columns=2)
     with pytest.raises(ValueError):
         RandomIntervalSegmenter(n_intervals=bad_interval).fit(X)
 
 
-# Helper function for checking generated intervals.
 @pytest.mark.parametrize(
     "random_state", list(np.random.randint(100, size=10))
 )  # run repeatedly
 @pytest.mark.parametrize("n_intervals", ["sqrt", "log", 0.1, 1, 3])
 def test_rand_intervals_fixed_n(random_state, n_intervals):
+    """Helper function for checking generated intervals."""
     n_timepoints = 30
     x = np.arange(n_timepoints)
 
@@ -73,6 +76,7 @@ def test_rand_intervals_fixed_n(random_state, n_intervals):
     "random_state", list(np.random.randint(100, size=10))
 )  # run repeatedly
 def test_rand_intervals_rand_n(random_state):
+    """Test random intervals."""
     n_timepoints = 30
     x = np.arange(n_timepoints)
 
@@ -93,6 +97,7 @@ def test_rand_intervals_rand_n(random_state):
 @pytest.mark.parametrize("max_length", [4, 5])
 @pytest.mark.parametrize("n_intervals", ["sqrt", "log", 0.1, 1, 3])
 def test_rand_intervals_fixed_n_min_max_length(n_intervals, min_length, max_length):
+    """Test random interval length."""
     n_timepoints = 30
     x = np.arange(n_timepoints)
 
diff --git a/sktime/transformations/panel/tests/test_sliding_window_segmenter_transformer.py b/sktime/transformations/panel/tests/test_sliding_window_segmenter_transformer.py
index 9866468ce..dfed8e46d 100644
--- a/sktime/transformations/panel/tests/test_sliding_window_segmenter_transformer.py
+++ b/sktime/transformations/panel/tests/test_sliding_window_segmenter_transformer.py
@@ -1,7 +1,13 @@
 # -*- coding: utf-8 -*-
+"""Tests for SlidingWindowSegmenter."""
+
+__authors__ = ["mloning"]
+
+
 import numpy as np
 import pandas as pd
 import pytest
+
 from sktime.transformations.panel.segment import SlidingWindowSegmenter
 from sktime.utils._testing.panel import _make_nested_from_array
 
@@ -11,6 +17,7 @@ from sktime.utils._testing.panel import _make_nested_from_array
 # correct input is meant to be a positive integer of 1 or more.
 @pytest.mark.parametrize("bad_window_length", ["str", 1.2, -1.2, -1, {}])
 def test_bad_input_args(bad_window_length):
+    """Test that bad inputs raise value error."""
     X = _make_nested_from_array(np.ones(10), n_instances=10, n_columns=1)
 
     if not isinstance(bad_window_length, int):
@@ -23,6 +30,7 @@ def test_bad_input_args(bad_window_length):
 
 # Check the transformer has changed the data correctly.
 def test_output_of_transformer():
+    """Test correct output of SlidingWindowSegmenter."""
     X = _make_nested_from_array(
         np.array([1, 2, 3, 4, 5, 6]), n_instances=1, n_columns=1
     )
@@ -66,6 +74,7 @@ def test_output_of_transformer():
     "time_series_length,window_length", [(5, 1), (10, 5), (15, 9), (20, 13), (25, 19)]
 )
 def test_output_dimensions(time_series_length, window_length):
+    """Test output dimensions of SlidingWindowSegmenter."""
     X = _make_nested_from_array(
         np.ones(time_series_length), n_instances=10, n_columns=1
     )
@@ -83,17 +92,8 @@ def test_output_dimensions(time_series_length, window_length):
     assert num_cols == time_series_length
 
 
-# Test that subsequence transformer fails when a multivariate ts
-# is fed into it.
-def test_fails_if_multivariate():
-    X = _make_nested_from_array(np.ones(5), n_instances=10, n_columns=5)
-
-    with pytest.raises(ValueError):
-        SlidingWindowSegmenter().fit(X).transform(X)
-
-
 def convert_list_to_dataframe(list_to_convert):
-    # Convert this into a panda's data frame
+    """Convert list into a pandas data frame."""
     df = pd.DataFrame()
     for i in range(len(list_to_convert)):
         inst = list_to_convert[i]
@@ -105,9 +105,7 @@ def convert_list_to_dataframe(list_to_convert):
 
 
 def check_if_dataframes_are_equal(df1, df2):
-    """
-    for some reason, this is how you check that two dataframes are equal.
-    """
+    """Check that pandas DataFrames are equal."""
     from pandas.testing import assert_frame_equal
 
     try:
diff --git a/sktime/transformations/tests/test_base.py b/sktime/transformations/tests/test_base.py
index 1c1a73e2c..a848fdcdc 100644
--- a/sktime/transformations/tests/test_base.py
+++ b/sktime/transformations/tests/test_base.py
@@ -16,9 +16,10 @@ __all__ = []
 import sys
 from inspect import isclass
 
+import pandas as pd
 import pytest
 
-from sktime.datatypes import check_is_scitype, mtype_to_scitype
+from sktime.datatypes import check_is_scitype, get_examples, mtype_to_scitype
 from sktime.transformations.compose import FitInTransform
 from sktime.transformations.panel.padder import PaddingTransformer
 from sktime.transformations.panel.tsfresh import (
@@ -36,6 +37,7 @@ from sktime.utils._testing.scenarios_transformers import (
     TransformerFitTransformSeriesMultivariate,
     TransformerFitTransformSeriesUnivariate,
 )
+from sktime.utils._testing.series import _make_series
 
 # other scenarios that might be needed later in development:
 # TransformerFitTransformPanelUnivariateWithClassY,
@@ -590,3 +592,35 @@ def test_vectorization_multivariate_and_hierarchical_empty_fit():
     # length of Xt should be number of hierarchy levels times number of time points
     assert len(Xt) == len(scenario.args["fit"]["X"])
     assert len(Xt.columns) == len(scenario.args["fit"]["X"].columns)
+
+
+def test_vectorize_reconstruct_unique_columns():
+    """Tests that vectorization on multivariate output yields unique columns.
+
+    Also test that the column names are as expected:
+    <variable>__<transformed> if multiple transformed variables per variable are present
+    <variable> if one variable is transformed to one output
+
+    Raises
+    ------
+    AssertionError if output columns are not as expected.
+    """
+    from sktime.transformations.series.detrend import Detrender
+    from sktime.transformations.series.theta import ThetaLinesTransformer
+
+    X = pd.DataFrame({"a": [1, 2], "b": [3, 4], "c": [5, 6]})
+    X_mi = get_examples("pd_multiindex_hier")[0]
+
+    t = ThetaLinesTransformer()
+
+    X_t_cols = t.fit_transform(X).columns
+
+    assert set(X_t_cols) == set(["a__0", "a__2", "b__0", "b__2", "c__0", "c__2"])
+
+    X_mi_cols = t.fit_transform(X_mi)
+    assert set(X_mi_cols) == set(["var_0__0", "var_0__2", "var_1__0", "var_1__2"])
+
+    X = _make_series(n_columns=2, n_timepoints=15)
+    t = Detrender.create_test_instance()
+    Xt = t.fit_transform(X)
+    assert set(Xt.columns) == set([0, 1])
diff --git a/sktime/utils/_testing/estimator_checks.py b/sktime/utils/_testing/estimator_checks.py
index c758f7666..6e44c5b03 100644
--- a/sktime/utils/_testing/estimator_checks.py
+++ b/sktime/utils/_testing/estimator_checks.py
@@ -135,7 +135,7 @@ def _make_fit_args(estimator, **kwargs):
     elif isinstance(estimator, BaseSeriesAnnotator):
         X = make_annotation_problem(**kwargs)
         return (X,)
-    elif isinstance(estimator, BaseClassifier):
+    elif isinstance(estimator, (BaseClassifier, BaseEarlyClassifier)):
         return make_classification_problem(**kwargs)
     elif isinstance(estimator, BaseRegressor):
         return make_regression_problem(**kwargs)
@@ -353,7 +353,7 @@ def _has_capability(est, method: str) -> bool:
         "predict_quantiles",
         "predict_var",
     ]:
-        ALWAYS_HAVE_PREDICT_PROBA = (BaseClassifier, BaseClusterer)
+        ALWAYS_HAVE_PREDICT_PROBA = (BaseClassifier, BaseEarlyClassifier, BaseClusterer)
         # all classifiers and clusterers implement predict_proba
         if method == "predict_proba" and isinstance(est, ALWAYS_HAVE_PREDICT_PROBA):
             return True
diff --git a/sktime/utils/_testing/forecasting.py b/sktime/utils/_testing/forecasting.py
index 0d15af855..d4645426a 100644
--- a/sktime/utils/_testing/forecasting.py
+++ b/sktime/utils/_testing/forecasting.py
@@ -138,6 +138,16 @@ def _assert_correct_pred_time_index(y_pred_index, cutoff, fh):
     y_pred_index.equals(expected)
 
 
+def _assert_correct_columns(y_pred, y_train):
+    """Check that forecast object has right column names."""
+    if isinstance(y_pred, pd.DataFrame) and isinstance(y_train, pd.DataFrame):
+        msg = (
+            "forecast must have same columns index as past data, "
+            f"expected {y_train.columns} but found {y_pred.columns}"
+        )
+        assert (y_pred.columns == y_train.columns).all(), msg
+
+
 def _make_fh(cutoff, steps, fh_type, is_relative):
     """Construct forecasting horizons for testing."""
     from sktime.forecasting.tests._config import INDEX_TYPE_LOOKUP
diff --git a/sktime/utils/_testing/scenarios_classification.py b/sktime/utils/_testing/scenarios_classification.py
index 6be95f445..bc2f70db1 100644
--- a/sktime/utils/_testing/scenarios_classification.py
+++ b/sktime/utils/_testing/scenarios_classification.py
@@ -6,13 +6,18 @@ Contains TestScenario concrete children to run in tests for classifiers/regressi
 
 __author__ = ["fkiraly"]
 
-__all__ = ["scenarios_classification", "scenarios_regression"]
+__all__ = [
+    "scenarios_classification",
+    "scenarios_early_classification",
+    "scenarios_regression",
+]
 
 from copy import deepcopy
 from inspect import isclass
 
 from sktime.base import BaseObject
 from sktime.classification.base import BaseClassifier
+from sktime.classification.early_classification import BaseEarlyClassifier
 from sktime.regression.base import BaseRegressor
 from sktime.utils._testing.hierarchical import _make_hierarchical
 from sktime.utils._testing.panel import _make_classification_y, _make_panel_X
@@ -73,10 +78,11 @@ class ClassifierTestScenario(TestScenario, BaseObject):
             else:
                 return obj.get_tag(tag_name)
 
-        regr_or_classf = (BaseClassifier, BaseRegressor)
+        regr_or_classf = (BaseClassifier, BaseEarlyClassifier, BaseRegressor)
 
-        # applicable only if obj inherits from BaseClassifier or BaseRegressor
-        #   currently we test both classifiers and regressors using these scenarios
+        # applicable only if obj inherits from BaseClassifier, BaseEarlyClassifier or
+        #   BaseRegressor. currently we test both classifiers and regressors using these
+        #   scenarios
         if not isinstance(obj, regr_or_classf) and not issubclass(obj, regr_or_classf):
             return False
 
@@ -173,6 +179,13 @@ scenarios_classification = [
     ClassifierFitPredictUnequalLength,
 ]
 
+# same scenarios used for early classification
+scenarios_early_classification = [
+    ClassifierFitPredict,
+    ClassifierFitPredictMultivariate,
+    ClassifierFitPredictUnequalLength,
+]
+
 # we use the same scenarios for regression, as in the old test suite
 scenarios_regression = [
     ClassifierFitPredict,
diff --git a/sktime/utils/_testing/scenarios_getter.py b/sktime/utils/_testing/scenarios_getter.py
index a61e8234c..728274e18 100644
--- a/sktime/utils/_testing/scenarios_getter.py
+++ b/sktime/utils/_testing/scenarios_getter.py
@@ -13,6 +13,7 @@ from sktime.registry import BASE_CLASS_LIST, BASE_CLASS_SCITYPE_LIST, scitype
 from sktime.utils._testing.scenarios_aligners import scenarios_aligners
 from sktime.utils._testing.scenarios_classification import (
     scenarios_classification,
+    scenarios_early_classification,
     scenarios_regression,
 )
 from sktime.utils._testing.scenarios_clustering import scenarios_clustering
@@ -26,6 +27,7 @@ from sktime.utils._testing.scenarios_transformers_pairwise import (
 scenarios = dict()
 scenarios["aligner"] = scenarios_aligners
 scenarios["classifier"] = scenarios_classification
+scenarios["early_classifier"] = scenarios_early_classification
 scenarios["clusterer"] = scenarios_clustering
 scenarios["forecaster"] = scenarios_forecasting
 scenarios["regressor"] = scenarios_regression