Repository: sktime/sktime

Issue:
[BUG] ForecastingBenchmark giving an error when the dataloader returns the tuple (y, X)
**Describe the bug**
<!--
A clear and concise description of what the bug is.
--> When using the ForecastingBenchmark with a data loader that returns a tuple (y, X) (where y is the endogenous data and X is the exogenous data), I encounter the following error:
`TypeError: Expected y dtype ['Series', 'Panel', 'Hierarchical']. Got <class 'tuple'> instead.`

**To Reproduce**
```python

from sktime.benchmarking.forecasting import ForecastingBenchmark
from sktime.performance_metrics.forecasting import MeanAbsolutePercentageError
from sktime.forecasting.naive import NaiveForecaster
from sktime.forecasting.exp_smoothing import ExponentialSmoothing
from sktime.forecasting.theta import ThetaForecaster
from sktime.split import ExpandingWindowSplitter

def load():
    return (y, X) # both are Hierarchical type dataframes

cv = ExpandingWindowSplitter(initial_window=60, step_length=24, fh=np.arange(1, 25))

benchmark = ForecastingBenchmark()

benchmark.add_estimator(NaiveForecaster(strategy='mean'))
benchmark.add_estimator(ExponentialSmoothing())
benchmark.add_estimator(ThetaForecaster())

benchmark.add_task(dataset_loader=load, cv_splitter=cv, scorers=[MeanAbsolutePercentageError()], task_id="M5_mape")

benchmark.run(output_file= 'results.txt')
```
![image](https://github.com/user-attachments/assets/e3b48aa3-d76a-4452-9a8a-3d5c5051df54)


<!-- Thanks for contributing! -->


Discussion:

Patch:
diff --git a/.all-contributorsrc b/.all-contributorsrc
index e44e1806a..7bad0c402 100644
--- a/.all-contributorsrc
+++ b/.all-contributorsrc
@@ -2652,8 +2652,10 @@
       "avatar_url": "https://avatars.githubusercontent.com/u/37468543?v=4",
       "profile": "https://github.com/meraldoantonio",
       "contributions": [
+        "bug",
+        "code",
         "doc",
-        "bug"
+        "test"
       ]
     },
     {
diff --git a/docs/source/api_reference/datasets.rst b/docs/source/api_reference/datasets.rst
index 6bcd61e83..9184d70dd 100644
--- a/docs/source/api_reference/datasets.rst
+++ b/docs/source/api_reference/datasets.rst
@@ -8,6 +8,7 @@ The ``datasets`` module contains:
 * loaders which fetch datasets from data repositories on the internet,
   and retrieve them as in-memory datasets in ``sktime`` compatible formats
 * loaders which fetch an individual dataset, usually for illustration purposes
+* toy data generators for didactic and illustrative purposes
 * utilities to write to, and load from, time series specific file formats
 
 Loaders from dataset repositories
@@ -64,7 +65,6 @@ Single time series
 Panels of time series
 ^^^^^^^^^^^^^^^^^^^^^
 
-
 .. automodule:: sktime.datasets
     :no-members:
     :no-inherited-members:
@@ -85,6 +85,25 @@ Panels of time series
     load_osuleaf
 
 
+Toy data generators
+-------------------
+
+Hierarchical time series data
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. automodule:: sktime.datasets
+    :no-members:
+    :no-inherited-members:
+
+.. currentmodule:: sktime.datasets
+
+.. autosummary::
+    :toctree: auto_generated/
+    :template: function.rst
+
+    load_hierarchical_sales_toydata
+
+
 Loading from and writing to files
 ---------------------------------
 
diff --git a/docs/source/api_reference/forecasting.rst b/docs/source/api_reference/forecasting.rst
index 73d8af577..13aa25386 100644
--- a/docs/source/api_reference/forecasting.rst
+++ b/docs/source/api_reference/forecasting.rst
@@ -320,6 +320,14 @@ All "ARIMA" and "Auto-ARIMA" models below include SARIMAX capability.
 
     VAR
 
+.. currentmodule:: sktime.forecasting.var_reduce
+
+.. autosummary::
+    :toctree: auto_generated/
+    :template: class.rst
+
+    VARReduce
+
 .. currentmodule:: sktime.forecasting.varmax
 
 .. autosummary::
diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst
index 66bcdb475..3b9fea703 100644
--- a/docs/source/changelog.rst
+++ b/docs/source/changelog.rst
@@ -19,6 +19,26 @@ For upcoming changes and next releases, see our `milestones <https://github.com/
 For our long-term plan, see our :ref:`roadmap`.
 
 
+Version 0.31.2 - 2024-08-13
+---------------------------
+
+Hotfix release for using ``make_reduction`` with not fully ``sklearn`` compliant
+tabular regressors such as from ``catboost``.
+
+For last non-maintenance content updates, see 0.31.1.
+
+Contents
+~~~~~~~~
+
+* [BUG] fix ``make_reduction`` type inference for non-sklearn estimators
+
+Notes
+^^^^^
+
+This is a hotfix for 0.31.1 release, fixing a regression. This release is not contained
+in the 0.32.0 or 0.32.1 releases.
+
+
 Version 0.32.1 - 2024-08-12
 ---------------------------
 
diff --git a/pyproject.toml b/pyproject.toml
index ed053590d..e8275d0c1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -96,6 +96,7 @@ all_extras = [
   'numba<0.61,>=0.53',
   'optuna<3.7',
   'pmdarima!=1.8.1,<3.0.0,>=1.8; python_version < "3.12"',
+  "polars[pandas]",
   'prophet>=1.1; python_version < "3.12"',
   "pycatch22<0.4.6",
   'pyod>=0.8; python_version < "3.11"',
@@ -140,6 +141,7 @@ all_extras_pandas2 = [
   'numba<0.61,>=0.53',
   'optuna<3.7',
   'pmdarima!=1.8.1,<3.0.0,>=1.8; python_version < "3.12"',
+  "polars[pandas]",
   'prophet>=1.1; python_version < "3.12"',
   "pycatch22<0.4.6",
   'pyod>=0.8; python_version < "3.11"',
diff --git a/sktime/base/_base.py b/sktime/base/_base.py
index c85c4ae4f..e4551a3ae 100644
--- a/sktime/base/_base.py
+++ b/sktime/base/_base.py
@@ -205,7 +205,7 @@ class BaseObject(_HTMLDocumentationLinkMixin, _BaseObject):
         from packaging.requirements import Requirement
 
         # pypi package names of soft dependencies that are not numpy 2 compatibleS
-        NOT_NP2_COMPATIBLE = ["prophet"]
+        NOT_NP2_COMPATIBLE = ["prophet", "pmdarima"]
 
         softdeps = self.get_class_tag("python_dependencies", [])
         if softdeps is None:
diff --git a/sktime/benchmarking/forecasting.py b/sktime/benchmarking/forecasting.py
index 3db224d7c..47eb852c5 100644
--- a/sktime/benchmarking/forecasting.py
+++ b/sktime/benchmarking/forecasting.py
@@ -72,14 +72,27 @@ def forecasting_validation(
     """
     y = dataset_loader()
     results = {}
-    scores_df = evaluate(
-        forecaster=estimator,
-        y=y,
-        cv=cv_splitter,
-        scoring=scorers,
-        backend=backend,
-        backend_params=backend_params,
-    )
+    if isinstance(y, tuple):
+        y, X = y
+        scores_df = evaluate(
+            forecaster=estimator,
+            y=y,
+            X=X,
+            cv=cv_splitter,
+            scoring=scorers,
+            backend=backend,
+            backend_params=backend_params,
+        )
+    else:
+        scores_df = evaluate(
+            forecaster=estimator,
+            y=y,
+            cv=cv_splitter,
+            scoring=scorers,
+            backend=backend,
+            backend_params=backend_params,
+        )
+
     for scorer in scorers:
         scorer_name = scorer.name
         for ix, row in scores_df.iterrows():
diff --git a/sktime/datasets/__init__.py b/sktime/datasets/__init__.py
index a2a11171e..19a0933ca 100644
--- a/sktime/datasets/__init__.py
+++ b/sktime/datasets/__init__.py
@@ -21,6 +21,7 @@ __all__ = [
     "load_electric_devices_segmentation",
     "load_acsf1",
     "load_macroeconomic",
+    "load_hierarchical_sales_toydata",
     "generate_example_long_table",
     "load_from_arff_to_dataframe",
     "load_from_long_to_dataframe",
@@ -49,6 +50,7 @@ from sktime.datasets._data_io import (
     make_multi_index_dataframe,
 )
 from sktime.datasets._fpp3_loaders import DATASET_NAMES_FPP3, load_fpp3
+from sktime.datasets._hierarchical_demo import load_hierarchical_sales_toydata
 from sktime.datasets._readers_writers.arff import (
     load_from_arff_to_dataframe,
     write_tabular_transformation_to_arff,
diff --git a/sktime/datasets/_hierarchical_demo.py b/sktime/datasets/_hierarchical_demo.py
new file mode 100644
index 000000000..0c57ae424
--- /dev/null
+++ b/sktime/datasets/_hierarchical_demo.py
@@ -0,0 +1,77 @@
+"""Demo data for hierarchical forecasting."""
+
+import numpy as np
+
+__author__ = ["marrov"]
+
+
+def load_hierarchical_sales_toydata():
+    """Return hierarchical sales toy data to demonstrate hierarchical forecasting.
+
+    Return data covers 5 years of monthly sales data for 2 product lines,
+    and 4 product groups.
+
+    Hierarchical structure:
+
+    - Product line 1: Food preparation
+        - Product group 1A: Hobs
+        - Product group 1B: Ovens
+    - Product line 2: Food preservation
+        - Product group 2A: Fridges
+        - Product group 2B: Freezers
+
+    Returns
+    -------
+    hierarchy : pd.DataFrame in pd_multiindex_hier mtype format
+        Product hierarchy with row MultiIndex "Product line", "Product group", "Date".
+        Column "Sales" contains total sales for each product group in monthly period.
+    """
+    # Get daily historic sales and rename columns and indexes according to hierarchy
+    from sktime.utils._testing.hierarchical import _make_hierarchical
+
+    n_years = 5
+    y = (
+        _make_hierarchical(
+            hierarchy_levels=(2, 4),
+            min_timepoints=365 * n_years,
+            max_timepoints=365 * n_years,
+            random_state=0,
+        )
+        .drop(
+            index=[
+                ("h0_0", "h1_2"),
+                ("h0_0", "h1_3"),
+                ("h0_1", "h1_0"),
+                ("h0_1", "h1_1"),
+            ]
+        )
+        .rename(
+            index={
+                "h0_0": "Food preparation",
+                "h0_1": "Food preservation",
+                "h1_0": "Hobs",
+                "h1_1": "Ovens",
+                "h1_2": "Fridges",
+                "h1_3": "Freezers",
+            }
+        )
+        .reset_index()
+        .rename(
+            columns={
+                "h0": "Product line",
+                "h1": "Product group",
+                "time": "Date",
+                "c0": "Sales",
+            }
+        )
+    )
+
+    # Set date as monthly as sales as int and aggregate date
+    y["Date"] = y["Date"].dt.to_period("M")
+    y = y.groupby(by=["Product line", "Product group", "Date"]).sum()
+
+    # Add noise to have different time series
+    noise = np.random.RandomState(seed=0).normal(1, 0.3, np.shape(y))
+    y = (y * noise).round(0)
+
+    return y
diff --git a/sktime/datatypes/_adapter/polars.py b/sktime/datatypes/_adapter/polars.py
index a14e3e1ed..257cb3652 100644
--- a/sktime/datatypes/_adapter/polars.py
+++ b/sktime/datatypes/_adapter/polars.py
@@ -3,6 +3,7 @@
 
 from sktime.datatypes._common import _req
 from sktime.datatypes._common import _ret as ret
+from sktime.datatypes._dtypekind import _get_feature_kind, _polars_dtype_to_kind
 
 
 def get_mi_cols(obj):
@@ -10,6 +11,52 @@ def get_mi_cols(obj):
     return [x for x in obj.columns if isinstance(x, str) and x.startswith("__index__")]
 
 
+def is_monotonically_increasing(obj):
+    """Check is polars frame columns(__index__) is monotonically increasing."""
+    index_cols = get_mi_cols(obj)
+
+    # Series with no index columns
+    if len(index_cols) == 0:
+        return True
+
+    # check for Series scitype
+    if len(index_cols) == 1:
+        if obj[index_cols[0]].is_sorted():
+            return True
+
+    # check for Panel and Hierarchical scitype
+    else:
+        import polars as pl
+
+        index_df = obj.with_columns(index_cols)
+        grouped = index_df.groupby(index_cols[:-1]).agg([pl.col(index_cols[-1])])
+        last_index_col = grouped.select([index_cols[-1]])
+        for val in last_index_col.iter_rows():
+            # iter rows returns a list of tuples
+            if not pl.Series(val[0]).is_sorted():
+                return False
+        return True
+
+    return False
+
+
+def _convert_period_index_to_datetime_index(obj):
+    """Convert PeriodIndex to DatatimeIndex as polars only supports DatetimeIndex."""
+    import pandas as pd
+
+    if isinstance(obj.index, pd.PeriodIndex):
+        obj.index = obj.index.to_timestamp(freq=obj.index.freq)
+
+    if isinstance(obj.index, pd.MultiIndex):
+        levels = obj.index.levels
+        if isinstance(levels[-1], pd.PeriodIndex):
+            new_levels = list(levels)
+            new_levels[-1] = new_levels[-1].to_timestamp(freq=new_levels[-1].freq)
+            obj.index = obj.index.set_levels(new_levels)
+
+    return obj
+
+
 def convert_pandas_to_polars(
     obj, schema_overrides=None, rechunk=True, nan_to_null=True, lazy=False
 ):
@@ -46,6 +93,10 @@ def convert_pandas_to_polars(
                 col[i] = f"__index__{name}"
         return col
 
+    obj = obj.copy()
+    obj = _convert_period_index_to_datetime_index(
+        obj
+    )  # Polars only supports DatetimeIndex
     index_names = obj.index.names
     index_names = polars_index_columns(index_names)
     obj.index.names = index_names
@@ -61,12 +112,14 @@ def convert_pandas_to_polars(
     return obj
 
 
-def convert_polars_to_pandas(obj):
+def convert_polars_to_pandas(obj, infer_freq=True):
     """Convert polars DataFrame to pandas DataFrame, preserving MultiIndex.
 
     Parameters
     ----------
     obj : polars.DataFrame, polars.LazyFrame
+    infer_freq : bool, optional (default=True)
+        Infer frequency and set freq attribute of DatetimeIndex and DatetimeIndex levels
 
     Returns
     -------
@@ -75,6 +128,7 @@ def convert_polars_to_pandas(obj):
         by converting columns with names __index__[indexname] to MultiIndex levels,
         and other columns identical to those of obj.
     """
+    import pandas as pd
     from polars.lazyframe.frame import LazyFrame
 
     # convert to DataFrame if LazyFrame
@@ -93,17 +147,31 @@ def convert_polars_to_pandas(obj):
                 names[i] = None
         return names
 
+    def set_freq(obj):
+        if isinstance(obj.index, pd.DatetimeIndex):
+            obj.index.freq = pd.infer_freq(obj.index)
+
+        if isinstance(obj.index, pd.MultiIndex):
+            levels = obj.index.levels
+            if isinstance(levels[-1], pd.DatetimeIndex):
+                obj.index.levels[-1].freq = pd.infer_freq(obj.index.levels[-1])
+        return obj
+
     pd_index_names = get_index_names(obj.columns)
 
     if len(pd_index_names) > 0:
         pl_index_names = get_mi_cols(obj)
         obj = obj.set_index(pl_index_names)
         obj.index.names = pd_index_names
+        if infer_freq:
+            obj = set_freq(obj)
 
     return obj
 
 
-def check_polars_frame(obj, return_metadata=False, var_name="obj", lazy=False):
+def check_polars_frame(
+    obj, return_metadata=False, var_name="obj", lazy=False, scitype="Table"
+):
     """Check polars frame, generic format."""
     import polars as pl
 
@@ -121,19 +189,76 @@ def check_polars_frame(obj, return_metadata=False, var_name="obj", lazy=False):
         return ret(False, msg, None, return_metadata)
 
     # we now know obj is a polars DataFrame or LazyFrame
+    if scitype not in ["Table", "Series", "Panel", "Hierarchical"]:
+        return RuntimeError(
+            'scitype arg of check_polars_frame must be one of strings "Table", '
+            f'"Series", "Panel", or "Hierarchical", but found {scitype}'
+        )
+
+    index_cols = []
+
+    if scitype in ["Series", "Panel", "Hierarchical"]:
+        index_cols = get_mi_cols(obj)
+        n_vars = len(index_cols)
+        scitypes_index = {
+            "Series": n_vars == 1 or n_vars == 0,
+            "Panel": n_vars == 2,
+            "Hierarchical": n_vars >= 3,
+        }
+
+        if not scitypes_index[scitype]:
+            cols_msg = (
+                f"{var_name} must have correct number of index columns for scitype, "
+                f"Series: 0 or 1, Panel: 2, Hierarchical: >= 3,"
+                f"found {len(index_cols)}, namely: {index_cols}"
+            )
+            return ret(False, cols_msg, None, return_metadata)
+
+        # check if index columns are monotonically increasing
+        if isinstance(obj, pl.DataFrame) and not is_monotonically_increasing(obj):
+            msg = (
+                f"The (time) index of {var_name} must be sorted monotonically "
+                f"increasing. Use {var_name}.sort() on columns representing "
+                f"index(__index__) to sort the index, or {var_name}.is_duplicated() "
+                f"to find duplicates."
+            )
+            return ret(False, msg, None, return_metadata)
+
+    # columns in polars are unique, no check required
+
     if _req("is_empty", return_metadata):
         metadata["is_empty"] = obj.width < 1
     if _req("is_univariate", return_metadata):
-        metadata["is_univariate"] = obj.width == 1
-    if _req("n_instances", return_metadata):
-        if hasattr(obj, "height"):
-            metadata["n_instances"] = obj.height
-        else:
-            metadata["n_instances"] = "NA"
+        metadata["is_univariate"] = obj.width - len(index_cols) == 1
     if _req("n_features", return_metadata):
-        metadata["n_features"] = obj.width
+        metadata["n_features"] = obj.width - len(index_cols)
     if _req("feature_names", return_metadata):
-        metadata["feature_names"] = obj.columns
+        feature_columns = [x for x in obj.columns if x not in index_cols]
+        metadata["feature_names"] = feature_columns
+    if _req("dtypekind_dfip", return_metadata):
+        index_cols_count = len(index_cols)
+        dtype_list = obj.dtypes[index_cols_count:]
+        metadata["dtypekind_dfip"] = _polars_dtype_to_kind(dtype_list)
+    if _req("feature_kind", return_metadata):
+        index_cols_count = len(index_cols)
+        dtype_list = obj.dtypes[index_cols_count:]
+        dtype_kind = _polars_dtype_to_kind(dtype_list)
+        metadata["feature_kind"] = _get_feature_kind(dtype_kind)
+
+    if scitype == "Table":
+        if _req("n_instances", return_metadata):
+            if hasattr(obj, "height"):
+                metadata["n_instances"] = obj.height
+            else:
+                metadata["n_instances"] = "NA"
+
+    if scitype in ["Panel", "Hierarchical"]:
+        if _req("n_instances", return_metadata):
+            if exp_type_str == "LazyFrame":
+                metadata["n_instances"] = "NA"
+            else:
+                instance_cols = index_cols[:-1]
+                metadata["n_instances"] = len(obj[instance_cols].unique())
 
     # check if there are any nans
     #   compute only if needed
diff --git a/sktime/datatypes/_dtypekind.py b/sktime/datatypes/_dtypekind.py
index 06e63504c..08dc3dc05 100644
--- a/sktime/datatypes/_dtypekind.py
+++ b/sktime/datatypes/_dtypekind.py
@@ -4,7 +4,7 @@ Enum values of DtypeKind are in accordance with the dataframe interchange protoc
 Utility functions to help convert dtypes to corresponding DtypeKind values.
 """
 
-__author__ = ["Abhay-Lejith"]
+__author__ = ["Abhay-Lejith", "pranavvp16"]
 
 from enum import IntEnum
 
@@ -76,6 +76,34 @@ def _pandas_dtype_to_kind(col_dtypes):
         )
 
 
+def _polars_dtype_to_kind(col_dtypes):
+    """Convert polars dtypes to enum DtypeKind values."""
+    import polars as pl
+
+    if isinstance(col_dtypes, list):
+        return [_polars_dtype_to_kind(x) for x in col_dtypes]
+
+    if col_dtypes.is_float():
+        return DtypeKind.FLOAT
+    elif col_dtypes.is_signed_integer():
+        return DtypeKind.INT
+    elif col_dtypes.is_unsigned_integer():
+        return DtypeKind.UINT
+    elif isinstance(col_dtypes, pl.Categorical) or isinstance(col_dtypes, pl.Object):
+        return DtypeKind.CATEGORICAL
+    elif isinstance(col_dtypes, pl.Boolean):
+        return DtypeKind.BOOL
+    elif col_dtypes.is_temporal():
+        return DtypeKind.DATETIME
+    elif isinstance(col_dtypes, pl.String):
+        return DtypeKind.STRING
+    else:
+        raise TypeError(
+            "Dtype of columns can be: INT, UINT, FLOAT, BOOL, STRING, "
+            f"DATETIME, CATEGORICAL. Found dtype: {col_dtypes}"
+        )
+
+
 # function for series scitype
 def _get_series_dtypekind(obj, mtype):
     if mtype in ["numpy", "xarray"]:
diff --git a/sktime/datatypes/_series/_check.py b/sktime/datatypes/_series/_check.py
index f80d1a0ff..888a38dda 100644
--- a/sktime/datatypes/_series/_check.py
+++ b/sktime/datatypes/_series/_check.py
@@ -384,6 +384,28 @@ if _check_soft_dependencies("dask", severity="none"):
 
     check_dict[("dask_series", "Series")] = check_dask_series
 
+
+if _check_soft_dependencies("polars", severity="none"):
+    from sktime.datatypes._adapter.polars import check_polars_frame
+
+    def check_polars_series(obj, return_metadata=False, var_name="obj"):
+        metadict = check_polars_frame(
+            obj=obj,
+            return_metadata=return_metadata,
+            var_name=var_name,
+            scitype="Series",
+        )
+
+        if isinstance(metadict, tuple) and metadict[0]:
+            # update dict with Series specific keys
+            if _req("is_equally_spaced", return_metadata):
+                metadict[2]["is_equally_spaced"] = "NA"
+
+        return metadict
+
+    check_dict[("pl.DataFrame", "Series")] = check_polars_series
+
+
 if _check_soft_dependencies("gluonts", severity="none"):
 
     def check_gluonTS_listDataset_series(obj, return_metadata=False, var_name="obj"):
diff --git a/sktime/datatypes/_series/_convert.py b/sktime/datatypes/_series/_convert.py
index 7e9e812e7..e1ac407ee 100644
--- a/sktime/datatypes/_series/_convert.py
+++ b/sktime/datatypes/_series/_convert.py
@@ -266,6 +266,57 @@ if _check_soft_dependencies("dask", severity="none"):
         "dask_series", "pd.DataFrame", convert_dict, mtype_universe=MTYPE_LIST_SERIES
     )
 
+
+if _check_soft_dependencies("polars", severity="none"):
+    from sktime.datatypes._adapter.polars import (
+        convert_pandas_to_polars,
+        convert_polars_to_pandas,
+    )
+
+    def convert_polars_to_uvs_as_series(obj, store=None):
+        pd_df = convert_polars_to_pandas(obj)
+        return convert_MvS_to_UvS_as_Series(pd_df, store=store)
+
+    convert_dict[("pl.DataFrame", "pd.Series", "Series")] = (
+        convert_polars_to_uvs_as_series
+    )
+
+    def convert_polars_to_mvs_as_series(obj, store=None):
+        return convert_polars_to_pandas(obj)
+
+    convert_dict[("pl.DataFrame", "pd.DataFrame", "Series")] = (
+        convert_polars_to_mvs_as_series
+    )
+
+    def convert_mvs_to_polars_as_series(obj, store=None):
+        return convert_pandas_to_polars(obj)
+
+    convert_dict[("pd.DataFrame", "pl.DataFrame", "Series")] = (
+        convert_mvs_to_polars_as_series
+    )
+
+    def convert_uvs_to_polars_as_series(obj, store=None):
+        return convert_pandas_to_polars(obj)
+
+    convert_dict[("pd.Series", "pl.DataFrame", "Series")] = (
+        convert_uvs_to_polars_as_series
+    )
+
+    def convert_polars_lazy_to_mvs_as_series(obj, store=None):
+        return convert_polars_to_pandas(obj)
+
+    convert_dict[("pl.LazyFrame", "pd.DataFrame", "Series")] = (
+        convert_polars_lazy_to_mvs_as_series
+    )
+
+    def convert_mvs_to_polars_lazy_as_series(obj, store=None):
+        return convert_pandas_to_polars(obj, lazy=True)
+
+    convert_dict[("pd.DataFrame", "pl.LazyFrame", "Series")] = (
+        convert_mvs_to_polars_lazy_as_series
+    )
+
+
 if _check_soft_dependencies("gluonts", severity="none"):
     from sktime.datatypes._adapter.gluonts import (
         convert_listDataset_to_pandas,
diff --git a/sktime/datatypes/_series/_examples.py b/sktime/datatypes/_series/_examples.py
index 0ad34481f..d1441e304 100644
--- a/sktime/datatypes/_series/_examples.py
+++ b/sktime/datatypes/_series/_examples.py
@@ -75,6 +75,13 @@ if _check_soft_dependencies("dask", severity="none"):
     example_dict[("dask_series", "Series", 0)] = df_dask
     example_dict_lossy[("dask_series", "Series", 0)] = False
 
+if _check_soft_dependencies("polars", severity="none"):
+    from polars import DataFrame
+
+    pl_df = DataFrame({"__index__0": [0, 1, 2, 3], "a": [1, 4, 0.5, -3]}, strict=False)
+
+    example_dict[("pl.DataFrame", "Series", 0)] = pl_df
+    example_dict_lossy[("pl.DataFrame", "Series", 0)] = False
 
 if _check_soft_dependencies("gluonts", severity="none"):
     from sktime.datatypes._adapter.gluonts import convert_pandas_to_listDataset
@@ -128,6 +135,19 @@ if _check_soft_dependencies("dask", severity="none"):
     example_dict[("dask_series", "Series", 1)] = df_dask
     example_dict_lossy[("dask_series", "Series", 1)] = False
 
+
+if _check_soft_dependencies("polars", severity="none"):
+    from polars import DataFrame
+
+    pl_df = DataFrame(
+        {"__index__0": [0, 1, 2, 3], "a": [1, 4, 0.5, -3], "b": [3, 7, 2, -3 / 7]},
+        strict=False,
+    )
+
+    example_dict[("pl.DataFrame", "Series", 1)] = pl_df
+    example_dict_lossy[("pl.DataFrame", "Series", 1)] = False
+
+
 if _check_soft_dependencies("gluonts", severity="none"):
     from sktime.datatypes._adapter.gluonts import convert_pandas_to_listDataset
 
@@ -181,6 +201,19 @@ if _check_soft_dependencies("dask", severity="none"):
     example_dict[("dask_series", "Series", 2)] = df_dask
     example_dict_lossy[("dask_series", "Series", 2)] = False
 
+
+if _check_soft_dependencies("polars", severity="none"):
+    from polars import DataFrame
+
+    pl_df = DataFrame(
+        {"__index__0": [0, 1, 2, 3], "a": [1, 4, 0.5, 3], "b": [3, 7, 2, 3 / 7]},
+        strict=False,
+    )
+
+    example_dict[("pl.DataFrame", "Series", 2)] = pl_df
+    example_dict_lossy[("pl.DataFrame", "Series", 2)] = False
+
+
 if _check_soft_dependencies("gluonts", severity="none"):
     from sktime.datatypes._adapter.gluonts import convert_pandas_to_listDataset
 
@@ -228,6 +261,16 @@ if _check_soft_dependencies("xarray", severity="none"):
     example_dict[("xr.DataArray", "Series", 3)] = da
     example_dict_lossy[("xr.DataArray", "Series", 3)] = False
 
+
+if _check_soft_dependencies("polars", severity="none"):
+    from polars import DataFrame
+
+    pl_df = DataFrame({"__index__0": [0, 1, 2, 3], "a": [1, 4, 0.5, 3]}, strict=False)
+
+    example_dict[("pl.DataFrame", "Series", 3)] = pl_df
+    example_dict_lossy[("pl.DataFrame", "Series", 3)] = False
+
+
 if _check_soft_dependencies("gluonts", severity="none"):
     from sktime.datatypes._adapter.gluonts import convert_pandas_to_listDataset
 
diff --git a/sktime/datatypes/_series/_registry.py b/sktime/datatypes/_series/_registry.py
index 0d5a0ccf1..a0d888410 100644
--- a/sktime/datatypes/_series/_registry.py
+++ b/sktime/datatypes/_series/_registry.py
@@ -32,6 +32,11 @@ MTYPE_REGISTER_SERIES = [
         "Series",
         "xdas representation of a uni- or multivariate series",
     ),
+    (
+        "pl.DataFrame",
+        "Series",
+        "pl.DataFrame representation of a uni- or multivariate series",
+    ),
     (
         "gluonts_ListDataset_series",
         "Series",
@@ -42,6 +47,7 @@ MTYPE_REGISTER_SERIES = [
 MTYPE_SOFT_DEPS_SERIES = {
     "xr.DataArray": "xarray",
     "dask_series": "dask",
+    "pl.DataFrame": "polars",
     "gluonts_ListDataset_series": "gluonts",
 }
 
diff --git a/sktime/forecasting/compose/_pipeline.py b/sktime/forecasting/compose/_pipeline.py
index 378618e55..ac0d30817 100644
--- a/sktime/forecasting/compose/_pipeline.py
+++ b/sktime/forecasting/compose/_pipeline.py
@@ -13,7 +13,6 @@ from sktime.forecasting.base._base import BaseForecaster
 from sktime.forecasting.base._delegate import _DelegatedForecaster
 from sktime.forecasting.base._fh import ForecastingHorizon
 from sktime.registry import scitype
-from sktime.utils.dependencies import _check_soft_dependencies
 from sktime.utils.validation.series import check_series
 from sktime.utils.warnings import warn
 
@@ -1787,8 +1786,10 @@ class ForecastX(BaseForecaster):
             instance.
             ``create_test_instance`` uses the first (or only) dictionary in ``params``
         """
+        from sktime.forecasting.arima import ARIMA
         from sktime.forecasting.compose import YfromX
         from sktime.forecasting.naive import NaiveForecaster
+        from sktime.utils.dependencies import _check_soft_dependencies
 
         fs, _ = YfromX.create_test_instances_and_names()
         fx = fs[0]
@@ -1797,9 +1798,8 @@ class ForecastX(BaseForecaster):
         params1 = {"forecaster_X": fx, "forecaster_y": fy}
 
         # example with probabilistic capability
-        if _check_soft_dependencies("pmdarima", severity="none"):
-            from sktime.forecasting.arima import ARIMA
-
+        # todo 0.33.0: check if numpy<2 is still needed
+        if _check_soft_dependencies(["pmdarima", "numpy<2"], severity="none"):
             fy_proba = ARIMA()
         else:
             fy_proba = NaiveForecaster()
diff --git a/sktime/forecasting/var_reduce.py b/sktime/forecasting/var_reduce.py
new file mode 100644
index 000000000..3aad8a2f5
--- /dev/null
+++ b/sktime/forecasting/var_reduce.py
@@ -0,0 +1,379 @@
+# copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
+"""Implements VARReduce, a VAR-like model combining tabularization with regression."""
+
+__author__ = ["meraldoantonio"]
+
+import numpy as np
+import pandas as pd
+
+from sktime.forecasting.base import BaseForecaster, ForecastingHorizon
+
+
+class VARReduce(BaseForecaster):
+    """
+    A flexible VAR-like forecaster that combines tabularization with regression.
+
+    The input data `Y_in` is a multivariate time series data containing `n` time series.
+    An example with `n` = 2:
+
+    +---------+-----+-----+
+    |  index  | ts1 | ts2 |
+    +---------+-----+-----+
+    |    1    | 11  | 6   |
+    |    2    | 12  | 7   |
+    |    3    | 13  | 8   |
+    |    4    | 14  | 9   |
+    |    5    | 15  | 10  |
+    +---------+-----+-----+
+
+    Fitting proceeds in two steps:
+    1. Tabularization:
+        For each time step and each time series within `Y_in`,
+        lagged values `X` are generated. The number of lagged values
+        are determined by the `lags` parameters.
+
+        Below is the `X` for the sample `Y_in` with `lags` = 2.
+        Note the absence of the earliest 2 timesteps as no
+        corresponding lag value is available.
+
+        +---------+----------+----------+----------+----------+
+        |  index  | ts1_lag1 | ts2_lag1 | ts1_lag2 | ts2_lag2 |
+        +---------+----------+----------+----------+----------+
+        |    3    | 12       | 7        | 11       | 6        |
+        |    4    | 13       | 8        | 12       | 7        |
+        |    5    | 14       | 9        | 13       | 8        |
+        +---------+----------+----------+----------+----------+
+
+    2. Regression:
+        The chosen regressor is fitted with `Y_in` as a target
+        and `X` as predictors. Care is taken to first remove the
+        first `lags` data points in `Y_in` as they do not have
+        corresponding indices in `X` (i.e. the first two data
+        points in the above example).
+
+    For forecasting, the last `lags` observations in `Y_in` are
+    reframed as lagged predictors `X_forecast` and passed to the
+    trained regressor to obtain the forecasts.
+    `X_forecast` is shown below.
+
+    +---------+----------+----------+----------+----------+
+    |  index  | ts1_lag1 | ts2_lag1 | ts1_lag2 | ts2_lag2 |
+    +---------+----------+----------+----------+----------+
+    |    6    | 15       | 10       | 14       | 9        |
+    +---------+----------+----------+----------+----------+
+
+    By default, `LinearRegression` is used, yielding results equivalent to
+    a traditional VAR model. Alternatively, any scikit-learn compatible regressor can
+    be used to introduce regularization and/or non-linearity.
+
+    For example:
+        - VARReduce(regressor = Ridge()) is equivalent to VAR with L2 regularization;
+        - VARReduce(regressor = Lasso()) is equivalent to VAR with L1 regularization.
+
+    These two models can be used to incorporate regularization and prevent overfitting
+    when the input data contain a large number of individual time series relative to
+    data points.
+
+    Parameters
+    ----------
+    lags : int, optional, default=1
+        The number of lagged values to include in the model.
+    regressor : object, optional (default=LinearRegression())
+        The regressor to use for fitting the model. Must be scikit-learn-compatible.
+
+    Attributes
+    ----------
+    coefficients_ : np.ndarray, shape (lags, num_series, num_series)
+        The estimated coefficients of the model;
+        only available if the regressor has `coef_` attribute
+    intercept_ : np.ndarray, shape (num_series,)
+        The intercept for each time series;
+        only available if the regressor has `coef_` attribute
+    num_series : int
+        The number of time series being modeled.
+    var_names : list of str
+        The names of the time series being modeled
+
+    Examples
+    --------
+    >>> from sktime.forecasting.var_reduce import VARReduce
+    >>> from sklearn.linear_model import Lasso
+    >>> from sktime.datasets import load_longley
+    >>> _, y = load_longley()
+    >>> forecaster = VARReduce(regressor=Lasso())  # doctest: +SKIP
+    >>> forecaster.fit(y)  # doctest: +SKIP
+    VARReduce(...)
+    >>> y_pred = forecaster.predict(fh=[1,2,3])  # doctest: +SKIP
+    """
+
+    _tags = {
+        "scitype:y": "both",
+        "authors": ["meraldoantonio"],
+        "y_inner_mtype": "pd.DataFrame",
+        "X_inner_mtype": "pd.DataFrame",
+        "ignores-exogeneous-X": True,
+        "requires-fh-in-fit": False,
+    }
+
+    def __init__(self, lags=1, regressor=None):
+        from sklearn.base import clone
+        from sklearn.linear_model import LinearRegression
+
+        self.regressor = regressor  # not used/modified
+        self.lags = lags
+        if regressor is None:
+            self.regressor_ = LinearRegression()
+        else:
+            self.regressor_ = clone(regressor)
+
+        assert hasattr(self.regressor_, "fit"), "Regressor must have 'fit'"
+        assert hasattr(self.regressor_, "predict"), "Regressor must have 'predict'"
+
+        self.coefficients_ = None
+        self.intercept_ = None
+        self.num_series = None
+        self.var_names = None
+        super().__init__()
+
+    def _prepare_for_fit(self, data, return_as_ndarray=True):
+        """
+        Prepare the data for VAR fitting though tabularization.
+
+        This function transforms the provided training data into
+        a tabular format suitable for regression. Specifically,
+        the predictors X consist of lagged values of the multivariate time series, while
+        the target variables Y are the current values of the multivariate time series.
+
+        Parameters
+        ----------
+        data : pd.DataFrame
+            The input multivariate time series data to be transformed
+            with shape (num_samples, num_series)
+        return_as_ndarray : bool, optional (default=True)
+            If True, returns the transformed data as numpy arrays.
+            If False, returns the transformed data as pandas DataFrames.
+
+        Returns
+        -------
+        X : np.ndarray or pd.DataFrame
+            The lagged values as predictors;
+            with shape (num_samples, num_series * lags)
+        Y : np.ndarray or pd.DataFrame
+            The original, unlagged data;
+            with shape (num_samples, num_series).
+        """
+        if isinstance(data, np.ndarray):
+            data = pd.DataFrame(data, columns=self.var_names)
+
+        # Generate the lagged values and rename the columns
+        df = pd.concat([data.shift(i) for i in range(self.lags + 1)], axis=1)
+        df.columns = [
+            f"{col}_lag{i}" for i in range(self.lags + 1) for col in data.columns
+        ]
+        df = df.dropna()
+
+        # Separate df into X and Y
+        Y = df[[f"{col}_lag0" for col in data.columns]]
+        Y.columns = data.columns
+        X = df.drop(columns=[f"{col}_lag0" for col in data.columns])
+
+        if return_as_ndarray:
+            return X.values, Y.values
+        else:
+            return X, Y
+
+    def _prepare_for_predict(self, data, return_as_ndarray=True):
+        """
+        Prepare the data for VAR prediction through tabularization.
+
+        This function extracts the last `lags` rows of the provided data
+        and linearizes them into a single row suitable for input into the regressor.
+        The number of lags is inferred from the `self.lags` attribute.
+
+        Parameters
+        ----------
+        data : pd.DataFrame
+            The input multivariate time series data to be transformed
+            with shape (num_samples, num_series)
+            Usually is training data
+
+        return_as_ndarray : bool, optional (default=True)
+            If True, returns the transformed data as a NumPy array.
+            If False, returns the transformed data as a pandas DataFrame
+
+        Returns
+        -------
+        np.ndarray or pd.DataFrame
+            A single row of the transformed data, with shape (1, num_series * lags).
+            If returned as a DataFrame, the columns are named according to the original
+            column names with the lag number appended (e.g., 'A_lag1', 'B_lag2').
+        """
+        # Take the last lags rows of the input data
+        if isinstance(data, np.ndarray):
+            data = pd.DataFrame(data, columns=self.var_names)
+        lags = self.lags
+        lagged_data = data.tail(lags)
+
+        # Invert the data
+        inverted_data = lagged_data[::-1]
+
+        # Linearize it
+        linearized_data = inverted_data.values.flatten().reshape(1, -1)
+
+        # Return as np.array or pd.DataFrame
+        if return_as_ndarray:
+            return np.array(linearized_data)
+        else:
+            # Generate column names
+            columns = []
+            for lag in range(1, lags + 1):
+                columns += [f"{col}_lag{lag}" for col in data.columns]
+            return pd.DataFrame(linearized_data, columns=columns)
+
+    def _fit(self, y, X=None, fh=None):
+        """Fit forecaster to training data.
+
+        Parameters
+        ----------
+        y : pd.DataFrame
+            Guaranteed to have a single column if scitype:y=="univariate".
+        fh : ForecastingHorizon, optional (default=None)
+            The forecasting horizon with the steps ahead to to predict.
+        X : pd.DataFrame, optional (default=None)
+            Exogenous time series to fit to; will be ignored
+
+        Returns
+        -------
+        self : reference to self
+        """
+        from sklearn.multioutput import MultiOutputRegressor
+
+        self.var_names = y.columns
+        self.num_series = y.shape[1]
+
+        X, Y = self._prepare_for_fit(y, return_as_ndarray=False)
+
+        native_multioutput_support = self.regressor_._get_tags().get(
+            "multioutput", False
+        )
+        if not native_multioutput_support:
+            self.regressor_ = MultiOutputRegressor(self.regressor_)
+        self.regressor_.fit(X, Y)
+
+        # if the model has `.intercept_` and `.coef_` attributes, extract them
+        if hasattr(self.regressor_, "intercept_") and hasattr(self.regressor_, "coef_"):
+            self.coefficients_ = self.regressor_.coef_.reshape(
+                (self.num_series, self.lags, self.num_series)
+            )
+            self.coefficients_ = np.transpose(self.coefficients_, (1, 0, 2))
+            self.intercept_ = self.regressor_.intercept_
+        else:
+            self.coefficients_ = None
+            self.intercept_ = None
+        return self
+
+    def _predict(self, fh, X=None):
+        """Forecast time series at future horizon.
+
+        Parameters
+        ----------
+        fh : ForecastingHorizon
+            The forecasters horizon with the steps ahead to to predict.
+            Default is one-step ahead forecast,
+            i.e. np.array([1])
+        X : pd.DataFrame, optional (default=None)
+            Exogenous variables are ignored.
+
+        Returns
+        -------
+        y_pred : pd.DataFrame
+            Series of predicted values
+        """
+        y_pred_outsample = pd.DataFrame()
+        y_pred_insample = pd.DataFrame()
+        fh_int = fh.to_relative(self.cutoff)
+
+        # ---- insample forecasts  -----
+        if fh_int.min() <= 0:
+            # Reproduce the original X we used for fitting
+            X, _ = self._prepare_for_fit(self._y, return_as_ndarray=False)
+
+            self._y_pred_insample = pd.DataFrame(
+                self.regressor_.predict(X),
+                index=X.index,
+                columns=self.var_names,
+            )
+
+            # Create a new `fh` object with only in-sample values
+            # and use it to filter previously-created in-sample predictions
+            fh_insample = fh_int.to_in_sample(cutoff=self.cutoff)
+            insample_index = fh_insample.to_absolute_index(cutoff=self.cutoff)
+            y_pred_insample = self._y_pred_insample.loc[insample_index]
+
+        # ---- outsample forecasts ----
+        if fh_int.max() > 0:
+            # Get the last available values for prediction
+            y_last = self._y.iloc[-self.lags :]
+
+            # Initialize a list to store out-of-sample predictions
+            y_pred_outsample = []
+
+            # Generate as many future forecast steps as the maximum number in fh_int
+            for _ in range(0, fh_int[-1]):
+                # Prepare X
+                X_last = self._prepare_for_predict(y_last, return_as_ndarray=False)
+                y_pred_step = self.regressor_.predict(X_last)
+                y_pred_step = pd.DataFrame(y_pred_step, columns=self.var_names)
+                y_pred_outsample.append(y_pred_step)
+
+                # Append the new predictions (y_pred_step) at the end of y_last
+                y_last = pd.concat([y_last, y_pred_step], axis=0)
+
+            # Concatenate the list of Series into the final DataFrame
+            y_pred_outsample = pd.concat(y_pred_outsample)
+
+            # Create an index for all generated forecasts
+            fh_outsample_all = ForecastingHorizon(
+                range(1, fh_int[-1] + 1), is_relative=True, freq=fh.freq
+            )
+            fh_outsample_all = fh_outsample_all.to_absolute(cutoff=self.cutoff)
+            y_pred_outsample.index = fh_outsample_all.to_absolute_index()
+
+            # Create a new `fh` object with only the required outsample indices
+            fh_outsample = fh_int.to_out_of_sample(cutoff=self.cutoff)
+            fh_outsample = fh_outsample.to_absolute(cutoff=self.cutoff)
+
+            # Use it to filter the previously generated outsample forecasts
+            outsample_index = fh_outsample.to_absolute_index()
+            y_pred_outsample = y_pred_outsample.loc[outsample_index]
+
+        # Concatenate filtered y_pred_insample and y_pred_outsample
+        y_pred = pd.concat([y_pred_insample, y_pred_outsample])
+        return y_pred
+
+    @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`
+        """
+        from sklearn.ensemble import RandomForestRegressor
+        from sklearn.svm import SVR
+
+        params1 = {"lags": 2, "regressor": RandomForestRegressor()}
+        params2 = {"lags": 3, "regressor": SVR(kernel="rbf")}
+
+        return [params1, params2]


Test patch:
diff --git a/sktime/datasets/tests/test_hierarchical_demo.py b/sktime/datasets/tests/test_hierarchical_demo.py
new file mode 100644
index 000000000..8914cca49
--- /dev/null
+++ b/sktime/datasets/tests/test_hierarchical_demo.py
@@ -0,0 +1,18 @@
+"""Test functions for load_hierarchical_sales_toydata."""
+
+__author__ = ["fkiraly"]
+
+from sktime.datatypes import check_raise
+
+
+def test_load_hierarchical_sales_toydata():
+    """Test that load_hierarchical_sales_toydata runs and returns expected format."""
+    from sktime.datasets import load_hierarchical_sales_toydata
+
+    df = load_hierarchical_sales_toydata()
+
+    check_raise(
+        df,
+        "pd_multiindex_hier",
+        var_name="return of load_hierarchical_sales_toydata",
+    )
diff --git a/sktime/datasets/tests/test_loaders.py b/sktime/datasets/tests/test_loaders.py
index af82ef1f5..4211fc90f 100644
--- a/sktime/datasets/tests/test_loaders.py
+++ b/sktime/datasets/tests/test_loaders.py
@@ -1,4 +1,4 @@
-"""Test functions for loose loaers."""
+"""Test functions for loose loaders."""
 
 __author__ = ["fkiraly"]
 
diff --git a/sktime/datatypes/_adapter/tests/test_polars.py b/sktime/datatypes/_adapter/tests/test_polars.py
index 2ce4d3c7d..5c5bcb323 100644
--- a/sktime/datatypes/_adapter/tests/test_polars.py
+++ b/sktime/datatypes/_adapter/tests/test_polars.py
@@ -3,6 +3,7 @@
 
 from copy import deepcopy
 
+import numpy as np
 import pandas as pd
 import pytest
 
@@ -32,7 +33,7 @@ PANDAS_FIXTURES = [pd_fixture_simple, pd_fixture_multiindex]
     reason="skip test if required soft dependency for polars not available",
 )
 def test_convert_polars_to_pandas_multiindex():
-    """Tests conversion from polars to pandas is correct for MultiIndex."""
+    """Tests conversion from polars to pandas."""
     from polars import from_pandas
 
     pd_fixture = pd.DataFrame(
@@ -57,6 +58,7 @@ def test_convert_polars_to_pandas_multiindex():
     reason="skip test if required soft dependency for polars not available",
 )
 def test_convert_pd_multiindex_to_polars():
+    """Tests conversion from pandas to polars is correct for MultiIndex."""
     pd_fixture = pd_fixture_multiindex.copy()
     pd_fixture_orig = deepcopy(pd_fixture)
     result = convert_pandas_to_polars(pd_fixture)
@@ -96,3 +98,60 @@ def test_convert_pd_polars_inverse_lazy(pd_fixture):
     pd_result = convert_polars_to_pandas(polars_result)
 
     assert pd_result.equals(pd_fixture)
+
+
+@pytest.mark.skipif(
+    not _check_soft_dependencies("polars", severity="none"),
+    reason="skip test if required soft dependency for polars not available",
+)
+def test_pandas_to_polars_with_pandas_period_index_and_inverse():
+    """Tests conversion from pandas to polars with pandas PeriodIndex."""
+    df = pd.DataFrame(
+        {"a": [1, 2, 3], "b": [1, 2, 3], "foo": [2, 3, 4], "bar": [3, 4, 5]}
+    )
+    df = df.set_index(pd.period_range("2021-01-01", periods=3))
+
+    pl_result = convert_pandas_to_polars(df)
+    pd_result = convert_polars_to_pandas(pl_result, infer_freq=True)
+
+    assert pl_result.columns == ["__index__0", "a", "b", "foo", "bar"]
+    assert list(pd_result.columns) == list(df.columns)
+    assert pd_result.index.freq == df.index.freq
+
+
+@pytest.mark.skipif(
+    not _check_soft_dependencies("polars", severity="none"),
+    reason="skip test if required soft dependency for polars not available",
+)
+def test_pd_multiindex_to_polars_with_pandas_period_index_and_inverse():
+    """Tests conversion from pandas to polars with pandas PeriodIndex for Multiindex."""
+    data = np.random.randn(6, 3)
+
+    arrays = [
+        ["A", "A", "B", "B", "C", "C"],
+        [1, 2, 1, 2, 1, 2],
+        pd.period_range("2023-01", periods=6, freq="M"),
+    ]
+
+    index = pd.MultiIndex.from_arrays(arrays, names=["1", "2", "period"])
+    pd_multiindex_period_index = pd.DataFrame(
+        data, index=index, columns=["foo", "bar", "a"]
+    )
+
+    pl_result = convert_pandas_to_polars(pd_multiindex_period_index)
+    pd_result = convert_polars_to_pandas(pl_result, infer_freq=True)
+
+    assert pl_result.columns == [
+        "__index__1",
+        "__index__2",
+        "__index__period",
+        "foo",
+        "bar",
+        "a",
+    ]
+    assert list(pd_result.columns) == list(pd_multiindex_period_index.columns)
+    assert pd_result.index.names == pd_multiindex_period_index.index.names
+    assert (
+        pd_result.index.levels[-1].freq
+        == pd_multiindex_period_index.index.levels[-1].freq
+    )
diff --git a/sktime/forecasting/base/tests/test_fh.py b/sktime/forecasting/base/tests/test_fh.py
index 9b5fd56f0..3e99663b5 100644
--- a/sktime/forecasting/base/tests/test_fh.py
+++ b/sktime/forecasting/base/tests/test_fh.py
@@ -733,9 +733,11 @@ def test_exponential_smoothing_case_with_naive():
 
 
 # TODO: Replace this long running test with fast unit test
+# todo 0.33.0: check whether numpy 2 bound is still necessary
 @pytest.mark.skipif(
     not run_test_module_changed(["sktime.forecasting.base", "sktime.datatypes"])
-    or not _check_estimator_deps(AutoARIMA, severity="none"),
+    or not _check_estimator_deps(AutoARIMA, severity="none")
+    or _check_soft_dependencies("numpy>=2.0", severity="none"),
     reason="run only if base module has changed or datatypes module has changed",
 )
 def test_auto_arima():
diff --git a/sktime/forecasting/tests/test_varreduce.py b/sktime/forecasting/tests/test_varreduce.py
new file mode 100644
index 000000000..c49dce566
--- /dev/null
+++ b/sktime/forecasting/tests/test_varreduce.py
@@ -0,0 +1,24 @@
+"""Tests the VARReduce model against statsmodels VAR."""
+
+import numpy as np
+import pytest
+
+from sktime.datasets import load_longley
+from sktime.forecasting.var import VAR
+from sktime.forecasting.var_reduce import VARReduce
+from sktime.tests.test_switch import run_test_for_class
+
+__author__ = ["meraldoantonio"]
+
+
+@pytest.mark.skipif(
+    not run_test_for_class([VAR, VARReduce]),
+    reason="run test only if softdeps are present and incrementally (if requested)",
+)
+def test_VAR_against_statsmodels():
+    """Compares default VARReduce to statmodels VAR."""
+    y = load_longley()[1]
+
+    y_var = VAR().fit_predict(y, fh=[1, 2, 3])
+    y_varr = VARReduce().fit_predict(y, fh=[1, 2, 3])
+    np.testing.assert_allclose(y_var.values, y_varr.values)
diff --git a/sktime/tests/test_softdeps.py b/sktime/tests/test_softdeps.py
index ccf7f1c09..6f0c7bb7a 100644
--- a/sktime/tests/test_softdeps.py
+++ b/sktime/tests/test_softdeps.py
@@ -90,8 +90,11 @@ def is_soft_dep_missing_message(msg):
     # special message for deep learning dependencies
     error_msg_dl = "required for deep learning"
     cond3 = error_msg_dl in msg
+    # message if environment marker not satisfied
+    error_msg_marker = "packaging marker"
+    cond4 = error_msg_marker in msg
 
-    return cond1 or cond2 or cond3
+    return cond1 or cond2 or cond3 or cond4
 
 
 @pytest.mark.parametrize("module", modules)