Repository: alan-turing-institute/sktime
Pull request #2763 (2022-06-08 05:35:03)

[ENH] statsmodels VARMAX adapter
#### Reference Issues/PRs
[Add VARMAX statespace model in sktime #2499](https://github.com/alan-turing-institute/sktime/issues/2499)

#### What does this implement/fix? Explain your changes.
Interfacing VARMAX from statsmodels, with corresponding unit tests

#### Does your contribution introduce a new dependency? If yes, which one?
No

#### What should a reviewer concentrate their feedback on?

#### Any other comments?
Questions:
- Is there a list of estimators I need to add VARMAX to?
- In var.py, there's code in predict() that I don't understand. I initially copied across this (with out-sample predictions using forecast and in-sample predictions using the diff of residuals from the fitted object). I'm not sure why? Instead, for this, I used the predict() method from statsmodels (returns In-sample predictions / Out-of-sample forecasts). Checking that I'm not missing anything here!
https://www.statsmodels.org/devel/generated/statsmodels.tsa.statespace.varmax.VARMAXResults.predict.html
- get_prediction() in statsmodels also gives confidence intervals, do we want to include this?
- In the example (provided in the doc string), it fails the test for stationarity and hence we get: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.
should we normalise this dataset in the example(s)? (same data in the VAR doc string example)

##### For new estimators
- [ ] I've added the estimator to the online documentation.
- [ ] I've updated the existing example notebooks or provided a new one to showcase how my estimator works.

@fkiraly @lmmentel @Lovkush-A 

KatieBuc: Some set-up issues to hopefully resolve with mentors next session:
- pre-commit extension not working with my VScode (need to check linting)
Extension installed.
Settings > Pre-commit-helper: Run On Save. Select, which hooks to run after saving a file.
Set to all hooks
terminal > "The terminal process failed to launch: Path to shell executable "C:\Users\n10907700\repos\sktime\pre-commit" does not exist."

- pytest discovery error in search bar of VScode (settings changed, restated. For now, tested via notebook)
![image](https://user-images.githubusercontent.com/34525024/172542529-232ccbbe-05d2-492c-b74b-3081f25db842.png)

@fkiraly @lmmentel @Lovkush-A 
Lovkush-A: > Some set-up issues to hopefully resolve with mentors next session:

unfortunately i cannot help (immediately anyway) with this. i have not used precommit with vscode, and i do not remember what i did to get pytest to work in vscode... :/

can you get pre-commit to work in terminal? i.e. Can you manage to get pre-commit to do stuff when you run `git commit ...` in terminal?
fkiraly: you can fix linting issues manually as well, either by manually running the formatting packages, or by typing out the suggested changes (which may be the most convenient option if there are only few)
fkiraly: @KatieBuc, your branch has a lot of linting (= code formatting) issues.

The best way to deal with these is to learn how to run the code formatters locally, and/or how to install precommit. See the developer docs here: https://www.sktime.org/en/stable/developer_guide/continuous_integration.html

Another way would be to go on "details" next to the failed checks on the GitHub web interface and see the logs.
lmmentel: I'm not sure if VSCode plugin installs `pre-commit` itself but it might be worth checking if you have `pre-commit` installed and if your plugin configuration is pointing the right location for `pre-commit` i.e. your virtual env.
fkiraly: > I'm not sure if VSCode plugin installs `pre-commit` itself 

I think it doesn´t, one has to do it manually
KatieBuc: Okay, thanks everyone, pre-commit can run locally

Lovkush-A: @KatieBuc Looking good! Making progress on the automated tests - got passed the code quality. 💪 

It looks like the parameter setting `{"order": (0, 0)}` is invalid (at least one of p and q should be non-zero, is what error message from statsmodel says) and that is causing a bunch of the automated tests to fail.
Lovkush-A: > * Is there a list of estimators I need to add VARMAX to?

I am not sure. @fkiraly?

> * In var.py, there's code in predict() that I don't understand. I initially copied across this (with out-sample predictions using forecast and in-sample predictions using the diff of residuals from the fitted object). I'm not sure why? Instead, for this, I used the predict() method from statsmodels (returns In-sample predictions / Out-of-sample forecasts). Checking that I'm not missing anything here!

- For out-sample, seems like using predict() is the only option.
- For in-sample, seems like there are two options: one is to use residuals and the other to use predict().  I do not know VARMAX or statsmodels well enough to say what the difference is between these too (or even if there is a difference).
- Question just to check: do you know what 'in sample' and 'out of sample' means in this context?
- Question: do you know how the residuals and the predictions are calculated, and whether they produce different answers?


> * get_prediction() in statsmodels also gives confidence intervals, do we want to include this?

Yes, I think so. If you look in [forecasting extension template](https://github.com/alan-turing-institute/sktime/blob/main/extension_templates/forecasting.py), there are method(s) one should implement for confidence intervals. Sounds like `_predict_interval` is the method to implement.

> hence we get: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.
> should we normalise this dataset in the example(s)?

I do not know VARMAX, so cannot say whether normalising dataset would avoid the error. Based on the error message, an alternative solution could be to explicitly provide `start_params`. 


fkiraly: >>Is there a list of estimators I need to add VARMAX to?

> I am not sure. @fkiraly?

In general not, the registry retrieves estimators automatically from the code base.
Maybe in `build_tools.azure.check_soft_dependencies`, but we're not using Azure anymore, so that list does not need being added to, right? @lmmentel ?

fkiraly: (very superficial look at CI failure traceback: looks like a `pandas` index related syntax error in `_predict`)
KatieBuc: ignore this last commit :)





fkiraly: ok, ignoring this last commit and definitely not looking into it, or observing, for instance, that it says linting but does not look like linting but addition of actual conversion/indexing logic
KatieBuc: > ok, ignoring this last commit and definitely not looking into it, or observing, for instance, that it says linting but does not look like linting but addition of actual conversion/indexing logic

I'm so glad you didn't notice that @fkiraly 

Okay - debugging - that was a wild ride. I got deep into the codebase but made it out alive. As it turns out, statsmodels have an indexing quirk.

This gives a quick example of it:
`
import pandas as pd
from sktime.utils._testing.series import _make_series

y_range = _make_series(
    n_columns=2, index_type='range', n_timepoints=5
)

print(y_range)

y_int = _make_series(
    n_columns=2, index_type='int', n_timepoints=5
)

print(y_int)

from statsmodels.tsa.api import VARMAX as _VARMAX

stats = _VARMAX(y_range)
stats_fit = stats.fit()
y_pred_stats = stats_fit.predict(start=0, end=10)
print(y_pred_stats)

stats = _VARMAX(y_int)
stats_fit = stats.fit()
y_pred_stats = stats_fit.predict(start=0, end=4)
print(y_pred_stats)

# ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.
stats = _VARMAX(y_int)
stats_fit = stats.fit()
y_pred_stats = stats_fit.predict(start=0, end=20)
print(y_pred_stats)

# ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.
stats = _VARMAX(y_int)
stats_fit = stats.fit()
y_pred_stats = stats_fit.predict(start=9, end=20)
print(y_pred_stats)
`
Basically, when `index_type='int'` and for the out-of-sample forecasting, statsmodels returns a new index reset at zero! Hence the new if statement.

I'm wondering if there's a better way (involving 'capturing' the ValueWarning somehow). This is passing all unit tests though.
KatieBuc: @Lovkush-A :)

One thing it's still complaining about is 
![image](https://user-images.githubusercontent.com/34525024/174510715-14d17e53-66e9-4e36-92cf-6c01e339ceb1.png)

But I've defini9tely added a blank line 
[weird docstring linting](https://github.com/alan-turing-institute/sktime/pull/2763/commits/c77fbe479d017446324ac1d22fad3ada638bf4d7)
[c77fbe4](https://github.com/alan-turing-institute/sktime/pull/2763/commits/c77fbe479d017446324ac1d22fad3ada638bf4d7)

![image](https://user-images.githubusercontent.com/34525024/174510787-ed617a11-5a29-466f-9329-492ad7914c34.png)

I'm not sure what this error wants from me!
fkiraly: as discussed in standup:
the message is not about a blank new line, but a blank character.

e.g., the first line should be `>>> from etc` instead of `>>>from etc`, and similar for the other lines.
fkiraly: The failure is in the first line, not the second one as far as I can tell.
`FAILED forecasting/varmax.py::sktime.forecasting.varmax.VARMAX `

means that the example, from the docstring, when run, is failing.

The easiest way to check what is happening is to copy all the lines starting with `>>> ` into a vs code python interpreter window and run them (it will run it despite the `>>> `). If they break, or the printout does not exactly equal example printouts in the docstring example, the doctest will fail.
KatieBuc: > looks good, could you kindly also add it to docs? You can check `VAR` e.g. of how to do that. Otherwise I can support :)

ohhhhh @fkiraly just showed me, and that answers my very first question! Okay, pushed
fkiraly: >In var.py, there's code in predict() that I don't understand. I initially copied across this (with out-sample predictions using forecast and in-sample predictions using the diff of residuals from the fitted object). I'm not sure why? Instead, for this, I used the predict() method from statsmodels (returns In-sample predictions / Out-of-sample forecasts). Checking that I'm not missing anything here!

It would look like `statsmodels` `varmax` can indeed deal with both in-sample and out-of-sample predictions. So, not doing anything fancy/custom but juts calling the `statsmodels` functionality that already exists seems the right thing to do.

> `get_prediction()` in `statsmodels` also gives confidence intervals, do we want to include this?

Are these confidence intervals (on the predictive mean), or predictive intervals? Which are not the same concept. We would wanti t only in the latter case, since that is what `predict_interval` is meant to return.

> In the example (provided in the doc string), it fails the test for stationarity and hence we get: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.
should we normalise this dataset in the example(s)? (same data in the VAR doc string example)

I suppose? Probably some non-stationarity test internal to `varmax` raises the alarm here.
Normalisation might not help, differencing typically does.

Patch:
diff --git a/docs/source/api_reference/forecasting.rst b/docs/source/api_reference/forecasting.rst
index b7a7a153e..f04d76e56 100644
--- a/docs/source/api_reference/forecasting.rst
+++ b/docs/source/api_reference/forecasting.rst
@@ -201,6 +201,14 @@ All "ARIMA" models below include SARIMAX capability.
 
     VAR
 
+.. currentmodule:: sktime.forecasting.varmax
+
+.. autosummary::
+    :toctree: auto_generated/
+    :template: class.rst
+
+    VARMAX
+
 Structural time series models
 -----------------------------
 
diff --git a/examples/01_forecasting.ipynb b/examples/01_forecasting.ipynb
index 51e6bc5b0..062192828 100644
--- a/examples/01_forecasting.ipynb
+++ b/examples/01_forecasting.ipynb
@@ -1128,7 +1128,7 @@
     "#  e.g., the non-symmetric verion\n",
     "# it also allows for inspection of metric properties\n",
     "#  e.g., are higher values better (answer: no)?\n",
-    "mape.greater_is_better"
+    "mape.get_tag(\"lower_is_better\")"
    ]
   },
   {
diff --git a/sktime/base/_base.py b/sktime/base/_base.py
index 9dea5fa86..a6abfe77f 100644
--- a/sktime/base/_base.py
+++ b/sktime/base/_base.py
@@ -630,14 +630,6 @@ class TagAliaserMixin:
     ensure toremove this class as a parent of BaseObject or BaseEstimator.
     """
 
-    # dictionary of aliases
-    # key = old tag; value = new tag, aliased by old tag
-    alias_dict = {"capability:early_prediction": ""}
-
-    # dictionary of removal version
-    # key = old tag; value = version in which tag will be removed, as string
-    deprecate_dict = {"capability:early_prediction": "0.13.0"}
-
     def __init__(self):
         super(TagAliaserMixin, self).__init__()
 
@@ -795,7 +787,7 @@ class TagAliaserMixin:
                 warnings.warn(msg, category=DeprecationWarning)
 
 
-class BaseEstimator(TagAliaserMixin, BaseObject):
+class BaseEstimator(BaseObject):
     """Base class for defining estimators in sktime.
 
     Extends sktime's BaseObject to include basic functionality for fittable estimators.
diff --git a/sktime/forecasting/base/_fh.py b/sktime/forecasting/base/_fh.py
index f871dfdd0..8d9bc6712 100644
--- a/sktime/forecasting/base/_fh.py
+++ b/sktime/forecasting/base/_fh.py
@@ -643,7 +643,7 @@ def _to_relative(fh: ForecastingHorizon, cutoff=None) -> ForecastingHorizon:
         # Out: Index([<0 * Hours>, <4 * Hours>, <8 * Hours>], dtype = 'object')
         # [v - periods[0] for v in periods]
         # Out: Index([<0 * Hours>, <2 * Hours>, <4 * Hours>], dtype='object')
-        # TODO: v0.12.1 OR 0.13.0: Check if this comment below can be removed,
+        # TODO: 0.14.0: Check if this comment below can be removed,
         # so check if pandas has released the fix to PyPI:
         # This bug was reported: https://github.com/pandas-dev/pandas/issues/45999
         # and fixed: https://github.com/pandas-dev/pandas/pull/46006
diff --git a/sktime/forecasting/model_selection/_tune.py b/sktime/forecasting/model_selection/_tune.py
index 412efa115..c5a42bb08 100644
--- a/sktime/forecasting/model_selection/_tune.py
+++ b/sktime/forecasting/model_selection/_tune.py
@@ -257,7 +257,7 @@ class BaseGridSearch(BaseForecaster):
 
         # Rank results, according to whether greater is better for the given scoring.
         results[f"rank_{scoring_name}"] = results.loc[:, f"mean_{scoring_name}"].rank(
-            ascending=not scoring.greater_is_better
+            ascending=scoring.get_tag("lower_is_better")
         )
 
         self.cv_results_ = results
@@ -274,7 +274,7 @@ class BaseGridSearch(BaseForecaster):
 
         # Sort values according to rank
         results = results.sort_values(
-            by=f"rank_{scoring_name}", ascending=not scoring.greater_is_better
+            by=f"rank_{scoring_name}", ascending=scoring.get_tag("lower_is_better")
         )
         # Select n best forecaster
         self.n_best_forecasters_ = []
diff --git a/sktime/forecasting/varmax.py b/sktime/forecasting/varmax.py
new file mode 100644
index 000000000..56b30b668
--- /dev/null
+++ b/sktime/forecasting/varmax.py
@@ -0,0 +1,405 @@
+# -*- coding: utf-8 -*-
+"""Vector Autoregressive Moving Average with eXogenous regressors model (VARMAX)."""
+__all__ = ["VARMAX"]
+__author__ = ["KatieBuc"]
+
+import warnings
+
+import pandas as pd
+from statsmodels.tsa.statespace.varmax import VARMAX as _VARMAX
+
+from sktime.forecasting.base.adapters import _StatsModelsAdapter
+
+
+class VARMAX(_StatsModelsAdapter):
+    r"""Wrapper for statsmodels VARMAX model.
+
+    Vector Autoregressive Moving Average with eXogenous regressors model (VARMAX)
+
+    Parameters
+    ----------
+    order : iterable
+        The (p,q) order of the model for the number of AR and MA parameters to
+        use.
+    trend : str{'n','c','t','ct'} or iterable, optional
+        Parameter controlling the deterministic trend polynomial :math:`A(t)`.
+        Can be specified as a string where 'c' indicates a constant (i.e. a
+        degree zero component of the trend polynomial), 't' indicates a
+        linear trend with time, and 'ct' is both. Can also be specified as an
+        iterable defining the non-zero polynomial exponents to include, in
+        increasing order. For example, `[1,1,0,1]` denotes
+        :math:`a + bt + ct^3`. Default is a constant trend component.
+    error_cov_type : {'diagonal', 'unstructured'}, optional
+        The structure of the covariance matrix of the error term, where
+        "unstructured" puts no restrictions on the matrix and "diagonal"
+        requires it to be a diagonal matrix (uncorrelated errors). Default is
+        "unstructured".
+    measurement_error : bool, optional
+        Whether or not to assume the endogenous observations `endog` were
+        measured with error. Default is False.
+    enforce_stationarity : bool, optional
+        Whether or not to transform the AR parameters to enforce stationarity
+        in the autoregressive component of the model. Default is True.
+    enforce_invertibility : bool, optional
+        Whether or not to transform the MA parameters to enforce invertibility
+        in the moving average component of the model. Default is True.
+    trend_offset : int, optional
+        The offset at which to start time trend values. Default is 1, so that
+        if `trend='t'` the trend is equal to 1, 2, ..., n_obs. Typically is only
+        set when the model created by extending a previous dataset.
+    start_params : array_like, optional
+        Initial guess of the solution for the loglikelihood maximization. If
+        None, the default is given by Model.start_params.
+    transformed : bool, optional
+        Whether or not start_params is already transformed. Default is True.
+    includes_fixed : bool, optional
+        If parameters were previously fixed with the fix_params method, this
+        argument describes whether or not start_params also includes the
+        fixed parameters, in addition to the free parameters. Default is False.
+    cov_type : str, optional
+        The `cov_type` keyword governs the method for calculating the
+        covariance matrix of parameter estimates. Can be one of:
+         - 'opg' for the outer product of gradient estimator
+         - 'oim' for the observed information matrix estimator, calculated
+            using the method of Harvey (1989)
+         - 'approx' for the observed information matrix estimator,
+            calculated using a numerical approximation of the Hessian matrix.
+         - 'robust' for an approximate (quasi-maximum likelihood) covariance
+            matrix that may be valid even in the presence of some
+            misspecifications. Intermediate calculations use the 'oim'
+            method.
+         - 'robust_approx' is the same as 'robust' except that the
+            intermediate calculations use the 'approx' method.
+         - 'none' for no covariance matrix calculation.
+        Default is 'opg' unless memory conservation is used to avoid computing the
+        loglikelihood values for each observation, in which case the default is
+        'approx'.
+    cov_kwds : dict or None, optional
+        A dictionary of arguments affecting covariance matrix computation.
+        opg, oim, approx, robust, robust_approx
+         - 'approx_complex_step' : bool, optional - If True, numerical
+            approximations are computed using complex-step methods. If False,
+            numerical approximations are computed using finite difference
+            methods. Default is True.
+         - 'approx_centered' : bool, optional - If True, numerical
+            approximations computed using finite difference methods use a
+            centered approximation. Default is False.
+    method : str, optional
+        The `method` determines which solver from `scipy.optimize`
+        is used, and it can be chosen from among the following strings:
+         - 'newton' for Newton-Raphson
+         - 'nm' for Nelder-Mead
+         - 'bfgs' for Broyden-Fletcher-Goldfarb-Shanno (BFGS)
+         - 'lbfgs' for limited-memory BFGS with optional box constraints
+         - 'powell' for modified Powell's method
+         - 'cg' for conjugate gradient
+         - 'ncg' for Newton-conjugate gradient
+         - 'basinhopping' for global basin-hopping solver
+        The explicit arguments in `fit` are passed to the solver,
+        with the exception of the basin-hopping solver. Each
+        solver has several optional arguments that are not the same across
+        solvers. See the notes section below (or scipy.optimize) for the
+        available arguments and for the list of explicit arguments that the
+        basin-hopping solver supports.
+    maxiter : int, optional
+        The maximum number of iterations to perform.
+    full_output : bool, optional
+        Set to True to have all available output in the Results object's
+        mle_retvals attribute. The output is dependent on the solver.
+        See LikelihoodModelResults notes section for more information.
+    disp : bool, optional
+        Set to True to print convergence messages.
+    callback : callable callback(xk), optional
+        Called after each iteration, as callback(xk), where xk is the
+        current parameter vector.
+    return_params : bool, optional
+        Whether or not to return only the array of maximizing parameters.
+        Default is False.
+    optim_score : {'harvey', 'approx'} or None, optional
+        The method by which the score vector is calculated. 'harvey' uses
+        the method from Harvey (1989), 'approx' uses either finite
+        difference or complex step differentiation depending upon the
+        value of `optim_complex_step`, and None uses the built-in gradient
+        approximation of the optimizer. Default is None. This keyword is
+        only relevant if the optimization method uses the score.
+    optim_complex_step : bool, optional
+        Whether or not to use complex step differentiation when
+        approximating the score; if False, finite difference approximation
+        is used. Default is True. This keyword is only relevant if
+        `optim_score` is set to 'harvey' or 'approx'.
+    optim_hessian : {'opg','oim','approx'}, optional
+        The method by which the Hessian is numerically approximated. 'opg'
+        uses outer product of gradients, 'oim' uses the information
+        matrix formula from Harvey (1989), and 'approx' uses numerical
+        approximation. This keyword is only relevant if the
+        optimization method uses the Hessian matrix.
+    low_memory : bool, optional
+        If set to True, techniques are applied to substantially reduce
+        memory usage. If used, some features of the results object will
+        not be available (including smoothed results and in-sample
+        prediction), although out-of-sample forecasting is possible.
+        Default is False.
+    dynamic : bool, int, str, or datetime, optional
+        Integer offset relative to `start` at which to begin dynamic
+        prediction. Can also be an absolute date string to parse or a
+        datetime type (these are not interpreted as offsets).
+        Prior to this observation, true endogenous values will be used for
+        prediction; starting with this observation and continuing through
+        the end of prediction, forecasted endogenous values will be used
+        instead.
+    information_set : str, optional
+        The information set to condition each prediction on. Default is
+        "predicted", which computes predictions of period t values
+        conditional on observed data through period t-1; these are
+        one-step-ahead predictions, and correspond with the typical
+        `fittedvalues` results attribute. Alternatives are "filtered",
+        which computes predictions of period t values conditional on
+        observed data through period t, and "smoothed", which computes
+        predictions of period t values conditional on the entire dataset
+        (including also future observations t+1, t+2, ...).
+    signal_only : bool, optional
+        Whether to compute predictions of only the "signal" component of
+        the observation equation. Default is False. For example, the
+        observation equation of a time-invariant model is
+        :math:`y_t = d + Z \alpha_t + \varepsilon_t`, and the "signal"
+        component is then :math:`Z \alpha_t`. If this argument is set to
+        True, then predictions of the "signal" :math:`Z \alpha_t` will be
+        returned. Otherwise, the default is for predictions of :math:`y_t`
+        to be returned.
+    suppress_warnings : bool, optional
+        Many warnings might be thrown inside of statsmodels. If
+        ``suppress_warnings`` is True, all of these warnings will be squelched.
+        Default is False.
+
+    Notes
+    -----
+    Generically, the VARMAX model is specified (see for example chapter 18 of
+    [1]_):
+    .. math::
+        y_t = A(t) + A_1 y_{t-1} + \dots + A_p y_{t-p} + B x_t + \epsilon_t +
+        M_1 \epsilon_{t-1} + \dots M_q \epsilon_{t-q}
+    where :math:`\epsilon_t \sim N(0, \Omega)`, and where :math:`y_t` is a
+    `k_endog x 1` vector. Additionally, this model allows considering the case
+    where the variables are measured with error.
+    Note that in the full VARMA(p,q) case there is a fundamental identification
+    problem in that the coefficient matrices :math:`\{A_i, M_j\}` are not
+    generally unique, meaning that for a given time series process there may
+    be multiple sets of matrices that equivalently represent it. See Chapter 12
+    of [1]_ for more information. Although this class can be used to estimate
+    VARMA(p,q) models, a warning is issued to remind users that no steps have
+    been taken to ensure identification in this case.
+
+    References
+    ----------
+    .. [1] Lütkepohl, Helmut. 2007.
+       New Introduction to Multiple Time Series Analysis.
+       Berlin: Springer.
+
+    Examples
+    --------
+    >>> from sktime.forecasting.varmax import VARMAX
+    >>> from sktime.datasets import load_macroeconomic
+    >>> from sktime.forecasting.model_selection import temporal_train_test_split
+    >>> y = load_macroeconomic()
+    >>> forecaster = VARMAX(suppress_warnings=True)
+    >>> forecaster.fit(y[['realgdp', 'unemp']])
+    VARMAX(...)
+    >>> y_pred = forecaster.predict(fh=[1,4,12])
+    """
+
+    _tags = {
+        "scitype:y": "multivariate",
+        "ignores-exogeneous-X": False,
+        "handles-missing-data": False,
+        "y_inner_mtype": "pd.DataFrame",
+        "X_inner_mtype": "pd.DataFrame",
+        "requires-fh-in-fit": False,
+        "X-y-must-have-same-index": True,
+        "enforce_index_type": None,
+        "capability:pred_int": False,
+    }
+
+    def __init__(
+        self,
+        order=(1, 0),
+        trend="c",
+        error_cov_type="unstructured",
+        measurement_error=False,
+        enforce_stationarity=True,
+        enforce_invertibility=True,
+        trend_offset=1,
+        start_params=None,
+        transformed=True,
+        includes_fixed=False,
+        cov_type=None,
+        cov_kwds=None,
+        method="lbfgs",
+        maxiter=50,
+        full_output=1,
+        disp=False,
+        callback=None,
+        return_params=False,
+        optim_score=None,
+        optim_complex_step=None,
+        optim_hessian=None,
+        flags=None,
+        low_memory=False,
+        dynamic=False,
+        information_set="predicted",
+        signal_only=False,
+        suppress_warnings=False,
+    ):
+        # Model parameters
+        self.order = order
+        self.trend = trend
+        self.error_cov_type = error_cov_type
+        self.measurement_error = measurement_error
+        self.enforce_stationarity = enforce_stationarity
+        self.enforce_invertibility = enforce_invertibility
+        self.trend_offset = trend_offset
+        self.start_params = start_params
+        self.transformed = transformed
+        self.includes_fixed = includes_fixed
+        self.cov_type = cov_type
+        self.cov_kwds = cov_kwds
+        self.method = method
+        self.maxiter = maxiter
+        self.full_output = full_output
+        self.disp = disp
+        self.callback = callback
+        self.return_params = return_params
+        self.optim_score = optim_score
+        self.optim_complex_step = optim_complex_step
+        self.optim_hessian = optim_hessian
+        self.flags = flags
+        self.low_memory = low_memory
+        self.dynamic = dynamic
+        self.information_set = information_set
+        self.signal_only = signal_only
+        self.suppress_warnings = suppress_warnings
+
+        super(VARMAX, self).__init__()
+
+    def _fit_forecaster(self, y, X=None):
+        """Fit forecaster to training data.
+
+        Writes to self:
+            Sets fitted model attributes ending in "_".
+
+        Parameters
+        ----------
+        y : array_like
+            The observed time-series process :math:`y`, shaped n_obs x k_endog.
+        X : array_like, optional (default=None)
+            Array of exogenous regressors, shaped n_obs x k.
+
+        Returns
+        -------
+        self : reference to self
+        """
+        if self.suppress_warnings:
+            warnings.filterwarnings("ignore")
+
+        self._forecaster = _VARMAX(
+            endog=y,
+            exog=X,
+            order=self.order,
+            trend=self.trend,
+            error_cov_type=self.error_cov_type,
+            measurement_error=self.measurement_error,
+            enforce_stationarity=self.enforce_stationarity,
+            enforce_invertibility=self.enforce_invertibility,
+            trend_offset=self.trend_offset,
+        )
+        self._fitted_forecaster = self._forecaster.fit(
+            start_params=self.start_params,
+            transformed=self.transformed,
+            includes_fixed=self.includes_fixed,
+            cov_type=self.cov_type,
+            cov_kwds=self.cov_kwds,
+            method=self.method,
+            maxiter=self.maxiter,
+            full_output=self.full_output,
+            disp=self.disp,
+            callback=self.callback,
+            return_params=self.return_params,
+            optim_score=self.optim_score,
+            optim_complex_step=self.optim_complex_step,
+            optim_hessian=self.optim_hessian,
+            flags=self.flags,
+            low_memory=self.low_memory,
+        )
+        return self
+
+    # defining `_predict`, instead of inheriting from `_StatsModelsAdapter`,
+    # for two reasons:
+    # 1. to pass in `dynamic`, `information_set` and `signal_only`
+    # 2. to deal with statsmodel integer indexing issue
+    def _predict(self, fh, X=None):
+        """
+        Wrap Statmodel's VARMAX forecast method.
+
+        Parameters
+        ----------
+        fh : ForecastingHorizon
+            The forecasters horizon with the steps ahead to to predict.
+            Default is one-step ahead forecast,
+            i.e. np.array([1])
+        X : pd.DataFrame, optional (default=None)
+            Exogenous variables.
+
+        Returns
+        -------
+        y_pred : np.ndarray
+            Returns series of predicted values.
+        """
+        start, end = fh.to_absolute_int(self._y.index[0], self.cutoff)[[0, -1]]
+
+        y_pred = self._fitted_forecaster.predict(
+            start=start,
+            end=end,
+            dynamic=self.dynamic,
+            information_set=self.information_set,
+            signal_only=self.signal_only,
+            exog=X,
+        )
+
+        # statsmodel returns zero-based index when index is of type int with the
+        # following warning
+        # ValueWarning: No supported index is available. Prediction results will be
+        # given with an integer index beginning at `start`...
+        # but only when out-of-sample forecasting, i.e. when forecasting horizon is
+        # greater than zero
+        if (type(self._y.index) == pd.core.indexes.numeric.Int64Index) & (
+            any(fh.to_relative(self.cutoff) > 0)
+        ):
+            y_pred.index = y_pred.index + self._y.index[0]
+
+        return y_pred.loc[fh.to_absolute(self.cutoff).to_pandas()]
+
+    @classmethod
+    def get_test_params(cls, parameter_set="default"):
+        """Return testing parameter settings for the estimator.
+
+        Parameters
+        ----------
+        parameter_set : str, default="default"
+            Name of the set of test parameters to return, for use in tests. If no
+            special parameters are defined for a value, will return `"default"` set.
+            There are currently no reserved values for forecasters.
+
+        Returns
+        -------
+        params : dict or list of dict, default = {}
+            Parameters to create testing instances of the class
+            Each dict are parameters to construct an "interesting" test instance, i.e.,
+            `MyClass(**params)` or `MyClass(**params[i])` creates a valid test instance.
+            `create_test_instance` uses the first (or only) dictionary in `params`
+        """
+        params = [
+            {"order": (1, 0)},
+            {"order": (0, 1)},
+            {"order": (1, 1)},
+        ]
+
+        return params
diff --git a/sktime/performance_metrics/base/_base.py b/sktime/performance_metrics/base/_base.py
index 594e01c88..401fffb8f 100644
--- a/sktime/performance_metrics/base/_base.py
+++ b/sktime/performance_metrics/base/_base.py
@@ -5,8 +5,6 @@
 __author__ = ["rnkuhns", "fkiraly"]
 __all__ = ["BaseMetric"]
 
-from warnings import warn
-
 from sktime.base import BaseObject
 
 
@@ -16,22 +14,9 @@ class BaseMetric(BaseObject):
     Extends sktime BaseObject.
     """
 
-    # todo: 0.13.0, remove the func/name args here and to all the metrics
-    def __init__(
-        self,
-        func=None,
-        name=None,
-    ):
-
-        if func is not None or name is not None:
-            warn(
-                "func and name parameters of BaseMetric are deprecated from 0.12.0"
-                "and will be removed in 0.13.0",
-                DeprecationWarning,
-            )
+    def __init__(self):
 
-        self.func = func
-        self.name = name
+        super(BaseMetric, self).__init__()
 
     def __call__(self, y_true, y_pred, **kwargs):
         """Calculate metric value using underlying metric function.
diff --git a/sktime/performance_metrics/forecasting/_classes.py b/sktime/performance_metrics/forecasting/_classes.py
index f4f814dba..3eea4e338 100644
--- a/sktime/performance_metrics/forecasting/_classes.py
+++ b/sktime/performance_metrics/forecasting/_classes.py
@@ -65,21 +65,6 @@ __all__ = [
 ]
 
 
-class _BaseForecastingErrorMetric(BaseMetric):
-    """Base class for defining forecasting error metrics in sktime.
-
-    Extends sktime's BaseMetric to the forecasting interface. Forecasting error
-    metrics measure the error (loss) between forecasts and true values. Lower
-    values are better.
-    """
-
-    greater_is_better = False
-
-    def __init__(self, func, name=None, multioutput="uniform_average"):
-        self.multioutput = multioutput
-        super().__init__(func, name=name)
-
-
 class BaseForecastingErrorMetric(BaseMetric):
     """Base class for defining forecasting error metrics in sktime.
 
@@ -113,30 +98,9 @@ class BaseForecastingErrorMetric(BaseMetric):
         "inner_implements_multilevel": False,
     }
 
-    # todo: 0.13.0, remove the greater_is_better property
-    @property
-    def greater_is_better(self):
-        """Whether greater is better for the metric."""
-        warn(
-            "The greater_is_better attribute is deprecated from 0.12.0 "
-            "and will be removed in 0.13.0. Use the lower_is_better tag instead, "
-            'e.g., my_metric.get_tag("lower_is_better")',
-            DeprecationWarning,
-        )
-        return not self.get_tag("lower_is_better", False)
-
-    def __init__(
-        self,
-        func=None,
-        name=None,
-        multioutput="uniform_average",
-        multilevel="uniform_average",
-    ):
+    def __init__(self, multioutput="uniform_average", multilevel="uniform_average"):
         self.multioutput = multioutput
         self.multilevel = multilevel
-        # todo: in 0.13.0, when func/name are removed, remove func/name args
-        #   temporarily, this is how we avoid overwrite with None
-        # super().__init__(func=func, name=name)
         self.name = type(self).__name__
 
     def __call__(self, y_true, y_pred, **kwargs):
@@ -441,15 +405,25 @@ class BaseForecastingErrorMetricFunc(BaseForecastingErrorMetric):
         return res
 
 
-class _BaseForecastingScoreMetric(_BaseForecastingErrorMetric):
-    """Base class for defining forecasting score metrics in sktime.
+class _DynamicForecastingErrorMetric(BaseForecastingErrorMetricFunc):
+    """Class for defining forecasting error metrics from a function dynamically."""
 
-    Extends sktime's BaseMetric to the forecasting interface. Forecasting score
-    metrics measure the agreement between forecasts and true values. Higher
-    values are better.
-    """
+    def __init__(
+        self,
+        func,
+        name=None,
+        multioutput="uniform_average",
+        multilevel="uniform_average",
+        lower_is_better=True,
+    ):
+        self.multioutput = multioutput
+        self.multilevel = multilevel
+        self.func = func
+        self.name = name
+        super().__init__()
 
-    _tags = {"lower_is_better": False}
+        if not lower_is_better:
+            self.set_tags(**{"lower_is_better": False})
 
 
 class _ScaledMetricTags:
@@ -463,9 +437,13 @@ class _ScaledMetricTags:
 
 
 def make_forecasting_scorer(
-    func, name=None, greater_is_better=False, multioutput="uniform_average"
+    func,
+    name=None,
+    greater_is_better=False,
+    multioutput="uniform_average",
+    multilevel="uniform_average",
 ):
-    """Create a metric class from metric functions.
+    """Create a metric class from a metric functions.
 
     Parameters
     ----------
@@ -483,16 +461,24 @@ def make_forecasting_scorer(
         If array-like, values used as weights to average the errors.
         If 'raw_values', returns a full set of errors in case of multioutput input.
         If 'uniform_average', errors of all outputs are averaged with uniform weight.
+    multilevel : {'raw_values', 'uniform_average'}
+        Defines how to aggregate metric for hierarchical data (with levels).
+        If 'uniform_average' (default), errors are mean-averaged across levels.
+        If 'raw_values', does not average errors across levels, hierarchy is retained.
 
     Returns
     -------
     scorer:
         Metric class that can be used as forecasting scorer.
     """
-    if greater_is_better:
-        return _BaseForecastingScoreMetric(func, name=name, multioutput=multioutput)
-    else:
-        return _BaseForecastingErrorMetric(func, name=name, multioutput=multioutput)
+    lower_is_better = not greater_is_better
+    return _DynamicForecastingErrorMetric(
+        func,
+        name=name,
+        multioutput=multioutput,
+        multilevel=multilevel,
+        lower_is_better=lower_is_better,
+    )
 
 
 class MeanAbsoluteScaledError(_ScaledMetricTags, BaseForecastingErrorMetricFunc):
diff --git a/sktime/performance_metrics/forecasting/probabilistic/_classes.py b/sktime/performance_metrics/forecasting/probabilistic/_classes.py
index d935874ab..210a1566d 100644
--- a/sktime/performance_metrics/forecasting/probabilistic/_classes.py
+++ b/sktime/performance_metrics/forecasting/probabilistic/_classes.py
@@ -8,12 +8,12 @@ import pandas as pd
 from sklearn.utils import check_array, check_consistent_length
 
 from sktime.datatypes import check_is_scitype, convert
-from sktime.performance_metrics.forecasting._classes import _BaseForecastingErrorMetric
+from sktime.performance_metrics.forecasting._classes import BaseForecastingErrorMetric
 
 # TODO: Rework tests now
 
 
-class _BaseProbaForecastingErrorMetric(_BaseForecastingErrorMetric):
+class _BaseProbaForecastingErrorMetric(BaseForecastingErrorMetric):
     """Base class for probabilistic forecasting error metrics in sktime.
 
     Extends sktime's BaseMetric to the forecasting interface. Forecasting error
@@ -39,16 +39,10 @@ class _BaseProbaForecastingErrorMetric(_BaseForecastingErrorMetric):
         "lower_is_better": True,
     }
 
-    def __init__(
-        self,
-        func=None,
-        name=None,
-        multioutput="uniform_average",
-        score_average=True,
-    ):
+    def __init__(self, multioutput="uniform_average", score_average=True):
         self.multioutput = multioutput
         self.score_average = score_average
-        super().__init__(func, name=name, multioutput=multioutput)
+        super().__init__(multioutput=multioutput)
 
     def __call__(self, y_true, y_pred, **kwargs):
         """Calculate metric value using underlying metric function.
@@ -398,19 +392,11 @@ class PinballLoss(_BaseProbaForecastingErrorMetric):
         "lower_is_better": True,
     }
 
-    def __init__(
-        self,
-        multioutput="uniform_average",
-        score_average=True,
-        alpha=None,
-    ):
-        name = "PinballLoss"
+    def __init__(self, multioutput="uniform_average", score_average=True, alpha=None):
         self.score_average = score_average
         self.alpha = self._check_alpha(alpha)
         self.metric_args = {"alpha": alpha}
-        super().__init__(
-            name=name, multioutput=multioutput, score_average=score_average
-        )
+        super().__init__(multioutput=multioutput, score_average=score_average)
 
     def _evaluate_by_index(self, y_true, y_pred, multioutput, **kwargs):
         """Logic for finding the metric evaluated at each index.
@@ -489,12 +475,9 @@ class EmpiricalCoverage(_BaseProbaForecastingErrorMetric):
     }
 
     def __init__(self, multioutput="uniform_average", score_average=True):
-        name = "EmpiricalCoverage"
         self.score_average = score_average
         self.multioutput = multioutput
-        super().__init__(
-            name=name, score_average=score_average, multioutput=multioutput
-        )
+        super().__init__(score_average=score_average, multioutput=multioutput)
 
     def _evaluate_by_index(self, y_true, y_pred, multioutput, **kwargs):
         """Logic for finding the metric evaluated at each index.
@@ -559,12 +542,9 @@ class ConstraintViolation(_BaseProbaForecastingErrorMetric):
     }
 
     def __init__(self, multioutput="uniform_average", score_average=True):
-        name = "ConstraintViolation"
         self.score_average = score_average
         self.multioutput = multioutput
-        super().__init__(
-            name=name, score_average=score_average, multioutput=multioutput
-        )
+        super().__init__(score_average=score_average, multioutput=multioutput)
 
     def _evaluate_by_index(self, y_true, y_pred, multioutput, **kwargs):
         """Logic for finding the metric evaluated at each index.
diff --git a/sktime/transformations/series/difference.py b/sktime/transformations/series/difference.py
index 49da24d2a..f13ebba39 100644
--- a/sktime/transformations/series/difference.py
+++ b/sktime/transformations/series/difference.py
@@ -2,11 +2,10 @@
 # -*- coding: utf-8 -*-
 # copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
 """Class to iteratively apply differences to a time series."""
-__author__ = ["RNKuhns"]
+__author__ = ["RNKuhns", "fkiraly"]
 __all__ = ["Differencer"]
 
 from typing import Union
-from warnings import warn
 
 import numpy as np
 import pandas as pd
@@ -61,9 +60,6 @@ def _inverse_diff(Z, lag):
     return Z
 
 
-# todo: deprecation in 0.13.0
-#   remove the drop_na *argument* and its handling
-#   change the default behaviour form "drop_na" to "fill_zero"
 class Differencer(BaseTransformer):
     """Apply iterative differences to a timeseries.
 
@@ -87,13 +83,7 @@ class Differencer(BaseTransformer):
         The lags used to difference the data.
         If a single `int` value is
 
-    drop_na : bool, default = True
-        deprecated from 0.12.0, to be removed in 0.13.0
-        Whether the differencer should drop the initial observations that
-        contain missing values as a result of the differencing operation(s).
-
-    na_handling : str, default = "drop_na"
-        default will change to "fill_zero" from 0.13.0
+    na_handling : str, default = "fill_zero"
         How to handle the NaNs that appear at the start of the series from differencing
         Example: there are only 3 differences in a series of length 4,
             differencing [a, b, c, d] gives [?, b-a, c-b, d-c]
@@ -128,27 +118,9 @@ class Differencer(BaseTransformer):
 
     VALID_NA_HANDLING_STR = ["drop_na", "keep_na", "fill_zero"]
 
-    def __init__(self, lags=1, drop_na=None, na_handling="drop_na"):
+    def __init__(self, lags=1, na_handling="fill_zero"):
         self.lags = lags
-        self.drop_na = drop_na
         self.na_handling = self._check_na_handling(na_handling)
-        # note: internally, we will use self._na_handling
-        #   because we must never change input param saves for sklearn compatibility
-        #   and because we need to "translate" the old "drop_na" arg
-        #   this could, in certain cases, overwrite the self. parameter
-        #   and cause sklearn compatibility issues due to the point mentioned
-
-        translate_arg = {True: "drop_na", False: "keep_na"}
-
-        if drop_na is not None:
-            warn(
-                f"the drop_na parameter is deprecated and will be removed in 0.13.0, "
-                f'use na_handling="{translate_arg[drop_na]}" instead',
-                DeprecationWarning,
-            )
-            self._na_handling = translate_arg[drop_na]
-        else:
-            self._na_handling = na_handling
 
         self._Z = None
         self._lags = None
@@ -159,7 +131,7 @@ class Differencer(BaseTransformer):
 
         # if the na_handling is "fill_zero" or "keep_na"
         #   then the returned indices are same to the passed indices
-        if self._na_handling in ["fill_zero", "keep_na"]:
+        if self.na_handling in ["fill_zero", "keep_na"]:
             self.set_tags(**{"transform-returns-same-time-index": True})
 
     def _check_na_handling(self, na_handling):
@@ -255,7 +227,7 @@ class Differencer(BaseTransformer):
         """
         Xt = _diff_transform(X, self._lags)
 
-        na_handling = self._na_handling
+        na_handling = self.na_handling
         if na_handling == "drop_na":
             Xt = Xt.iloc[self._cumulative_lags[-1] :]
         elif na_handling == "fill_zero":
diff --git a/sktime/transformations/series/impute.py b/sktime/transformations/series/impute.py
index fb32f61f4..90f4b4ea1 100644
--- a/sktime/transformations/series/impute.py
+++ b/sktime/transformations/series/impute.py
@@ -7,8 +7,6 @@ __author__ = ["aiwalter"]
 __all__ = ["Imputer"]
 
 
-from warnings import warn
-
 import numpy as np
 from sklearn.utils import check_random_state
 
@@ -170,16 +168,6 @@ class Imputer(BaseTransformer):
         """
         X = X.copy()
 
-        # TODO v0.13.0: Remove this if statement and warning
-        if self.method in ["drift", "mean", "median", "random"]:
-            warn(
-                """Imputer methods "drift", "mean", "median", "random" have been
-                moved to the fit() method, so usage via transform() is
-                deprecated. To still fit on the transform data only, please use the
-                new FitInTransform transformer.
-                """,
-                DeprecationWarning,
-            )
         # replace missing_values with np.nan
         if self.missing_values:
             X = X.replace(to_replace=self.missing_values, value=np.nan)
@@ -190,7 +178,7 @@ class Imputer(BaseTransformer):
         if self.method == "random":
             for col in X.columns:
                 X[col] = X[col].apply(
-                    lambda i: self._get_random(col) if np.isnan(i) else i
+                    lambda i: self._get_random(col) if np.isnan(i) else i  # noqa: B023
                 )
         elif self.method == "constant":
             X = X.fillna(value=self.value)
diff --git a/sktime/transformations/series/summarize.py b/sktime/transformations/series/summarize.py
index 6c25f0ffc..9ac75e86c 100644
--- a/sktime/transformations/series/summarize.py
+++ b/sktime/transformations/series/summarize.py
@@ -6,8 +6,6 @@
 __author__ = ["mloning", "RNKuhns", "danbartl", "grzegorzrut"]
 __all__ = ["SummaryTransformer", "WindowSummarizer"]
 
-import warnings
-
 import pandas as pd
 from joblib import Parallel, delayed
 
@@ -152,19 +150,6 @@ class WindowSummarizer(BaseTransformer):
     >>> transformer = WindowSummarizer(**kwargs)
     >>> y_transformed = transformer.fit_transform(y)
 
-        Example where we transform on a different, later test set:
-    >>> y = load_airline()
-    >>> y_train, y_test = temporal_train_test_split(y)
-    >>> kwargs = {
-    ...     "lag_config": {
-    ...         "lag": ["lag", [[1, 0]]],
-    ...         "mean": ["mean", [[3, 0], [12, 0]]],
-    ...         "std": ["std", [[4, 0]]],
-    ...     }
-    ... }
-    >>> transformer = WindowSummarizer(**kwargs)
-    >>> y_test_transformed = transformer.fit(y_train).transform(y_test)
-
         Example with transforming multiple columns of exogeneous features
     >>> y, X = load_longley()
     >>> y_train, y_test, X_train, X_test = temporal_train_test_split(y, X)
@@ -218,15 +203,12 @@ class WindowSummarizer(BaseTransformer):
 
     def __init__(
         self,
-        lag_config=None,
         lag_feature=None,
         n_jobs=-1,
         target_cols=None,
         truncate=None,
     ):
 
-        # self._converter_store_X = dict()
-        self.lag_config = lag_config
         self.lag_feature = lag_feature
         self.n_jobs = n_jobs
         self.target_cols = target_cols
@@ -275,47 +257,32 @@ class WindowSummarizer(BaseTransformer):
             self._target_cols = self.target_cols
 
         # Convert lag config dictionary to pandas dataframe
-        if self.lag_config is not None:
-            func_dict = pd.DataFrame(self.lag_config).T.reset_index()
-            func_dict.rename(
-                columns={"index": "name", 0: "summarizer", 1: "window"},
-                inplace=True,
-            )
-            func_dict = func_dict.explode("window")
-            func_dict["window"] = func_dict["window"].apply(lambda x: [x[1] + 1, x[0]])
-            func_dict.drop("name", inplace=True, axis=1)
-            warnings.warn(
-                "Specifying lag features via lag_config is deprecated since 0.12.0,"
-                + " and will be removed in 0.13.0. Please use the lag_feature notation"
-                + " (see the documentation for the new notation)."
-            )
+        if self.lag_feature is None:
+            func_dict = pd.DataFrame(
+                {
+                    "lag": [1],
+                }
+            ).T.reset_index()
         else:
-            if self.lag_feature is None:
-                func_dict = pd.DataFrame(
-                    {
-                        "lag": [1],
-                    }
-                ).T.reset_index()
-            else:
-                func_dict = pd.DataFrame.from_dict(
-                    self.lag_feature, orient="index"
-                ).reset_index()
+            func_dict = pd.DataFrame.from_dict(
+                self.lag_feature, orient="index"
+            ).reset_index()
 
-            func_dict = pd.melt(
-                func_dict, id_vars="index", value_name="window", ignore_index=False
-            )
-            func_dict.sort_index(inplace=True)
-            func_dict.drop("variable", axis=1, inplace=True)
-            func_dict.rename(
-                columns={"index": "summarizer"},
-                inplace=True,
-            )
-            func_dict = func_dict.dropna(axis=0, how="any")
-            # Identify lags (since they can follow special notation)
-            lags = func_dict["summarizer"] == "lag"
-            # Convert lags to default list notation with window_length 1
-            boost_lag = func_dict.loc[lags, "window"].apply(lambda x: [int(x), 1])
-            func_dict.loc[lags, "window"] = boost_lag
+        func_dict = pd.melt(
+            func_dict, id_vars="index", value_name="window", ignore_index=False
+        )
+        func_dict.sort_index(inplace=True)
+        func_dict.drop("variable", axis=1, inplace=True)
+        func_dict.rename(
+            columns={"index": "summarizer"},
+            inplace=True,
+        )
+        func_dict = func_dict.dropna(axis=0, how="any")
+        # Identify lags (since they can follow special notation)
+        lags = func_dict["summarizer"] == "lag"
+        # Convert lags to default list notation with window_length 1
+        boost_lag = func_dict.loc[lags, "window"].apply(lambda x: [int(x), 1])
+        func_dict.loc[lags, "window"] = boost_lag
         self.truncate_start = func_dict["window"].apply(lambda x: x[0] + x[1]).max()
         self._func_dict = func_dict
 


Test patch:
diff --git a/sktime/forecasting/tests/test_varmax.py b/sktime/forecasting/tests/test_varmax.py
new file mode 100644
index 000000000..1c39e829f
--- /dev/null
+++ b/sktime/forecasting/tests/test_varmax.py
@@ -0,0 +1,66 @@
+# -*- coding: utf-8 -*-
+"""Tests the VARMAX model."""
+__author__ = ["KatieBuc"]
+
+import numpy as np
+import pandas as pd
+from numpy.testing import assert_allclose
+from statsmodels.tsa.api import VARMAX as _VARMAX
+
+from sktime.forecasting.base import ForecastingHorizon
+from sktime.forecasting.model_selection import temporal_train_test_split
+from sktime.forecasting.varmax import VARMAX
+
+np.random.seed(13455)
+index = pd.date_range(start="2020-01", end="2021-12", freq="M")
+df = pd.DataFrame(
+    np.random.randint(0, 100, size=(23, 3)),
+    columns=list("ABC"),
+    index=pd.PeriodIndex(index),
+)
+
+
+def test_VARMAX_against_statsmodels():
+    """Compares Sktime's and Statsmodel's VARMAX.
+
+    with default variables.
+    """
+    train, _ = temporal_train_test_split(df.astype("float64"))
+    y = train[["A", "B"]]
+
+    sktime_model = VARMAX()
+    fh = ForecastingHorizon([1, 3, 4, 5, 7, 9])
+    sktime_model.fit(y)
+    y_pred = sktime_model.predict(fh=fh)
+
+    stats = _VARMAX(y)
+    stats_fit = stats.fit()
+    start, end = len(train) + fh[0] - 1, len(train) + fh[-1] - 1
+    y_pred_stats = stats_fit.predict(start=start, end=end)
+    y_pred_stats = y_pred_stats.loc[fh.to_absolute(train.index[-1]).to_pandas()]
+
+    assert_allclose(y_pred, y_pred_stats)
+
+
+def test_VARMAX_against_statsmodels_with_exog():
+    """Compares Sktime's and Statsmodel's VARMAX.
+
+    with exogenous input.
+    """
+    train, test = temporal_train_test_split(df.astype("float64"))
+    y_train, X_train = train[["A", "B"]], train[["C"]]
+    _, X_test = test[["A", "B"]], test[["C"]]
+    fh = ForecastingHorizon([1, 2, 3, 4, 5, 6])
+    assert len(fh) == len(X_test)
+
+    sktime_model = VARMAX()
+    sktime_model.fit(y_train, X=X_train)
+    y_pred = sktime_model.predict(fh=fh, X=X_test)
+
+    stats = _VARMAX(y_train, exog=X_train)
+    stats_fit = stats.fit()
+    start, end = len(train) + fh[0] - 1, len(train) + fh[-1] - 1
+    y_pred_stats = stats_fit.predict(start=start, end=end, exog=X_test)
+    y_pred_stats = y_pred_stats.loc[fh.to_absolute(train.index[-1]).to_pandas()]
+
+    assert_allclose(y_pred, y_pred_stats)