Repository: sktime/sktime
Pull request #5613 (2023-12-09 14:34:23)

[MNT] 0.25.0 deprecations and change actions
Carries out deprecations and change actions scheduled for 0.25.0:

* remove deprecated `forecasting.model_selection._split` module
* in forecasting `evaluate`, change order of return data frame
* in forecasting `evaluate`, remove `compute` arg. `kwargs` are still accepted, with warning message. Bumped removal to 0.26.0 for smoother transition.
* bump checks


Patch:
diff --git a/README.md b/README.md
index e9e621732..b7d2cb4b6 100644
--- a/README.md
+++ b/README.md
@@ -169,8 +169,8 @@ flexible choice of soft dependencies is unavailable via `conda`)
 ``` python
 from sktime.datasets import load_airline
 from sktime.forecasting.base import ForecastingHorizon
-from sktime.forecasting.model_selection import temporal_train_test_split
 from sktime.forecasting.theta import ThetaForecaster
+from sktime.split import temporal_train_test_split
 from sktime.performance_metrics.forecasting import mean_absolute_percentage_error
 
 y = load_airline()
diff --git a/docs/source/api_reference/transformations.rst b/docs/source/api_reference/transformations.rst
index 25513ffd2..0bbc48474 100644
--- a/docs/source/api_reference/transformations.rst
+++ b/docs/source/api_reference/transformations.rst
@@ -665,6 +665,19 @@ Bootstrap transformations
     SplitterBootstrapTransformer
     STLBootstrapTransformer
 
+Panel-to-Series transformers
+----------------------------
+
+These transformers create a single series from a panel.
+
+.. currentmodule:: sktime.transformations.merger
+
+.. autosummary::
+    :toctree: auto_generated/
+    :template: class.rst
+
+    Merger
+
 Outlier detection, changepoint detection
 ----------------------------------------
 
diff --git a/docs/source/get_started.rst b/docs/source/get_started.rst
index cce05ec13..5295cf594 100644
--- a/docs/source/get_started.rst
+++ b/docs/source/get_started.rst
@@ -75,9 +75,9 @@ Forecasting
 
     >>> from sktime.datasets import load_airline
     >>> from sktime.forecasting.base import ForecastingHorizon
-    >>> from sktime.forecasting.model_selection import temporal_train_test_split
     >>> from sktime.forecasting.theta import ThetaForecaster
     >>> from sktime.performance_metrics.forecasting import mean_absolute_percentage_error
+    >>> from sktime.split import temporal_train_test_split
 
     >>> y = load_airline()
     >>> y_train, y_test = temporal_train_test_split(y)
diff --git a/docs/source/user_guide/introduction.rst b/docs/source/user_guide/introduction.rst
index 672dca8a6..70f6c7fd8 100644
--- a/docs/source/user_guide/introduction.rst
+++ b/docs/source/user_guide/introduction.rst
@@ -31,8 +31,8 @@ can simply write:
     from sktime.datasets import load_airline
     from sktime.forecasting.compose import make_reduction
     from sklearn.ensemble import RandomForestRegressor
-    from sktime.forecasting.model_selection import temporal_train_test_split
     from sktime.performance_metrics.forecasting import MeanAbsolutePercentageError
+    from sktime.split import temporal_train_test_split
 
     y = load_airline()
     y_train, y_test = temporal_train_test_split(y)
diff --git a/extension_templates/classification.py b/extension_templates/classification.py
index 532743cd5..884d7f77f 100644
--- a/extension_templates/classification.py
+++ b/extension_templates/classification.py
@@ -71,9 +71,11 @@ class MyTimeSeriesClassifier(BaseClassifier):
     # these are the default values, only add if different to these.
     _tags = {
         "X_inner_mtype": "numpy3D",  # which type do _fit/_predict accept, usually
-        # this is either "numpy3D" or "nested_univ" (nested pd.DataFrame). Other
+        "y_inner_mtype": "numpy1D",  # which type do _fit/_predict return, usually
+        # this is either "numpy3D", "pd-multiindex" or "nested_univ" (nested df). Other
         # types are allowable, see datatypes/panel/_registry.py for options.
-        "capability:multivariate": False,
+        "capability:multivariate": False,  # ability to handle multivariate X
+        "capability:multioutput": False,  # ability to predict multiple columns in y
         "capability:unequal_length": False,
         "capability:missing_values": False,
         "capability:train_estimate": False,
@@ -129,13 +131,19 @@ class MyTimeSeriesClassifier(BaseClassifier):
         ----------
         X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
             if self.get_tag("X_inner_mtype") = "numpy3D":
-                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
-            if self.get_tag("X_inner_mtype") = "nested_univ":
-                pd.DataFrame with each column a dimension, each cell a pd.Series
+            3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "pd-multiindex:":
+            pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
             for list of other mtypes, see datatypes.SCITYPE_REGISTER
             for specifications, see examples/AA_datatypes_and_datasets.ipynb
-        y : 1D np.array of int, of shape [n_instances] - class labels for fitting
-            indices correspond to instance indices in X
+        y : guaranteed to be of a type in self.get_tag("y_inner_mtype")
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            class labels for fitting
+            if self.get_tag("capaility:multioutput") = False, guaranteed to be 1D
+            if self.get_tag("capaility:multioutput") = True, guaranteed to be 2D
 
         Returns
         -------
@@ -152,31 +160,32 @@ class MyTimeSeriesClassifier(BaseClassifier):
         #   3. read from self in _fit,  4. pass to interfaced_model.fit in _fit
 
     # todo: implement this, mandatory
-    def _predict(self, X) -> np.ndarray:
+    def _predict(self, X):
         """Predict labels for sequences in X.
 
         private _predict containing the core logic, called from predict
 
-        State required:
-            Requires state to be "fitted".
-
-        Accesses in self:
-            Fitted model attributes ending in "_"
-
         Parameters
         ----------
         X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
             if self.get_tag("X_inner_mtype") = "numpy3D":
-                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
-            if self.get_tag("X_inner_mtype") = "nested_univ":
-                pd.DataFrame with each column a dimension, each cell a pd.Series
+            3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "pd-multiindex:":
+            pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
             for list of other mtypes, see datatypes.SCITYPE_REGISTER
             for specifications, see examples/AA_datatypes_and_datasets.ipynb
 
         Returns
         -------
-        y : 1D np.array of int, of shape [n_instances] - predicted class labels
+        y : should be of mtype in self.get_tag("y_inner_mtype")
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            predicted class labels
             indices correspond to instance indices in X
+            if self.get_tag("capaility:multioutput") = False, should be 1D
+            if self.get_tag("capaility:multioutput") = True, should be 2D
         """
 
         # implement here
@@ -200,9 +209,11 @@ class MyTimeSeriesClassifier(BaseClassifier):
         ----------
         X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
             if self.get_tag("X_inner_mtype") = "numpy3D":
-                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
-            if self.get_tag("X_inner_mtype") = "nested_univ":
-                pd.DataFrame with each column a dimension, each cell a pd.Series
+            3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "pd-multiindex:":
+            pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
             for list of other mtypes, see datatypes.SCITYPE_REGISTER
             for specifications, see examples/AA_datatypes_and_datasets.ipynb
 
diff --git a/sktime/classification/base.py b/sktime/classification/base.py
index 7d8110f49..74728f65d 100644
--- a/sktime/classification/base.py
+++ b/sktime/classification/base.py
@@ -31,9 +31,11 @@ import pandas as pd
 from sktime.base import BaseEstimator
 from sktime.datatypes import (
     MTYPE_LIST_PANEL,
+    MTYPE_LIST_TABLE,
+    VectorizedDF,
     check_is_error_msg,
     check_is_scitype,
-    convert_to,
+    convert,
 )
 from sktime.utils.sklearn import is_sklearn_transformer
 from sktime.utils.validation import check_n_jobs
@@ -60,7 +62,9 @@ class BaseClassifier(BaseEstimator, ABC):
     _tags = {
         "object_type": "classifier",  # type of object
         "X_inner_mtype": "numpy3D",  # which type do _fit/_predict, support for X?
+        "y_inner_mtype": "numpy1D",  # which type do _fit/_predict, support for y?
         #    it should be either "numpy3D" or "nested_univ" (nested pd.DataFrame)
+        "capability:multioutput": False,  # whether classifier supports multioutput
         "capability:multivariate": False,
         "capability:unequal_length": False,
         "capability:missing_values": False,
@@ -93,6 +97,8 @@ class BaseClassifier(BaseEstimator, ABC):
         # required for compatibility with some sklearn interfaces
         # i.e. CalibratedClassifierCV
         self._estimator_type = "classifier"
+        self._is_vectorized = False
+        self._converter_store_y = {}
 
         super().__init__()
         _check_estimator_deps(self)
@@ -136,35 +142,60 @@ class BaseClassifier(BaseEstimator, ABC):
             return NotImplemented
 
     def fit(self, X, y):
-        """Fit time series classifier to training data.
+        """
+        Fit time series classifier to training data.
+
+        State change:
+            Changes state to "fitted".
+
+        Writes to self:
+            Sets self.is_fitted to True.
+            Sets fitted model attributes ending in "_".
 
         Parameters
         ----------
-        X : 3D np.array (any number of dimensions, equal length series)
-                of shape [n_instances, n_dimensions, series_length]
-            or 2D np.array (univariate, equal length series)
-                of shape [n_instances, series_length]
-            or pd.DataFrame with each column a dimension, each cell a pd.Series
-                (any number of dimensions, equal or unequal length series)
+        X : sktime compatible time series panel data container, Panel scitype, e.g.,
+            pd-multiindex: pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
+            numpy3D: 3D np.array (any number of dimensions, equal length series)
+            of shape [n_instances, n_dimensions, series_length]
             or of any other supported Panel mtype
-                for list of mtypes, see datatypes.SCITYPE_REGISTER
-                for specifications, see examples/AA_datatypes_and_datasets.ipynb
-        y : 1D np.array of int, of shape [n_instances] - class labels for fitting
-            indices correspond to instance indices in X
+            for list of mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
+        y : sktime compatible tabular data container, Table scitype
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            class labels for fitting
+            0-th indices correspond to instance indices in X
+            1-st indices (if applicable) correspond to multioutput vector indices in X
+            supported sktime types: np.ndarray (1D, 2D), pd.Series, pd.DataFrame
 
         Returns
         -------
         self : Reference to self.
-
-        Notes
-        -----
-        Changes state by creating a fitted model that updates attributes
-        ending in "_" and sets is_fitted flag to True.
         """
         # reset estimator at the start of fit
         self.reset()
 
+        # fit timer start
         start = int(round(time.time() * 1000))
+
+        # check and convert y for multioutput vectorization
+        y, y_metadata = self._check_y(y)
+        self._y_metadata = y_metadata
+        self._is_vectorized = isinstance(y, VectorizedDF)
+
+        if self._is_vectorized:
+            self._vectorize("fit", X=X, y=y)
+            # fit timer end
+            self.fit_time_ = int(round(time.time() * 1000)) - start
+            # this should happen last: fitted state is set to True
+            self._is_fitted = True
+            return self
+
+        # no vectorization needed, proceed with normal fit
+
         # convenience conversions to allow user flexibility:
         # if X is 2D array, convert to 3D, if y is Series, convert to numpy
         X, y = self._internal_convert(X, y)
@@ -174,6 +205,7 @@ class BaseClassifier(BaseEstimator, ABC):
         missing = X_metadata["has_nans"]
         multivariate = not X_metadata["is_univariate"]
         unequal = not X_metadata["is_equal_length"]
+        X_mtype = X_metadata["mtype"]
         self._X_metadata = X_metadata
 
         # Check this classifier can handle characteristics
@@ -194,7 +226,7 @@ class BaseClassifier(BaseEstimator, ABC):
             return self
 
         # Convert data as dictated by the classifier tags
-        X = self._convert_X(X)
+        X = self._convert_X(X, X_mtype)
         multithread = self.get_tag("capability:multithreading")
         if multithread:
             try:
@@ -208,32 +240,42 @@ class BaseClassifier(BaseEstimator, ABC):
         self._fit(X, y)
         self.fit_time_ = int(round(time.time() * 1000)) - start
 
-        # this should happen last
+        # this should happen last: fitted state is set to True
         self._is_fitted = True
         return self
 
-    def predict(self, X) -> np.ndarray:
+    def predict(self, X):
         """Predicts labels for sequences in X.
 
         Parameters
         ----------
-        X : 3D np.array (any number of dimensions, equal length series)
-                of shape [n_instances, n_dimensions, series_length]
-            or 2D np.array (univariate, equal length series)
-                of shape [n_instances, series_length]
-            or pd.DataFrame with each column a dimension, each cell a pd.Series
-                (any number of dimensions, equal or unequal length series)
+        X : sktime compatible time series panel data container, Panel scitype, e.g.,
+            pd-multiindex: pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
+            numpy3D: 3D np.array (any number of dimensions, equal length series)
+            of shape [n_instances, n_dimensions, series_length]
             or of any other supported Panel mtype
-                for list of mtypes, see datatypes.SCITYPE_REGISTER
-                for specifications, see examples/AA_datatypes_and_datasets.ipynb
+            for list of mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
 
         Returns
         -------
-        y : 1D np.array of int, of shape [n_instances] - predicted class labels
-            indices correspond to instance indices in X
+        y_pred : sktime compatible tabular data container, Table scitype
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            predicted class labels
+            0-th indices correspond to instance indices in X
+            1-st indices (if applicable) correspond to multioutput vector indices in X
+            1D np.npdarray, if y univariate (one dimension)
+            otherwise, same type as y passed in fit
         """
         self.check_is_fitted()
 
+        # enter vectorized mode if needed
+        if self._is_vectorized:
+            return self._vectorize("predict", X=X)
+
         # boilerplate input checks for predict-like methods
         X = self._check_convert_X_for_predict(X)
 
@@ -241,33 +283,42 @@ class BaseClassifier(BaseEstimator, ABC):
         if len(self._class_dictionary) == 1:
             return self._single_class_y_pred(X, method="predict")
 
-        # call internal _predict_proba
-        return self._predict(X)
+        # call internal _predict, convert output
+        y_pred_inner = self._predict(X)
+        y_pred = self._convert_output_y(y_pred_inner)
+        return y_pred
 
-    def predict_proba(self, X) -> np.ndarray:
+    def predict_proba(self, X):
         """Predicts labels probabilities for sequences in X.
 
         Parameters
         ----------
-        X : 3D np.array (any number of dimensions, equal length series)
-                of shape [n_instances, n_dimensions, series_length]
-            or 2D np.array (univariate, equal length series)
-                of shape [n_instances, series_length]
-            or pd.DataFrame with each column a dimension, each cell a pd.Series
-                (any number of dimensions, equal or unequal length series)
+        X : sktime compatible time series panel data container, Panel scitype, e.g.,
+            pd-multiindex: pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
+            numpy3D: 3D np.array (any number of dimensions, equal length series)
+            of shape [n_instances, n_dimensions, series_length]
             or of any other supported Panel mtype
-                for list of mtypes, see datatypes.SCITYPE_REGISTER
-                for specifications, see examples/AA_datatypes_and_datasets.ipynb
+            for list of mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
 
         Returns
         -------
-        y : 2D array of shape [n_instances, n_classes] - predicted class probabilities
-            1st dimension indices correspond to instance indices in X
-            2nd dimension indices correspond to possible labels (integers)
-            (i, j)-th entry is predictive probability that i-th instance is of class j
+        y_pred : 2D np.array of int, of shape [n_instances, n_classes]
+            predicted class label probabilities
+            0-th indices correspond to instance indices in X
+            1-st indices correspond to class index, in same order as in self.classes_
+            entries are predictive class probabilities, summing to 1
         """
         self.check_is_fitted()
 
+        # enter vectorized mode if needed
+        if self._is_vectorized:
+            return self._vectorize("predict_proba", X=X)
+
+        self.check_is_fitted()
+
         # boilerplate input checks for predict-like methods
         X = self._check_convert_X_for_predict(X)
 
@@ -278,7 +329,7 @@ class BaseClassifier(BaseEstimator, ABC):
         # call internal _predict_proba
         return self._predict_proba(X)
 
-    def fit_predict(self, X, y, cv=None, change_state=True) -> np.ndarray:
+    def fit_predict(self, X, y, cv=None, change_state=True):
         """Fit and predict labels for sequences in X.
 
         Convenience method to produce in-sample predictions and
@@ -292,17 +343,22 @@ class BaseClassifier(BaseEstimator, ABC):
 
         Parameters
         ----------
-        X : 3D np.array (any number of dimensions, equal length series)
-                of shape [n_instances, n_dimensions, series_length]
-            or 2D np.array (univariate, equal length series)
-                of shape [n_instances, series_length]
-            or pd.DataFrame with each column a dimension, each cell a pd.Series
-                (any number of dimensions, equal or unequal length series)
+        X : sktime compatible time series panel data container, Panel scitype, e.g.,
+            pd-multiindex: pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
+            numpy3D: 3D np.array (any number of dimensions, equal length series)
+            of shape [n_instances, n_dimensions, series_length]
             or of any other supported Panel mtype
-                for list of mtypes, see datatypes.SCITYPE_REGISTER
-                for specifications, see examples/AA_datatypes_and_datasets.ipynb
-        y : 1D np.array of int, of shape [n_instances] - class labels for fitting
-            indices correspond to instance indices in X
+            for list of mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
+        y : sktime compatible tabular data container, Table scitype
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            class labels for fitting
+            0-th indices correspond to instance indices in X
+            1-st indices (if applicable) correspond to multioutput vector indices in X
+            supported sktime types: np.ndarray (1D, 2D), pd.Series, pd.DataFrame
         cv : None, int, or sklearn cross-validation object, optional, default=None
             None : predictions are in-sample, equivalent to fit(X, y).predict(X)
             cv : predictions are equivalent to fit(X_train, y_train).predict(X_test)
@@ -320,14 +376,59 @@ class BaseClassifier(BaseEstimator, ABC):
 
         Returns
         -------
-        y : 1D np.array of int, of shape [n_instances] - predicted class labels
-            indices correspond to instance indices in X
-            if cv is passed, -1 indicates entries not seen in union of test sets
+        y_pred : sktime compatible tabular data container, Table scitype
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            predicted class labels
+            0-th indices correspond to instance indices in X
+            1-st indices (if applicable) correspond to multioutput vector indices in X
+            1D np.npdarray, if y univariate (one dimension)
+            otherwise, same type as y passed in fit
         """
         return self._fit_predict_boilerplate(
             X=X, y=y, cv=cv, change_state=change_state, method="predict"
         )
 
+    def _vectorize(self, methodname, **kwargs):
+        """Vectorized/iterated loop over method of BaseClassifier.
+
+        Uses classifiers_ attribute to store one classifier per loop index.
+        """
+        y = kwargs.get("y")
+        X = kwargs.get("X")
+        if X is not None:
+            kwargs.pop("X")
+        if y is not None:
+            kwargs.pop("y")
+            self._y_vec = y
+        classifiers_ = self._y_vec.vectorize_est(
+            self,
+            method="clone",
+        )
+        if methodname == "fit":
+            self.classifiers_ = self._y_vec.vectorize_est(
+                classifiers_,
+                method=methodname,
+                args={"y": y},
+                X=X,
+            )
+            return self
+        else:
+            if self.classifiers_ is not None:
+                classifiers_ = self.classifiers_
+            y_pred = self._y_vec.vectorize_est(
+                classifiers_,
+                method=methodname,
+                # return_type="list",
+                X=X,
+                args={"y": y} if y is not None else {},
+                **kwargs,  # contains X inside
+            )
+            y_pred = pd.DataFrame(
+                {str(i): y_pred[col].values[0] for i, col in enumerate(y_pred.columns)}
+            )
+            return y_pred
+
     def _fit_predict_boilerplate(self, X, y, cv, change_state, method):
         """Boilerplate logic for fit_predict and fit_predict_proba."""
         from sklearn.model_selection import KFold
@@ -355,6 +456,7 @@ class BaseClassifier(BaseEstimator, ABC):
         missing = X_metadata["has_nans"]
         multivariate = not X_metadata["is_univariate"]
         unequal = not X_metadata["is_equal_length"]
+        X_mtype = X_metadata["mtype"]
         # Check this classifier can handle characteristics
         self._check_capabilities(missing, multivariate, unequal)
 
@@ -364,15 +466,17 @@ class BaseClassifier(BaseEstimator, ABC):
 
         # Convert data to format easily usable for applying cv
         if isinstance(X, np.ndarray):
-            X = convert_to(
+            X = convert(
                 X,
+                from_type=X_mtype,
                 to_type="numpy3D",
                 as_scitype="Panel",
                 store_behaviour="freeze",
             )
         else:
-            X = convert_to(
+            X = convert(
                 X,
+                from_type=X_mtype,
                 to_type="nested_univ",
                 as_scitype="Panel",
                 store_behaviour="freeze",
@@ -400,25 +504,36 @@ class BaseClassifier(BaseEstimator, ABC):
 
         return y_pred
 
-    def fit_predict_proba(self, X, y, cv=None, change_state=True) -> np.ndarray:
+    def fit_predict_proba(self, X, y, cv=None, change_state=True):
         """Fit and predict labels probabilities for sequences in X.
 
         Convenience method to produce in-sample predictions and
         cross-validated out-of-sample predictions.
 
+        Writes to self, if change_state=True:
+            Sets self.is_fitted to True.
+            Sets fitted model attributes ending in "_".
+
+        Does not update state if change_state=False.
+
         Parameters
         ----------
-        X : 3D np.array (any number of dimensions, equal length series)
-                of shape [n_instances, n_dimensions, series_length]
-            or 2D np.array (univariate, equal length series)
-                of shape [n_instances, series_length]
-            or pd.DataFrame with each column a dimension, each cell a pd.Series
-                (any number of dimensions, equal or unequal length series)
+        X : sktime compatible time series panel data container, Panel scitype, e.g.,
+            pd-multiindex: pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
+            numpy3D: 3D np.array (any number of dimensions, equal length series)
+            of shape [n_instances, n_dimensions, series_length]
             or of any other supported Panel mtype
-                for list of mtypes, see datatypes.SCITYPE_REGISTER
-                for specifications, see examples/AA_datatypes_and_datasets.ipynb
-        y : 1D np.array of int, of shape [n_instances] - class labels for fitting
-            indices correspond to instance indices in X
+            for list of mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
+        y : sktime compatible tabular data container, Table scitype
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            class labels for fitting
+            0-th indices correspond to instance indices in X
+            1-st indices (if applicable) correspond to multioutput vector indices in X
+            supported sktime types: np.ndarray (1D, 2D), pd.Series, pd.DataFrame
         cv : None, int, or sklearn cross-validation object, optional, default=None
             None : predictions are in-sample, equivalent to fit(X, y).predict(X)
             cv : predictions are equivalent to fit(X_train, y_train).predict(X_test)
@@ -434,10 +549,11 @@ class BaseClassifier(BaseEstimator, ABC):
 
         Returns
         -------
-        y : 2D array of shape [n_instances, n_classes] - predicted class probabilities
-            1st dimension indices correspond to instance indices in X
-            2nd dimension indices correspond to possible labels (integers)
-            (i, j)-th entry is predictive probability that i-th instance is of class j
+        y_pred : 2D np.array of int, of shape [n_instances, n_classes]
+            predicted class label probabilities
+            0-th indices correspond to instance indices in X
+            1-st indices correspond to class index, in same order as in self.classes_
+            entries are predictive class probabilities, summing to 1
         """
         return self._fit_predict_boilerplate(
             X=X, y=y, cv=cv, change_state=change_state, method="predict_proba"
@@ -460,17 +576,22 @@ class BaseClassifier(BaseEstimator, ABC):
 
         Parameters
         ----------
-        X : 3D np.array (any number of dimensions, equal length series)
-                of shape [n_instances, n_dimensions, series_length]
-            or 2D np.array (univariate, equal length series)
-                of shape [n_instances, series_length]
-            or pd.DataFrame with each column a dimension, each cell a pd.Series
-                (any number of dimensions, equal or unequal length series)
+        X : sktime compatible time series panel data container, e.g.,
+            pd-multiindex: pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
+            numpy3D: 3D np.array (any number of dimensions, equal length series)
+            of shape [n_instances, n_dimensions, series_length]
             or of any other supported Panel mtype
-                for list of mtypes, see datatypes.SCITYPE_REGISTER
-                for specifications, see examples/AA_datatypes_and_datasets.ipynb
-        y : 1D np.ndarray of int, of shape [n_instances] - class labels (ground truth)
-            indices correspond to instance indices in X
+            for list of mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
+        y : sktime compatible tabular data container, Table scitype
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            class labels for fitting
+            0-th indices correspond to instance indices in X
+            1-st indices (if applicable) correspond to multioutput vector indices in X
+            supported sktime types: np.ndarray (1D, 2D), pd.Series, pd.DataFrame
 
         Returns
         -------
@@ -512,27 +633,30 @@ class BaseClassifier(BaseEstimator, ABC):
 
         Abstract method, must be implemented.
 
+        Writes to self:
+            Sets fitted model attributes ending in "_".
+
         Parameters
         ----------
         X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
             if self.get_tag("X_inner_mtype") = "numpy3D":
-                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
-            if self.get_tag("X_inner_mtype") = "nested_univ":
-                pd.DataFrame with each column a dimension, each cell a pd.Series
+            3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "pd-multiindex:":
+            pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
             for list of other mtypes, see datatypes.SCITYPE_REGISTER
             for specifications, see examples/AA_datatypes_and_datasets.ipynb
-        y : 1D np.array of int, of shape [n_instances] - class labels for fitting
-            indices correspond to instance indices in X
+        y : guaranteed to be of a type in self.get_tag("y_inner_mtype")
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            class labels for fitting
+            if self.get_tag("capaility:multioutput") = False, guaranteed to be 1D
+            if self.get_tag("capaility:multioutput") = True, guaranteed to be 2D
 
         Returns
         -------
-        self :
-            Reference to self.
-
-        Notes
-        -----
-        Changes state by creating a fitted model that updates attributes
-        ending in "_" and sets is_fitted flag to True.
+        self : Reference to self.
         """
         ...
 
@@ -545,16 +669,23 @@ class BaseClassifier(BaseEstimator, ABC):
         ----------
         X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
             if self.get_tag("X_inner_mtype") = "numpy3D":
-                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
-            if self.get_tag("X_inner_mtype") = "nested_univ":
-                pd.DataFrame with each column a dimension, each cell a pd.Series
+            3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "pd-multiindex:":
+            pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
             for list of other mtypes, see datatypes.SCITYPE_REGISTER
             for specifications, see examples/AA_datatypes_and_datasets.ipynb
 
         Returns
         -------
-        y : 1D np.array of int, of shape [n_instances] - predicted class labels
+        y : should be of mtype in self.get_tag("y_inner_mtype")
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            predicted class labels
             indices correspond to instance indices in X
+            if self.get_tag("capaility:multioutput") = False, should be 1D
+            if self.get_tag("capaility:multioutput") = True, should be 2D
         """
         y_proba = self._predict_proba(X)
         y_pred = y_proba.argmax(axis=1)
@@ -572,9 +703,11 @@ class BaseClassifier(BaseEstimator, ABC):
         ----------
         X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
             if self.get_tag("X_inner_mtype") = "numpy3D":
-                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
-            if self.get_tag("X_inner_mtype") = "nested_univ":
-                pd.DataFrame with each column a dimension, each cell a pd.Series
+            3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "pd-multiindex:":
+            pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
             for list of other mtypes, see datatypes.SCITYPE_REGISTER
             for specifications, see examples/AA_datatypes_and_datasets.ipynb
 
@@ -617,10 +750,11 @@ class BaseClassifier(BaseEstimator, ABC):
         missing = X_metadata["has_nans"]
         multivariate = not X_metadata["is_univariate"]
         unequal = not X_metadata["is_equal_length"]
+        X_mtype = X_metadata["mtype"]
         # Check this classifier can handle characteristics
         self._check_capabilities(missing, multivariate, unequal)
         # Convert data as dictated by the classifier tags
-        X = self._convert_X(X)
+        X = self._convert_X(X, X_mtype=X_mtype)
 
         return X
 
@@ -670,29 +804,121 @@ class BaseClassifier(BaseEstimator, ABC):
             else:
                 raise ValueError(msg)
 
-    def _convert_X(self, X):
+    def _convert_X(self, X, X_mtype):
         """Convert equal length series from DataFrame to numpy array or vice versa.
 
         Parameters
         ----------
-        self : this classifier
-        X : pd.DataFrame or np.ndarray. Input attribute data
+        X : input data for the classifier
+        X_mtype : str, a Panel mtype string, e.g., "pd_multiindex", "numpy3D"
 
         Returns
         -------
         X : input X converted to type in "X_inner_mtype" tag
-                usually a pd.DataFrame (nested) or 3D np.ndarray
+            usually a pd.DataFrame (nested) or 3D np.ndarray
             Checked and possibly converted input data
         """
         inner_type = self.get_tag("X_inner_mtype")
         # convert pd.DataFrame
-        X = convert_to(
+        X = convert(
             X,
+            from_type=X_mtype,
             to_type=inner_type,
             as_scitype="Panel",
         )
         return X
 
+    def _check_y(self, y=None):
+        """Check and coerce X/y for fit/transform functions.
+
+        Parameters
+        ----------
+        y : pd.DataFrame, pd.Series or np.ndarray
+
+        Returns
+        -------
+        y : object of sktime compatible time series type
+            can be Series, Panel, Hierarchical
+        y_metadata : dict
+            metadata of y, retured by check_is_scitype
+        """
+        if y is None:
+            return None
+
+        capa_multioutput = self.get_tag("capability:multioutput")
+        y_inner_mtype = self.get_tag("y_inner_mtype")
+
+        y_valid, y_msg, y_metadata = check_is_scitype(
+            y, "Table", return_metadata=["is_univariate"]
+        )
+
+        if not y_valid:
+            allowed_msg = (
+                f"In classification, y must be of a supported type, "
+                f"for instance 1D or 2D numpy arrays, pd.DataFrame, or pd.Series. "
+                f"Allowed compatible mtype format specifications are:"
+                f" {MTYPE_LIST_TABLE} ."
+            )
+            check_is_error_msg(
+                y_msg, var_name="y", allowed_msg=allowed_msg, raise_exception=True
+            )
+
+        y_uni = y_metadata["is_univariate"]
+        y_mtype = y_metadata["mtype"]
+
+        requires_vectorization = not capa_multioutput and not y_uni
+
+        if requires_vectorization:
+            y_df = convert(
+                y,
+                from_type=y_mtype,
+                to_type="pd_DataFrame_Table",
+                as_scitype="Table",
+                store=self._converter_store_y,
+            )
+            y_vec = VectorizedDF([y_df], iterate_cols=True)
+            return y_vec, y_metadata
+
+        y_inner = convert(
+            y,
+            from_type=y_mtype,
+            to_type=y_inner_mtype,
+            as_scitype="Table",
+            store=self._converter_store_y,
+        )
+
+        return y_inner, y_metadata
+
+    def _convert_output_y(self, y):
+        """Convert output y to original format.
+
+        Parameters
+        ----------
+        y : np.ndarray or pd.DataFrame
+
+        Returns
+        -------
+        y : np.ndarray or pd.DataFrame
+        """
+        # for consistency with legacy behaviour:
+        # output is coerced to numpy1D in case of univariate output
+        if not self._y_metadata["is_univariate"]:
+            output_mtype = self._y_metadata["mtype"]
+            converter_store = self._converter_store_y
+        else:
+            output_mtype = "numpy1D"
+            converter_store = None
+
+        y = convert(
+            y,
+            from_type=self.get_tag("y_inner_mtype"),
+            to_type=output_mtype,
+            as_scitype="Table",
+            store=converter_store,
+            store_behaviour="freeze",
+        )
+        return y
+
     def _check_classifier_input(
         self, X, y=None, enforce_min_instances=1, return_metadata=True
     ):
@@ -744,9 +970,10 @@ class BaseClassifier(BaseEstimator, ABC):
         # Check y if passed
         if y is not None:
             # Check y valid input
-            if not isinstance(y, (pd.Series, np.ndarray)):
+            if not isinstance(y, (pd.Series, pd.DataFrame, np.ndarray)):
                 raise ValueError(
-                    f"y must be a np.array or a pd.Series, but found type: {type(y)}"
+                    "y must be a np.array or a pd.Series or pd.DataFrame, but found ",
+                    f"type: {type(y)}",
                 )
             # Check matching number of labels
             n_labels = y.shape[0]
@@ -756,9 +983,9 @@ class BaseClassifier(BaseEstimator, ABC):
                     f"{n_labels}"
                 )
             if isinstance(y, np.ndarray):
-                if y.ndim > 1:
+                if y.ndim > 2:
                     raise ValueError(
-                        f"np.ndarray y must be 1-dimensional, "
+                        f"np.ndarray y must be 1-dimensional or 2-dimensional, "
                         f"but found {y.ndim} dimensions"
                     )
             # warn if only a single class label is seen
@@ -794,9 +1021,6 @@ class BaseClassifier(BaseEstimator, ABC):
             # done here, but touches a lot of files, so will get this to work first.
             if X.ndim == 2:
                 X = X.reshape(X.shape[0], 1, X.shape[1])
-        if y is not None and isinstance(y, pd.Series):
-            # y should be a numpy array, although we allow Series for user convenience
-            y = pd.Series.to_numpy(y)
         if y is None:
             return X
         return X, y
diff --git a/sktime/classification/dictionary_based/_boss.py b/sktime/classification/dictionary_based/_boss.py
index cfcfb6ca9..88889a05a 100644
--- a/sktime/classification/dictionary_based/_boss.py
+++ b/sktime/classification/dictionary_based/_boss.py
@@ -666,6 +666,7 @@ class IndividualBOSS(BaseClassifier):
         new_boss.n_classes_ = self.n_classes_
         new_boss.classes_ = self.classes_
         new_boss._class_dictionary = self._class_dictionary
+        new_boss._y_metadata = self._y_metadata
         new_boss._is_fitted = True
 
         return new_boss
diff --git a/sktime/classification/early_classification/_probability_threshold.py b/sktime/classification/early_classification/_probability_threshold.py
index 0464bfd3e..4697bf9de 100644
--- a/sktime/classification/early_classification/_probability_threshold.py
+++ b/sktime/classification/early_classification/_probability_threshold.py
@@ -20,7 +20,7 @@ from sktime.classification.interval_based import CanonicalIntervalForest
 from sktime.utils.validation.panel import check_X
 
 
-# TODO: fix this in 0.25.0
+# TODO: fix this in 0.26.0
 # base class should have been changed to BaseEarlyClassifier
 class ProbabilityThresholdEarlyClassifier(BaseClassifier):
     """Probability Threshold Early Classifier.
diff --git a/sktime/classification/early_classification/base.py b/sktime/classification/early_classification/base.py
index 335cb38b7..3e6642569 100644
--- a/sktime/classification/early_classification/base.py
+++ b/sktime/classification/early_classification/base.py
@@ -55,7 +55,9 @@ class BaseEarlyClassifier(BaseEstimator, ABC):
     _tags = {
         "object_type": "early_classifier",  # type of object
         "X_inner_mtype": "numpy3D",  # which type do _fit/_predict, support for X?
+        "y_inner_mtype": "numpy1D",  # which type do _fit/_predict, support for y?
         #    it should be either "numpy3D" or "nested_univ" (nested pd.DataFrame)
+        "capability:multioutput": False,  # whether classifier supports multioutput
         "capability:multivariate": False,
         "capability:unequal_length": False,
         "capability:missing_values": False,
@@ -89,6 +91,8 @@ class BaseEarlyClassifier(BaseEstimator, ABC):
         """
         self.state_info = None
 
+        self._converter_store_y = {}
+
         super().__init__()
 
     def fit(self, X, y):
@@ -603,22 +607,37 @@ class BaseEarlyClassifier(BaseEstimator, ABC):
         _check_capabilities = BaseClassifier._check_capabilities
         return _check_capabilities(self, missing, multivariate, unequal)
 
-    def _convert_X(self, X):
+    def _convert_X(self, X, X_mtype):
         """Convert equal length series from DataFrame to numpy array or vice versa.
 
         Parameters
         ----------
-        self : this classifier
-        X : pd.DataFrame or np.ndarray. Input attribute data
+        X : input data for the classifier
+        X_mtype : str, a Panel mtype string, e.g., "pd_multiindex", "numpy3D"
 
         Returns
         -------
         X : input X converted to type in "X_inner_mtype" tag
-                usually a pd.DataFrame (nested) or 3D np.ndarray
+            usually a pd.DataFrame (nested) or 3D np.ndarray
             Checked and possibly converted input data
         """
         _convert_X = BaseClassifier._convert_X
-        return _convert_X(self, X)
+        return _convert_X(self, X, X_mtype)
+
+    def _check_y(self, y=None):
+        """Check and coerce X/y for fit/transform functions.
+
+        Parameters
+        ----------
+        y : pd.DataFrame, pd.Series or np.ndarray
+
+        Returns
+        -------
+        y : object of sktime compatible time series type
+            can be Series, Panel, Hierarchical
+        """
+        _check_y = BaseClassifier._check_y
+        return _check_y(self, y)
 
     def _check_classifier_input(
         self, X, y=None, enforce_min_instances=1, return_metadata=True
diff --git a/sktime/datatypes/_table/_check.py b/sktime/datatypes/_table/_check.py
index 6b25cf333..4a1d96d88 100644
--- a/sktime/datatypes/_table/_check.py
+++ b/sktime/datatypes/_table/_check.py
@@ -65,11 +65,6 @@ def check_pddataframe_table(obj, return_metadata=False, var_name="obj"):
     if _req("has_nans", return_metadata):
         metadata["has_nans"] = obj.isna().values.any()
 
-    # check that no dtype is object
-    if "object" in obj.dtypes.values:
-        msg = f"{var_name} should not have column of 'object' dtype"
-        return _ret(False, msg, None, return_metadata)
-
     return _ret(True, None, metadata, return_metadata)
 
 
@@ -92,11 +87,6 @@ def check_pdseries_table(obj, return_metadata=False, var_name="obj"):
     if _req("n_instances", return_metadata):
         metadata["n_instances"] = len(index)
 
-    # check that dtype is not object
-    if "object" == obj.dtypes:
-        msg = f"{var_name} should not be of 'object' dtype"
-        return _ret(False, msg, None, return_metadata)
-
     # check whether index is equally spaced or if there are any nans
     #   compute only if needed
     if _req("has_nans", return_metadata):
diff --git a/sktime/forecasting/all/__init__.py b/sktime/forecasting/all/__init__.py
index 2d72d78b8..7ac4e1a8e 100644
--- a/sktime/forecasting/all/__init__.py
+++ b/sktime/forecasting/all/__init__.py
@@ -11,13 +11,6 @@ import pandas as pd
 from sktime.datasets import load_airline, load_longley, load_lynx, load_shampoo_sales
 from sktime.forecasting.base import ForecastingHorizon
 from sktime.forecasting.model_evaluation import evaluate
-from sktime.forecasting.model_selection import (
-    CutoffSplitter,
-    ExpandingWindowSplitter,
-    SingleWindowSplitter,
-    SlidingWindowSplitter,
-    temporal_train_test_split,
-)
 from sktime.performance_metrics.forecasting import (
     GeometricMeanRelativeAbsoluteError,
     GeometricMeanRelativeSquaredError,
@@ -58,6 +51,13 @@ from sktime.performance_metrics.forecasting import (
     relative_loss,
 )
 from sktime.registry import all_estimators
+from sktime.split import (
+    CutoffSplitter,
+    ExpandingWindowSplitter,
+    SingleWindowSplitter,
+    SlidingWindowSplitter,
+    temporal_train_test_split,
+)
 from sktime.transformations.series.detrend import Deseasonalizer, Detrender
 from sktime.utils.plotting import plot_series
 
diff --git a/sktime/forecasting/base/_fh.py b/sktime/forecasting/base/_fh.py
index ff5341c96..7a71c6e9a 100644
--- a/sktime/forecasting/base/_fh.py
+++ b/sktime/forecasting/base/_fh.py
@@ -806,7 +806,7 @@ def _to_relative(fh: ForecastingHorizon, cutoff=None) -> ForecastingHorizon:
             absolute = _coerce_to_period(absolute, freq=fh.freq)
             cutoff = _coerce_to_period(cutoff, freq=fh.freq)
 
-        # TODO: 0.25.0:
+        # TODO: 0.26.0:
         # Check at every minor release whether lower pandas bound >=0.15.0
         # if yes, can remove the workaround in the "else" condition and the check
         #
diff --git a/sktime/forecasting/compose/_pipeline.py b/sktime/forecasting/compose/_pipeline.py
index b802dfa38..cae9a4c1e 100644
--- a/sktime/forecasting/compose/_pipeline.py
+++ b/sktime/forecasting/compose/_pipeline.py
@@ -170,7 +170,7 @@ class _Pipeline(_HeterogenousMetaEstimator, BaseForecaster):
                         if len(levels) == 1:
                             levels = levels[0]
                         yt[ix] = y.xs(ix, level=levels, axis=1)
-                        # todo 0.25.0 - check why this cannot be easily removed
+                        # todo 0.26.0 - check why this cannot be easily removed
                         # in theory, we should get rid of the "Coverage" case treatment
                         # (the legacy naming convention was removed in 0.23.0)
                         # deal with the "Coverage" case, we need to get rid of this
diff --git a/sktime/forecasting/model_evaluation/_functions.py b/sktime/forecasting/model_evaluation/_functions.py
index f2082a86d..a0a62c01c 100644
--- a/sktime/forecasting/model_evaluation/_functions.py
+++ b/sktime/forecasting/model_evaluation/_functions.py
@@ -318,8 +318,6 @@ def _evaluate_window(x, meta):
         return result
 
 
-# todo 0.25.0: remove compute argument and docstring
-# todo 0.25.0: remove kwargs and docstring
 def evaluate(
     forecaster,
     cv,
@@ -330,7 +328,6 @@ def evaluate(
     return_data: bool = False,
     error_score: Union[str, int, float] = np.nan,
     backend: Optional[str] = None,
-    compute: bool = None,
     cv_X=None,
     backend_params: Optional[dict] = None,
     **kwargs,
@@ -422,9 +419,6 @@ def evaluate(
         "threading" is unlikely to see speed ups due to the GIL and the serialization
         backend (``cloudpickle``) for "dask" and "loky" is generally more robust
         than the standard ``pickle`` library used in "multiprocessing".
-    compute : bool, default=True, deprecated and will be removed in 0.25.0.
-        If backend="dask", whether returned DataFrame is computed.
-        If set to True, returns `pd.DataFrame`, otherwise `dask.dataframe.DataFrame`.
     cv_X : sktime BaseSplitter descendant, optional
         determines split of ``X`` into test and train folds
         default is ``X`` being split to identical ``loc`` indices as ``y``
@@ -525,10 +519,10 @@ def evaluate(
                 "installed, but dask is not present in the python environment"
             )
 
-    # todo 0.25.0: remove kwargs and this warning
+    # todo 0.26.0: remove kwargs and this warning
     if kwargs != {}:
         warnings.warn(
-            "in evaluate, kwargs will no longer be supported from sktime 0.25.0. "
+            "in evaluate, kwargs are no longer supported. "
             "to pass configuration arguments to the parallelization backend, "
             "use backend_params instead. "
             f"The following kwargs were found: {kwargs.keys()}, pass these as "
@@ -537,28 +531,8 @@ def evaluate(
             stacklevel=2,
         )
 
-    # todo 0.25.0: remove compute argument and logic, and remove this warning
-    if compute is not None:
-        warnings.warn(
-            "the compute argument of evaluate is deprecated and will be removed "
-            "in sktime 0.25.0. For the same behaviour in the future, "
-            'use backend="dask_lazy"',
-            DeprecationWarning,
-            stacklevel=2,
-        )
-    if compute is None:
-        compute = True
-    if backend == "dask" and not compute:
-        backend = "dask_lazy"
-
     _check_strategy(strategy)
     cv = check_cv(cv, enforce_start_with_window=True)
-    # TODO: remove lines(four lines below) and 599-612 in v0.25.0
-    if isinstance(scoring, list):
-        raise_warn, num = True, len(scoring)
-    else:
-        raise_warn, num = False, 1
-    # removal until here
     scoring = _check_scores(scoring)
 
     ALLOWED_SCITYPES = ["Series", "Panel", "Hierarchical"]
@@ -669,21 +643,4 @@ def evaluate(
     # final formatting of results DataFrame
     results = results.reset_index(drop=True)
 
-    # TODO: remove 16 lines below and 451-455 in v0.25.0
-    if raise_warn:
-        warnings.warn(
-            "Starting v0.25.0 model_evaluation.evaluate module will rearrange "
-            "all metric columns to the left of its output result DataFrame. "
-            "Please use loc references when addressing the columns. You can "
-            "safely ignore this warning if you don't use evaluate function directly.",
-            DeprecationWarning,
-            stacklevel=2,
-        )
-        columns = results.columns.to_list()
-        non_first_metrics = []
-        for _ in range(1, num):
-            metric = columns.pop(1)
-            non_first_metrics.append(metric)
-        results = results.reindex(columns=columns + non_first_metrics)
-    #  removal until here
     return results
diff --git a/sktime/forecasting/model_selection/__init__.py b/sktime/forecasting/model_selection/__init__.py
index c55f73960..311a2157e 100644
--- a/sktime/forecasting/model_selection/__init__.py
+++ b/sktime/forecasting/model_selection/__init__.py
@@ -4,29 +4,11 @@
 
 __author__ = ["mloning", "kkoralturk"]
 __all__ = [
-    "CutoffSplitter",
-    "SameLocSplitter",
-    "SingleWindowSplitter",
-    "SlidingWindowSplitter",
-    "temporal_train_test_split",
-    "ExpandingGreedySplitter",
-    "ExpandingWindowSplitter",
-    "TestPlusTrainSplitter",
     "ForecastingGridSearchCV",
     "ForecastingRandomizedSearchCV",
     "ForecastingSkoptSearchCV",
 ]
 
-from sktime.forecasting.model_selection._split import (
-    CutoffSplitter,
-    ExpandingGreedySplitter,
-    ExpandingWindowSplitter,
-    SameLocSplitter,
-    SingleWindowSplitter,
-    SlidingWindowSplitter,
-    TestPlusTrainSplitter,
-    temporal_train_test_split,
-)
 from sktime.forecasting.model_selection._tune import (
     ForecastingGridSearchCV,
     ForecastingRandomizedSearchCV,
diff --git a/sktime/forecasting/model_selection/_split.py b/sktime/forecasting/model_selection/_split.py
deleted file mode 100644
index 3d77d9d77..000000000
--- a/sktime/forecasting/model_selection/_split.py
+++ /dev/null
@@ -1,42 +0,0 @@
-#!/usr/bin/env python3 -u
-# copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
-"""Implement dataset splitting for model evaluation and selection."""
-
-# This module is deprecated. Below the deprecation cycle :
-# - Until the 0.25.0 release, imports are allowed from this location but raise a
-# DeprecationWarning.
-# - In the 0.25.0 release, this file (sktime/forecasting/model_selection/_split.py)
-# will be removed.
-# todo 0.25.0 : Please ensure this file is deleted.
-
-__all__ = [
-    "ExpandingGreedySplitter",
-    "ExpandingWindowSplitter",
-    "SlidingWindowSplitter",
-    "CutoffSplitter",
-    "SingleWindowSplitter",
-    "SameLocSplitter",
-    "temporal_train_test_split",
-    "TestPlusTrainSplitter",
-]
-
-import warnings
-
-from sktime.split import (
-    CutoffSplitter,
-    ExpandingGreedySplitter,
-    ExpandingWindowSplitter,
-    SameLocSplitter,
-    SingleWindowSplitter,
-    SlidingWindowSplitter,
-    TestPlusTrainSplitter,
-    temporal_train_test_split,
-)
-
-warnings.warn(
-    "Please import sktime splitters from `sktime.split`. Importing splitters from "
-    "`sktime.forecasting.model_selection` is deprecated and will be removed in release "
-    "0.25.0",
-    DeprecationWarning,
-    stacklevel=2,
-)
diff --git a/sktime/forecasting/model_selection/_tune.py b/sktime/forecasting/model_selection/_tune.py
index fd0fcd585..5c11825bf 100644
--- a/sktime/forecasting/model_selection/_tune.py
+++ b/sktime/forecasting/model_selection/_tune.py
@@ -544,9 +544,7 @@ class ForecastingGridSearchCV(BaseGridSearch):
     --------
     >>> from sktime.datasets import load_shampoo_sales
     >>> from sktime.forecasting.model_selection import ForecastingGridSearchCV
-    >>> from sktime.split import (
-    ...     ExpandingWindowSplitter,
-    ...     ExpandingWindowSplitter)
+    >>> from sktime.split import ExpandingWindowSplitter
     >>> from sktime.forecasting.naive import NaiveForecaster
     >>> y = load_shampoo_sales()
     >>> fh = [1,2,3]
diff --git a/sktime/forecasting/trend/_pwl_trend_forecaster.py b/sktime/forecasting/trend/_pwl_trend_forecaster.py
index c3d947660..e59bc3fb9 100644
--- a/sktime/forecasting/trend/_pwl_trend_forecaster.py
+++ b/sktime/forecasting/trend/_pwl_trend_forecaster.py
@@ -60,8 +60,8 @@ class ProphetPiecewiseLinearTrendForecaster(_ProphetAdapter):
     --------
     >>> from sktime.datasets import load_airline
     >>> from sktime.forecasting.trend import ProphetPiecewiseLinearTrendForecaster
-    >>> from sktime.forecasting.model_selection import temporal_train_test_split
     >>> from sktime.forecasting.base import ForecastingHorizon
+    >>> from sktime.split import temporal_train_test_split
     >>> y =load_airline().to_timestamp(freq='M')
     >>> y_train, y_test = temporal_train_test_split(y)
     >>> fh = ForecastingHorizon(y.index, is_relative=False)
diff --git a/sktime/registry/_tags.py b/sktime/registry/_tags.py
index 1290a351d..b40bb47b3 100644
--- a/sktime/registry/_tags.py
+++ b/sktime/registry/_tags.py
@@ -492,6 +492,12 @@ ESTIMATOR_TAG_REGISTER = [
         "int",
         "sample size used in approximating other statistics if not available",
     ),
+    (
+        "capability:multioutput",
+        ["classifier", "regressor"],  # might need to add "early_classifier" here
+        "bool",
+        "can the estimator handle multioutput data?",
+    ),
 ]
 
 ESTIMATOR_TAG_TABLE = pd.DataFrame(ESTIMATOR_TAG_REGISTER)
diff --git a/sktime/regression/base.py b/sktime/regression/base.py
index 903cc0cf7..4c0351608 100644
--- a/sktime/regression/base.py
+++ b/sktime/regression/base.py
@@ -31,9 +31,11 @@ import pandas as pd
 from sktime.base import BaseEstimator
 from sktime.datatypes import (
     MTYPE_LIST_PANEL,
+    MTYPE_LIST_TABLE,
+    VectorizedDF,
     check_is_error_msg,
     check_is_scitype,
-    convert_to,
+    convert,
 )
 from sktime.utils.sklearn import is_sklearn_transformer
 from sktime.utils.validation import check_n_jobs
@@ -57,7 +59,9 @@ class BaseRegressor(BaseEstimator, ABC):
     _tags = {
         "object_type": "regressor",  # type of object
         "X_inner_mtype": "numpy3D",  # which type do _fit/_predict, support for X?
+        "y_inner_mtype": "numpy1D",  # which type do _fit/_predict, support for y?
         #    it should be either "numpy3D" or "nested_univ" (nested pd.DataFrame)
+        "capability:multioutput": False,  # whether regressor supports multioutput
         "capability:multivariate": False,
         "capability:unequal_length": False,
         "capability:missing_values": False,
@@ -84,6 +88,9 @@ class BaseRegressor(BaseEstimator, ABC):
         # required for compatibility with some sklearn interfaces
         # i.e. CalibratedRegressorCV
         self._estimator_type = "regressor"
+        self._is_vectorized = False
+        self._is_timed = False
+        self._converter_store_y = {}
 
         super().__init__()
 
@@ -128,19 +135,31 @@ class BaseRegressor(BaseEstimator, ABC):
     def fit(self, X, y):
         """Fit time series regressor to training data.
 
+        State change:
+            Changes state to "fitted".
+
+        Writes to self:
+            Sets self.is_fitted to True.
+            Sets fitted model attributes ending in "_".
+
         Parameters
         ----------
-        X : 3D np.array (any number of dimensions, equal length series)
-                of shape [n_instances, n_dimensions, series_length]
-            or 2D np.array (univariate, equal length series)
-                of shape [n_instances, series_length]
-            or pd.DataFrame with each column a dimension, each cell a pd.Series
-                (any number of dimensions, equal or unequal length series)
+        X : sktime compatible time series panel data container, Panel scitype, e.g.,
+            pd-multiindex: pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
+            numpy3D: 3D np.array (any number of dimensions, equal length series)
+            of shape [n_instances, n_dimensions, series_length]
             or of any other supported Panel mtype
-                for list of mtypes, see datatypes.SCITYPE_REGISTER
-                for specifications, see examples/AA_datatypes_and_datasets.ipynb
-        y : 1D np.array of float, of shape [n_instances] - regression labels for fitting
-            indices correspond to instance indices in X
+            for list of mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
+        y : sktime compatible tabular data container, Table scitype
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            class labels for fitting
+            0-th indices correspond to instance indices in X
+            1-st indices (if applicable) correspond to multioutput vector indices in X
+            supported sktime types: np.ndarray (1D, 2D), pd.Series, pd.DataFrame
 
         Returns
         -------
@@ -151,7 +170,26 @@ class BaseRegressor(BaseEstimator, ABC):
         Changes state by creating a fitted model that updates attributes
         ending in "_" and sets is_fitted flag to True.
         """
+        self.reset()
+
+        # fit timer start
         start = int(round(time.time() * 1000))
+
+        # check and convert y for multioutput vectorization
+        y, y_metadata = self._check_y(y)
+        self._y_metadata = y_metadata
+        self._is_vectorized = isinstance(y, VectorizedDF)
+
+        if self._is_vectorized:
+            self._vectorize("fit", X=X, y=y)
+            # fit timer end
+            self.fit_time_ = int(round(time.time() * 1000)) - start
+            # this should happen last: fitted state is set to True
+            self._is_fitted = True
+            return self
+
+        # no vectorization needed, proceed with normal fit
+
         # convenience conversions to allow user flexibility:
         # if X is 2D array, convert to 3D, if y is Series, convert to numpy
         X, y = _internal_convert(X, y)
@@ -162,10 +200,14 @@ class BaseRegressor(BaseEstimator, ABC):
         missing = X_metadata["has_nans"]
         multivariate = not X_metadata["is_univariate"]
         unequal = not X_metadata["is_equal_length"]
+        X_mtype = X_metadata["mtype"]
+        self._X_metadata = X_metadata
+
         # Check this regressor can handle characteristics
         self._check_capabilities(missing, multivariate, unequal)
+
         # Convert data as dictated by the regressor tags
-        X = self._convert_X(X)
+        X = self._convert_X(X, X_mtype)
         multithread = self.get_tag("capability:multithreading")
         if multithread:
             try:
@@ -177,7 +219,8 @@ class BaseRegressor(BaseEstimator, ABC):
 
         self._fit(X, y)
         self.fit_time_ = int(round(time.time() * 1000)) - start
-        # this should happen last
+
+        # this should happen last: fitted state is set to True
         self._is_fitted = True
         return self
 
@@ -186,44 +229,64 @@ class BaseRegressor(BaseEstimator, ABC):
 
         Parameters
         ----------
-        X : 3D np.array (any number of dimensions, equal length series)
-                of shape [n_instances, n_dimensions, series_length]
-            or 2D np.array (univariate, equal length series)
-                of shape [n_instances, series_length]
-            or pd.DataFrame with each column a dimension, each cell a pd.Series
-                (any number of dimensions, equal or unequal length series)
+        X : sktime compatible time series panel data container, Panel scitype, e.g.,
+            pd-multiindex: pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
+            numpy3D: 3D np.array (any number of dimensions, equal length series)
+            of shape [n_instances, n_dimensions, series_length]
             or of any other supported Panel mtype
-                for list of mtypes, see datatypes.SCITYPE_REGISTER
-                for specifications, see examples/AA_datatypes_and_datasets.ipynb
+            for list of mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
 
         Returns
         -------
-        y : 1D np.array of float, of shape [n_instances] - predicted regression labels
-            indices correspond to instance indices in X
+        y_pred : sktime compatible tabular data container, Table scitype
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            predicted class labels
+            0-th indices correspond to instance indices in X
+            1-st indices (if applicable) correspond to multioutput vector indices in X
+            1D np.npdarray, if y univariate (one dimension)
+            otherwise, same type as y passed in fit
         """
         self.check_is_fitted()
 
+        # enter vectorized mode if needed
+        if self._is_vectorized:
+            return self._vectorize("predict", X=X)
+
         # boilerplate input checks for predict-like methods
         X = self._check_convert_X_for_predict(X)
 
-        return self._predict(X)
+        # call internal _predict, convert output
+        y_pred_inner = self._predict(X)
+        y_pred = self._convert_output_y(y_pred_inner)
+        return y_pred
 
-    def score(self, X, y) -> float:
+    def score(self, X, y, multioutput="uniform_average") -> float:
         """Scores predicted labels against ground truth labels on X.
 
         Parameters
         ----------
-        X : 3D np.array (any number of dimensions, equal length series)
-                of shape [n_instances, n_dimensions, series_length]
-            or 2D np.array (univariate, equal length series)
-                of shape [n_instances, series_length]
-            or pd.DataFrame with each column a dimension, each cell a pd.Series
-                (any number of dimensions, equal or unequal length series)
+        X : sktime compatible time series panel data container, Panel scitype, e.g.,
+            pd-multiindex: pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
+            numpy3D: 3D np.array (any number of dimensions, equal length series)
+            of shape [n_instances, n_dimensions, series_length]
             or of any other supported Panel mtype
-                for list of mtypes, see datatypes.SCITYPE_REGISTER
-                for specifications, see examples/AA_datatypes_and_datasets.ipynb
-        y : 1D np.array of float, of shape [n_instances] - regression labels (gnr truth)
-            indices correspond to instance indices in X
+            for list of mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
+        y : 2D np.array of int, of shape [n_instances, n_dimensions] - regression labels
+            for fitting indices correspond to instance indices in X
+            or 1D np.array of int, of shape [n_instances] - regression labels for
+            fitting indices correspond to instance indices in X
+        multioutput : str, optional (default="uniform_average")
+            {"raw_values", "uniform_average", "variance_weighted"}, array-like of shape
+            (n_outputs,) or None, default="uniform_average".
+            Defines aggregating of multiple output scores. Array-like value defines
+            weights used to average scores.
 
         Returns
         -------
@@ -233,7 +296,47 @@ class BaseRegressor(BaseEstimator, ABC):
 
         self.check_is_fitted()
 
-        return r2_score(y, self.predict(X), normalize=True)
+        return r2_score(y, self.predict(X), normalize=True, multioutput=multioutput)
+
+    def _vectorize(self, methodname, **kwargs):
+        """Vectorized/iterated loop over method of BaseClassifier.
+
+        Uses regressors_ attribute to store one regressor per loop index.
+        """
+        y = kwargs.get("y")
+        X = kwargs.get("X")
+        if X is not None:
+            kwargs.pop("X")
+        if y is not None:
+            kwargs.pop("y")
+            self._y_vec = y
+        regressors_ = self._y_vec.vectorize_est(
+            self,
+            method="clone",
+        )
+        if methodname == "fit":
+            self.regressors_ = self._y_vec.vectorize_est(
+                regressors_,
+                method=methodname,
+                args={"y": y},
+                X=X,
+            )
+            return self
+        else:
+            regressors_ = self.regressors_
+            y_pred = self._y_vec.vectorize_est(
+                regressors_,
+                method=methodname,
+                # return_type="list",
+                X=X,
+                args={"y": y} if y is not None else {},
+                **kwargs,  # contains X inside
+            )
+            y_pred = pd.DataFrame(
+                {str(i): y_pred[col].values[0] for i, col in enumerate(y_pred.columns)}
+            )
+            return y_pred
+        # add code for score
 
     @abstractmethod
     def _fit(self, X, y):
@@ -245,13 +348,21 @@ class BaseRegressor(BaseEstimator, ABC):
         ----------
         X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
             if self.get_tag("X_inner_mtype") = "numpy3D":
-                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "pd-multiindex":
+            pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
             if self.get_tag("X_inner_mtype") = "nested_univ":
-                pd.DataFrame with each column a dimension, each cell a pd.Series
+            pd.DataFrame with each column a dimension, each cell a pd.Series
             for list of other mtypes, see datatypes.SCITYPE_REGISTER
             for specifications, see examples/AA_datatypes_and_datasets.ipynb
-        y : 1D np.array of float, of shape [n_instances] - regression labels for fitting
-            indices correspond to instance indices in X
+        y : guaranteed to be of a type in self.get_tag("y_inner_mtype")
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            class labels for fitting
+            if self.get_tag("capaility:multioutput") = False, guaranteed to be 1D
+            if self.get_tag("capaility:multioutput") = True, guaranteed to be 2D
 
         Returns
         -------
@@ -310,13 +421,105 @@ class BaseRegressor(BaseEstimator, ABC):
         missing = X_metadata["has_nans"]
         multivariate = not X_metadata["is_univariate"]
         unequal = not X_metadata["is_equal_length"]
-        # Check this regressor can handle characteristics
+        X_mtype = X_metadata["mtype"]
+        # Check this classifier can handle characteristics
         self._check_capabilities(missing, multivariate, unequal)
-        # Convert data as dictated by the regressor tags
-        X = self._convert_X(X)
+        # Convert data as dictated by the classifier tags
+        X = self._convert_X(X, X_mtype=X_mtype)
 
         return X
 
+    def _check_y(self, y=None):
+        """Check and coerce X/y for fit/transform functions.
+
+        Parameters
+        ----------
+        y : pd.DataFrame, pd.Series or np.ndarray
+
+        Returns
+        -------
+        y : object of sktime compatible time series type
+            can be Series, Panel, Hierarchical
+        y_metadata : dict
+            metadata of y, retured by check_is_scitype
+        """
+        if y is None:
+            return None
+
+        capa_multioutput = self.get_tag("capability:multioutput")
+        y_inner_mtype = self.get_tag("y_inner_mtype")
+
+        y_valid, y_msg, y_metadata = check_is_scitype(
+            y, "Table", return_metadata=["is_univariate"]
+        )
+
+        if not y_valid:
+            allowed_msg = (
+                f"In regression, y must be of a supported type, "
+                f"for instance 1D or 2D numpy arrays, pd.DataFrame, or pd.Series. "
+                f"Allowed compatible mtype format specifications are:"
+                f" {MTYPE_LIST_TABLE} ."
+            )
+            check_is_error_msg(
+                y_msg, var_name="y", allowed_msg=allowed_msg, raise_exception=True
+            )
+
+        y_uni = y_metadata["is_univariate"]
+        y_mtype = y_metadata["mtype"]
+
+        requires_vectorization = not capa_multioutput and not y_uni
+
+        if requires_vectorization:
+            y_df = convert(
+                y,
+                from_type=y_mtype,
+                to_type="pd_DataFrame_Table",
+                as_scitype="Table",
+                store=self._converter_store_y,
+            )
+            y_vec = VectorizedDF([y_df], iterate_cols=True)
+            return y_vec, y_metadata
+
+        y_inner = convert(
+            y,
+            from_type=y_mtype,
+            to_type=y_inner_mtype,
+            as_scitype="Table",
+            store=self._converter_store_y,
+        )
+
+        return y_inner, y_metadata
+
+    def _convert_output_y(self, y):
+        """Convert output y to original format.
+
+        Parameters
+        ----------
+        y : np.ndarray or pd.DataFrame
+
+        Returns
+        -------
+        y : np.ndarray or pd.DataFrame
+        """
+        # for consistency with legacy behaviour:
+        # output is coerced to numpy1D in case of univariate output
+        if not self._y_metadata["is_univariate"]:
+            output_mtype = self._y_metadata["mtype"]
+            converter_store = self._converter_store_y
+        else:
+            output_mtype = "numpy1D"
+            converter_store = None
+
+        y = convert(
+            y,
+            from_type=self.get_tag("y_inner_mtype"),
+            to_type=output_mtype,
+            as_scitype="Table",
+            store=converter_store,
+            store_behaviour="freeze",
+        )
+        return y
+
     def _check_capabilities(self, missing, multivariate, unequal):
         """Check whether this regressor can handle the data characteristics.
 
@@ -363,24 +566,25 @@ class BaseRegressor(BaseEstimator, ABC):
             else:
                 raise ValueError(msg)
 
-    def _convert_X(self, X):
+    def _convert_X(self, X, X_mtype):
         """Convert equal length series from DataFrame to numpy array or vice versa.
 
         Parameters
         ----------
-        self : this regressor
-        X : pd.DataFrame or np.ndarray. Input attribute data
+        X : input data for the classifier
+        X_mtype : str, a Panel mtype string, e.g., "pd_multiindex", "numpy3D"
 
         Returns
         -------
         X : input X converted to type in "X_inner_mtype" tag
-                usually a pd.DataFrame (nested) or 3D np.ndarray
+            usually a pd.DataFrame (nested) or 3D np.ndarray
             Checked and possibly converted input data
         """
         inner_type = self.get_tag("X_inner_mtype")
         # convert pd.DataFrame
-        X = convert_to(
+        X = convert(
             X,
+            from_type=X_mtype,
             to_type=inner_type,
             as_scitype="Panel",
         )
@@ -421,7 +625,7 @@ def _check_regressor_input(
     )
     # raise informative error message if X is in wrong format
     allowed_msg = (
-        f"Allowed scitypes for classifiers are Panel mtypes, "
+        f"Allowed scitypes for regressors are Panel mtypes, "
         f"for instance a pandas.DataFrame with MultiIndex and last(-1) "
         f"level an sktime compatible time index. "
         f"Allowed compatible mtype format specifications are: {MTYPE_LIST_PANEL} ."
@@ -441,12 +645,10 @@ def _check_regressor_input(
     # Check y if passed
     if y is not None:
         # Check y valid input
-        if not isinstance(y, (pd.Series, np.ndarray)) and not (
-            isinstance(y, pd.DataFrame) and y.shape[1] == 1
-        ):
+        if not isinstance(y, (pd.Series, pd.DataFrame, np.ndarray)):
             raise ValueError(
-                f"y must be a np.array, pd.Series or a 1-D pd.DataFrame,"
-                f"but found type: {type(y)}"
+                "y must be a np.array or a pd.Series or pd.DataFrame, but found ",
+                f"type: {type(y)}",
             )
         # Check matching number of labels
         n_labels = y.shape[0]
@@ -456,9 +658,9 @@ def _check_regressor_input(
                 f"{n_labels}"
             )
         if isinstance(y, np.ndarray):
-            if y.ndim > 1:
+            if y.ndim > 2:
                 raise ValueError(
-                    f"np.ndarray y must be 1-dimensional, "
+                    f"np.ndarray y must be 1-dimensional or 2-dimensional, "
                     f"but found {y.ndim} dimensions"
                 )
     return X_metadata
diff --git a/sktime/split/sameloc.py b/sktime/split/sameloc.py
index 040c748d3..c94e1bb6c 100644
--- a/sktime/split/sameloc.py
+++ b/sktime/split/sameloc.py
@@ -141,10 +141,7 @@ class SameLocSplitter(BaseSplitter):
             `create_test_instance` uses the first (or only) dictionary in `params`
         """
         from sktime.datasets import load_airline
-        from sktime.forecasting.model_selection import (
-            ExpandingWindowSplitter,
-            SingleWindowSplitter,
-        )
+        from sktime.split import ExpandingWindowSplitter, SingleWindowSplitter
 
         y = load_airline()
         y_temp = y[:60]
diff --git a/sktime/transformations/base.py b/sktime/transformations/base.py
index 0dff77c85..860eb91a3 100644
--- a/sktime/transformations/base.py
+++ b/sktime/transformations/base.py
@@ -1186,10 +1186,20 @@ class BaseTransformer(BaseEstimator):
         #   skipped for output_scitype = "Primitives"
         #       since then the output always is a pd.DataFrame
         if case == "case 2: higher scitype supported" and output_scitype == "Series":
-            Xt = convert_to(
-                Xt,
-                to_type=["pd-multiindex", "numpy3D", "df-list", "pd_multiindex_hier"],
-            )
+            if self.get_tags()["scitype:transform-input"] == "Panel":
+                # Conversion from Series to Panel done for being compatible with
+                # algorithm. Thus, the returned Series should stay a Series.
+                pass
+            else:
+                Xt = convert_to(
+                    Xt,
+                    to_type=[
+                        "pd-multiindex",
+                        "numpy3D",
+                        "df-list",
+                        "pd_multiindex_hier",
+                    ],
+                )
             Xt = convert_to_scitype(Xt, to_scitype=X_input_scitype)
 
         # now, in all cases, Xt is in the right scitype,
@@ -1203,7 +1213,6 @@ class BaseTransformer(BaseEstimator):
         if output_scitype == "Series":
             # output mtype is input mtype
             X_output_mtype = X_input_mtype
-
             # exception to this: if the transformer outputs multivariate series,
             #   we cannot convert back to pd.Series, do pd.DataFrame instead then
             #   this happens only for Series, not Panel
@@ -1229,6 +1238,13 @@ class BaseTransformer(BaseEstimator):
                     )
                 if X_input_mtype == "pd.Series" and not metadata["is_univariate"]:
                     X_output_mtype = "pd.DataFrame"
+            elif self.get_tags()["scitype:transform-input"] == "Panel":
+                # Input has always to be Panel
+                X_output_mtype = "pd.DataFrame"
+            else:
+                # Input can be Panel or Hierachical, since it is supported
+                # by the used mtype
+                output_scitype = X_input_scitype
                 # Xt_mtype = metadata["mtype"]
             # else:
             #     Xt_mtype = X_input_mtype
@@ -1241,10 +1257,10 @@ class BaseTransformer(BaseEstimator):
             #     store=_converter_store_X,
             #     store_behaviour="freeze",
             # )
-            Xt = convert_to(
+            return convert_to(
                 Xt,
                 to_type=X_output_mtype,
-                as_scitype=X_input_scitype,
+                as_scitype=output_scitype,
                 store=_converter_store_X,
                 store_behaviour="freeze",
             )
@@ -1260,7 +1276,7 @@ class BaseTransformer(BaseEstimator):
                 # else this is only zeros and should be reset to RangeIndex
                 else:
                     Xt = Xt.reset_index(drop=True)
-            Xt = convert_to(
+            return convert_to(
                 Xt,
                 to_type="pd_DataFrame_Table",
                 as_scitype="Table",
diff --git a/sktime/transformations/merger.py b/sktime/transformations/merger.py
new file mode 100644
index 000000000..10c7024fb
--- /dev/null
+++ b/sktime/transformations/merger.py
@@ -0,0 +1,134 @@
+# copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
+"""Implements a merger for panel data."""
+import numpy as np
+
+__author__ = ["benHeid"]
+
+
+from sktime.transformations.base import BaseTransformer
+
+
+class Merger(BaseTransformer):
+    """Aggregates Panel data containing overlapping windows of one time series.
+
+    The input data contains multiple overlapping time series elements that could
+    arranged as follows:
+    xxxx.....
+    .xxxx....
+    ..xxxx...
+    ...xxxx..
+    ....xxxx.
+    .....xxxx
+    ......xxxx
+    .......xxxx
+    ........xxxx
+    .........xxxx
+    The merger aggregates the data by aligning the time series windows as shown above
+    and applying a aggregation function to the overlapping data points.
+    The aggregation function can be one of "mean" or "median". I.e., the `mean` or
+    `median` of each column is calculated, resulting in a univariate time series.
+
+    Parameters
+    ----------
+    method : {`median`, `mean`}, default="median"
+        The method to use for aggregation. Can be one of "mean" or "median".
+    stride : int, default=0
+        The stride to use for the aggregation. The stride determines the number of
+        shifts between consecutive instances. A stride of 0 means no shift. A
+        stride of 1 means that the time series is aggregated as above.
+
+    Examples
+    --------
+    >>> from sktime.transformations.merger import Merger
+    >>> from sktime.utils._testing.panel import _make_panel
+    >>> y = _make_panel(n_instances=10, n_columns=3, n_timepoints=5)
+    >>> result = Merger(method="median").fit_transform(y)
+    >>> result.shape
+    (3, 5)
+
+    >>> from sktime.transformations.merger import Merger
+    >>> from sktime.utils._testing.panel import _make_panel
+    >>> y = _make_panel(n_instances=10, n_columns=3, n_timepoints=5)
+    >>> result = Merger(method="median", stride=1).fit_transform(y)
+    >>> result.shape
+    (14, 3)
+    """
+
+    _tags = {
+        "scitype:transform-input": "Panel",
+        "scitype:transform-output": "Series",
+        "X_inner_mtype": "numpy3D",
+        "fit_is_empty": True,
+    }
+
+    def __init__(self, method="median", stride=0):
+        self.stride = stride
+        if method not in ["median", "mean"]:
+            raise ValueError(f"{method} must be 'mean' or 'median'.")
+        self.method = method
+        super().__init__()
+
+    def _transform(self, X=None, y=None):
+        """Merge the Panel data by aligning them temporally.
+
+        Parameters
+        ----------
+        X : pd.DataFrame
+            The input panel data.
+        y : pd.Series
+            ignored
+
+        Returns
+        -------
+        returns a single time series
+        """
+        horizon = X.shape[-1]
+
+        if self.method == "mean":
+            result = np.nanmean(self._align_temporal(horizon, X), axis=0)
+        elif self.method == "median":
+            result = np.nanmedian(self._align_temporal(horizon, X), axis=0)
+        else:
+            raise ValueError(f"{self.method} must be 'mean' or 'median'.")
+        return result
+
+    def _align_temporal(self, horizon, x):
+        x = x.astype(float)
+        if self.stride > 0:
+            x = np.insert(
+                x, np.arange(1, x.shape[0]).repeat(self.stride - 1), np.nan, axis=0
+            )
+        elif self.stride == 0:
+            return x
+        r = []
+        for i in range(horizon):
+            _res = np.concatenate(
+                [
+                    np.full(fill_value=np.nan, shape=(i, x.shape[1])),
+                    x[:, :, i],
+                    np.full(fill_value=np.nan, shape=((horizon - 1 - i), x.shape[1])),
+                ]
+            )
+            r.append(_res)
+        return np.stack(r)
+
+    @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`
+        """
+        return [{"method": "mean"}, {"method": "median"}]


Test patch:
diff --git a/sktime/classification/tests/test_all_classifiers.py b/sktime/classification/tests/test_all_classifiers.py
index d84533eac..e8f362266 100644
--- a/sktime/classification/tests/test_all_classifiers.py
+++ b/sktime/classification/tests/test_all_classifiers.py
@@ -4,6 +4,7 @@ __author__ = ["mloning", "TonyBagnall", "fkiraly"]
 
 
 import numpy as np
+import pandas as pd
 import pytest
 
 from sktime.classification.tests._expected_outputs import (
@@ -211,3 +212,25 @@ class TestAllClassifiers(ClassifierFixtureGenerator, QuickTester):
 
         with pytest.warns(UserWarning, match=error_msg):
             estimator_instance.fit(X, y)
+
+    def test_multioutput(self, estimator_instance):
+        """Test multioutput classification for all classifiers.
+
+        All classifiers should follow the same interface,
+        those that do not genuinely should vectorize/broadcast over y.
+        """
+        n_instances = 20
+        X, y = make_classification_problem(n_instances=n_instances)
+        y_mult = pd.DataFrame({"a": y, "b": y})
+
+        estimator_instance.fit(X, y_mult)
+        y_pred = estimator_instance.predict(X)
+
+        assert isinstance(y_pred, pd.DataFrame)
+        assert y_pred.shape == y_mult.shape
+
+        vectorized = estimator_instance.get_tag("capability:multioutput")
+        if vectorized:
+            assert hasattr(estimator_instance, "classifiers_")
+            assert isinstance(estimator_instance.classifiers_, pd.DataFrame)
+            assert estimator_instance.classifiers_.shape == (1, 2)
diff --git a/sktime/classification/tests/test_base.py b/sktime/classification/tests/test_base.py
index 9713be6f0..525c29125 100644
--- a/sktime/classification/tests/test_base.py
+++ b/sktime/classification/tests/test_base.py
@@ -120,8 +120,6 @@ class _DummyConvertPandas(BaseClassifier):
 multivariate_message = r"multivariate series"
 missing_message = r"missing values"
 unequal_message = r"unequal length series"
-incorrect_X_data_structure = r"must be a np.array or a pd.Series"
-incorrect_y_data_structure = r"must be 1-dimensional"
 
 
 def test_base_classifier_fit():
@@ -144,8 +142,11 @@ def test_base_classifier_fit():
     test_X3 = _create_example_dataframe(cases=cases, dimensions=1, length=length)
     test_X4 = _create_example_dataframe(cases=cases, dimensions=3, length=length)
     test_y1 = np.random.randint(0, 2, size=(cases))
+    test_y2 = pd.DataFrame({"0": [1] * cases, "1": [0] * cases})
     result = dummy.fit(test_X1, test_y1)
     assert result is dummy
+    result = dummy.fit(test_X3, test_y2)
+    assert result is dummy
     with pytest.raises(ValueError, match=multivariate_message):
         result = dummy.fit(test_X2, test_y1)
     assert result is dummy
@@ -154,15 +155,6 @@ def test_base_classifier_fit():
     with pytest.raises(ValueError, match=multivariate_message):
         result = dummy.fit(test_X4, test_y1)
     assert result is dummy
-    # Raise a specific error if y is in a 2D matrix (1,cases)
-    test_y2 = np.array([test_y1])
-    # What if y is in a 2D matrix (cases,1)?
-    test_y2 = np.array([test_y1]).transpose()
-    with pytest.raises(ValueError, match=incorrect_y_data_structure):
-        result = dummy.fit(test_X1, test_y2)
-    # Pass a data fram
-    with pytest.raises(ValueError, match=incorrect_X_data_structure):
-        result = dummy.fit(test_X1, test_X3)
 
 
 TF = [True, False]
@@ -228,26 +220,26 @@ def test_convert_input():
     test_X1 = np.random.uniform(-1, 1, size=(cases, length))
     test_X2 = np.random.uniform(-1, 1, size=(cases, 2, length))
     tester = _DummyClassifier()
-    tempX = tester._convert_X(test_X2)
+    tempX = tester._convert_X(test_X2, "numpy3D")
     assert tempX.shape[0] == cases and tempX.shape[1] == 2 and tempX.shape[2] == length
     instance_list = []
     for _ in range(0, cases):
         instance_list.append(pd.Series(np.random.randn(10)))
     test_X3 = _create_example_dataframe(cases=cases, dimensions=1, length=length)
     test_X4 = _create_example_dataframe(cases=cases, dimensions=3, length=length)
-    tempX = tester._convert_X(test_X3)
+    tempX = tester._convert_X(test_X3, "nested_univ")
     assert tempX.shape[0] == cases and tempX.shape[1] == 1 and tempX.shape[2] == length
-    tempX = tester._convert_X(test_X4)
+    tempX = tester._convert_X(test_X4, "nested_univ")
     assert tempX.shape[0] == cases and tempX.shape[1] == 3 and tempX.shape[2] == length
     tester = _DummyConvertPandas()
-    tempX = tester._convert_X(test_X2)
+    tempX = tester._convert_X(test_X2, "numpy3D")
     assert isinstance(tempX, pd.DataFrame)
     assert tempX.shape[0] == cases
     assert tempX.shape[1] == 2
     test_y1 = np.random.randint(0, 1, size=(cases))
     test_y1 = pd.Series(test_y1)
     tempX, tempY = _internal_convert(test_X1, test_y1)
-    assert isinstance(tempY, np.ndarray)
+    assert isinstance(tempY, pd.Series)
     assert isinstance(tempX, np.ndarray)
     assert tempX.ndim == 3
 
diff --git a/sktime/forecasting/compose/tests/test_pipeline.py b/sktime/forecasting/compose/tests/test_pipeline.py
index d3b650dca..299cd7229 100644
--- a/sktime/forecasting/compose/tests/test_pipeline.py
+++ b/sktime/forecasting/compose/tests/test_pipeline.py
@@ -485,7 +485,7 @@ def test_forecastx_logic():
 )
 def test_forecastx_fit_behavior():
     from sktime.forecasting.compose import ForecastX
-    from sktime.forecasting.model_selection import temporal_train_test_split
+    from sktime.split import temporal_train_test_split
 
     y, X = load_longley()
     y_train, y_test, X_train, X_test = temporal_train_test_split(y, X)
diff --git a/sktime/forecasting/trend/tests/test_pwl_trend.py b/sktime/forecasting/trend/tests/test_pwl_trend.py
index a1e4d938c..b5ab77109 100644
--- a/sktime/forecasting/trend/tests/test_pwl_trend.py
+++ b/sktime/forecasting/trend/tests/test_pwl_trend.py
@@ -11,11 +11,11 @@ import pytest
 
 from sktime.datasets import load_airline
 from sktime.forecasting.base import ForecastingHorizon
-from sktime.forecasting.model_selection import temporal_train_test_split
 from sktime.forecasting.trend import (
     PolynomialTrendForecaster,
     ProphetPiecewiseLinearTrendForecaster,
 )
+from sktime.split import temporal_train_test_split
 from sktime.tests.test_switch import run_test_for_class
 
 
diff --git a/sktime/regression/tests/test_all_regressors.py b/sktime/regression/tests/test_all_regressors.py
index 3dfea1768..c3e8b1c4a 100644
--- a/sktime/regression/tests/test_all_regressors.py
+++ b/sktime/regression/tests/test_all_regressors.py
@@ -4,10 +4,12 @@ __author__ = ["mloning", "TonyBagnall", "fkiraly"]
 
 
 import numpy as np
+import pandas as pd
 import pytest
 
 from sktime.datatypes import check_is_scitype
 from sktime.tests.test_all_estimators import BaseFixtureGenerator, QuickTester
+from sktime.utils._testing.panel import make_regression_problem
 from sktime.utils._testing.scenarios_classification import (
     ClassifierFitPredictMultivariate,
 )
@@ -80,3 +82,25 @@ class TestAllRegressors(RegressorFixtureGenerator, QuickTester):
         assert isinstance(y_pred, np.ndarray)
         assert y_pred.shape == (X_new_instances,)
         assert np.issubdtype(y_pred.dtype, np.floating)
+
+    def test_multioutput(self, estimator_instance):
+        """Test multioutput regression for all classifiers.
+
+        All classifiers should follow the same interface,
+        those that do not genuinely should vectorize/broadcast over y.
+        """
+        n_instances = 20
+        X, y = make_regression_problem(n_instances=n_instances)
+        y_mult = pd.DataFrame({"a": y, "b": y})
+
+        estimator_instance.fit(X, y_mult)
+        y_pred = estimator_instance.predict(X)
+
+        assert isinstance(y_pred, pd.DataFrame)
+        assert y_pred.shape == y_mult.shape
+
+        vectorized = estimator_instance.get_tag("capability:multioutput")
+        if vectorized:
+            assert hasattr(estimator_instance, "regressors_")
+            assert isinstance(estimator_instance.regressors_, pd.DataFrame)
+            assert estimator_instance.regressors_.shape == (1, 2)
diff --git a/sktime/split/temporal_train_test_split.py b/sktime/split/temporal_train_test_split.py
index 2509b0ec9..55e7e9195 100644
--- a/sktime/split/temporal_train_test_split.py
+++ b/sktime/split/temporal_train_test_split.py
@@ -91,8 +91,8 @@ def temporal_train_test_split(
 
     Examples
     --------
-    >>> from sktime.forecasting.model_selection import temporal_train_test_split
     >>> from sktime.datasets import load_airline, load_osuleaf
+    >>> from sktime.split import temporal_train_test_split
     >>> from sktime.utils._testing.panel import _make_panel
     >>> # univariate time series
     >>> y = load_airline()
diff --git a/sktime/split/testplustrain.py b/sktime/split/testplustrain.py
index 73981a9c6..a803a52c9 100644
--- a/sktime/split/testplustrain.py
+++ b/sktime/split/testplustrain.py
@@ -131,10 +131,7 @@ class TestPlusTrainSplitter(BaseSplitter):
             `MyClass(**params)` or `MyClass(**params[i])` creates a valid test instance.
             `create_test_instance` uses the first (or only) dictionary in `params`
         """
-        from sktime.forecasting.model_selection import (
-            ExpandingWindowSplitter,
-            SingleWindowSplitter,
-        )
+        from sktime.split import ExpandingWindowSplitter, SingleWindowSplitter
 
         cv_1 = ExpandingWindowSplitter(fh=[2, 4], initial_window=24, step_length=12)
         cv_2 = SingleWindowSplitter(fh=[2, 4], window_length=24)
diff --git a/sktime/transformations/tests/test_all_transformers.py b/sktime/transformations/tests/test_all_transformers.py
index 493fb63d3..b6a9699ac 100644
--- a/sktime/transformations/tests/test_all_transformers.py
+++ b/sktime/transformations/tests/test_all_transformers.py
@@ -78,6 +78,8 @@ class TestAllTransformers(TransformerFixtureGenerator, QuickTester):
                 return "Panel"
             if X_scitype in ["Panel", "Hierarchical"]:
                 return "Hierarchical"
+        if trafo_input == "Panel" and trafo_output == "Series":
+            return "Series"
 
     def test_fit_transform_output(self, estimator_instance, scenario):
         """Test that transform output is of expected scitype."""