Repository: sktime/sktime
Pull request #8242 (2025-05-05 21:15:33)

[MNT] [Dependabot](deps): Update pytest-timeout requirement from <2.4,>=2.1 to >=2.1,<2.5
Updates the requirements on [pytest-timeout](https://github.com/pytest-dev/pytest-timeout) to permit the latest version.
<details>
<summary>Commits</summary>
<ul>
<li><a href="https://github.com/pytest-dev/pytest-timeout/commit/ddabc934535081a5bf9ba7c9ca5b494aeaf8f665"><code>ddabc93</code></a> Add python 3.13</li>
<li><a href="https://github.com/pytest-dev/pytest-timeout/commit/617b1b6f0ab74517b92a4966b6db9e509d987972"><code>617b1b6</code></a> Prepare release</li>
<li><a href="https://github.com/pytest-dev/pytest-timeout/commit/2e96621ee1b61057438c0f8c5b158eddc31654c5"><code>2e96621</code></a> Change error message to clearly relate it to pytest-timeout (<a href="https://redirect.github.com/pytest-dev/pytest-timeout/issues/180">#180</a>)</li>
<li><a href="https://github.com/pytest-dev/pytest-timeout/commit/733b26b5d3b6f1c86876339f336cd827dc7d1a0c"><code>733b26b</code></a> Update changelog</li>
<li><a href="https://github.com/pytest-dev/pytest-timeout/commit/fa5fd451afa3e7e676c7a31a9e7b058dabf0988b"><code>fa5fd45</code></a> Add support sys.monitoring</li>
<li><a href="https://github.com/pytest-dev/pytest-timeout/commit/3d41ba8042ab1d07613572d95853dcc1c3602880"><code>3d41ba8</code></a> Update README.rst with pytest-xdist 3.0.2 news</li>
<li><a href="https://github.com/pytest-dev/pytest-timeout/commit/935f51cf0931dbbee8162e89b53289f1f43df499"><code>935f51c</code></a> Add missing Python 3.12 trove classifier</li>
<li><a href="https://github.com/pytest-dev/pytest-timeout/commit/04432f5b10eec31433482e1b191c0d71f19acff5"><code>04432f5</code></a> Some fixes to readme, bump version</li>
<li><a href="https://github.com/pytest-dev/pytest-timeout/commit/48179d984d4a716b426dbf04e6b290d620518743"><code>48179d9</code></a> Prep release</li>
<li><a href="https://github.com/pytest-dev/pytest-timeout/commit/38c5f24e73cb73ea8a368024bf171921e8e500e7"><code>38c5f24</code></a> Tweak docs a little</li>
<li>Additional commits viewable in <a href="https://github.com/pytest-dev/pytest-timeout/compare/2.1.0...2.4.0">compare view</a></li>
</ul>
</details>
<br />


Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting `@dependabot rebase`.

[//]: # (dependabot-automerge-start)
[//]: # (dependabot-automerge-end)

---

<details>
<summary>Dependabot commands and options</summary>
<br />

You can trigger Dependabot actions by commenting on this PR:
- `@dependabot rebase` will rebase this PR
- `@dependabot recreate` will recreate this PR, overwriting any edits that have been made to it
- `@dependabot merge` will merge this PR after your CI passes on it
- `@dependabot squash and merge` will squash and merge this PR after your CI passes on it
- `@dependabot cancel merge` will cancel a previously requested merge and block automerging
- `@dependabot reopen` will reopen this PR if it is closed
- `@dependabot close` will close this PR and stop Dependabot recreating it. You can achieve the same result by closing it manually
- `@dependabot show <dependency name> ignore conditions` will show all of the ignore conditions of the specified dependency
- `@dependabot ignore this major version` will close this PR and stop Dependabot creating any more for this major version (unless you reopen the PR or upgrade to it yourself)
- `@dependabot ignore this minor version` will close this PR and stop Dependabot creating any more for this minor version (unless you reopen the PR or upgrade to it yourself)
- `@dependabot ignore this dependency` will close this PR and stop Dependabot creating any more for this dependency (unless you reopen the PR or upgrade to it yourself)


</details>


Patch:
diff --git a/.all-contributorsrc b/.all-contributorsrc
index c132b8812..8e3017487 100644
--- a/.all-contributorsrc
+++ b/.all-contributorsrc
@@ -3158,7 +3158,7 @@
       ]
     },
     {
-      "login": "Garve",
+      "login": "RobKuebler",
       "name": "Robert Kübler",
       "avatar_url": "https://avatars.githubusercontent.com/u/932327?v=4",
       "profile": "https://linkedin.com/in/robert-kuebler",
diff --git a/docs/source/api_reference/forecasting.rst b/docs/source/api_reference/forecasting.rst
index 971275ed5..7ff4815b8 100644
--- a/docs/source/api_reference/forecasting.rst
+++ b/docs/source/api_reference/forecasting.rst
@@ -549,6 +549,14 @@ Pre-trained and foundation models
 
     HFTransformersForecaster
 
+.. currentmodule:: sktime.forecasting.hf_momentfm_forecaster
+
+.. autosummary::
+    :toctree: auto_generated/
+    :template: class.rst
+
+    MomentFMForecaster
+
 .. currentmodule:: sktime.forecasting.chronos
 
 .. autosummary::
diff --git a/pyproject.toml b/pyproject.toml
index d6d91b184..8e53516af 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -266,7 +266,7 @@ docs = [
 tests = [
   "pytest>=7.4,<8.4",
   "pytest-randomly>=3.15,<3.17",
-  "pytest-timeout>=2.1,<2.4",
+  "pytest-timeout>=2.1,<2.5",
   "pytest-xdist>=3.3,<3.7",
 ]
 
@@ -299,6 +299,8 @@ dl = [
   'gluonts>=0.14.3; python_version < "3.12"',
   'einops>0.7.0; python_version < "3.12"',
   'huggingface-hub>=0.23.0; python_version < "3.12"',
+  'accelerate',
+  'tqdm',
   'hydra-core; python_version < "3.13"',
 ]
 mlflow = [
diff --git a/sktime/forecasting/dummy.py b/sktime/forecasting/dummy.py
index eb190cf87..602ace83d 100644
--- a/sktime/forecasting/dummy.py
+++ b/sktime/forecasting/dummy.py
@@ -1,7 +1,7 @@
 # copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
 """Dummy forecasters."""
 
-__author__ = ["fkiraly"]
+__author__ = ["fkiraly", "RobKuebler"]
 
 import pandas as pd
 
@@ -61,7 +61,7 @@ class ForecastKnownValues(BaseForecaster):
     _tags = {
         # packaging info
         # --------------
-        "authors": ["fkiraly"],
+        "authors": ["fkiraly", "RobKuebler"],
         # estimator type
         # --------------
         "y_inner_mtype": "pd.DataFrame",
@@ -155,16 +155,21 @@ class ForecastKnownValues(BaseForecaster):
         fh_abs = fh.to_absolute_index(self.cutoff)
 
         try:
+            idx = self._y.index
+            if isinstance(idx, pd.MultiIndex):
+                unique_levels = idx.droplevel(-1).unique()
+                fh_abs = pd.MultiIndex.from_tuples(
+                    ((*level, time) for level in unique_levels for time in fh_abs),
+                    name=idx.names,
+                )
+
             y_pred = self._y_known.reindex(fh_abs, **reindex_params)
             y_pred = y_pred.reindex(self._y.columns, axis=1, **reindex_params)
         # TypeError happens if indices are incompatible types
         except TypeError:
-            if self.fill_value is None:
-                y_pred = pd.DataFrame(index=fh_abs, columns=self._y.columns)
-            else:
-                y_pred = pd.DataFrame(
-                    self.fill_value, index=fh_abs, columns=self._y.columns
-                )
+            y_pred = pd.DataFrame(
+                self.fill_value, index=fh_abs, columns=self._y.columns
+            )
 
         return y_pred
 
diff --git a/sktime/forecasting/hf_momentfm_forecaster.py b/sktime/forecasting/hf_momentfm_forecaster.py
new file mode 100644
index 000000000..d3b08e1d0
--- /dev/null
+++ b/sktime/forecasting/hf_momentfm_forecaster.py
@@ -0,0 +1,785 @@
+"""Interface for the momentfm deep learning time series forecaster."""
+
+import warnings
+
+import numpy as np
+import pandas as pd
+from skbase.utils.dependencies import _check_soft_dependencies
+
+from sktime.forecasting.base import ForecastingHorizon, _BaseGlobalForecaster
+from sktime.split import temporal_train_test_split
+
+if _check_soft_dependencies("torch", severity="none"):
+    from torch.cuda import empty_cache
+    from torch.utils.data import Dataset
+else:
+
+    class Dataset:
+        """Dummy class if torch is unavailable."""
+
+        pass
+
+
+if _check_soft_dependencies("accelerate", severity="none"):
+    pass
+
+if _check_soft_dependencies("transformers", severity="none"):
+    from sktime.libs.momentfm import MOMENTPipeline
+
+
+class MomentFMForecaster(_BaseGlobalForecaster):
+    """
+    Interface for forecasting with the deep learning time series model momentfm.
+
+    MomentFM is a collection of open source foundation models for the general
+    purpose of time series analysis. The Moment Foundation Model is a pre-trained
+    model that is capable of accomplishing various time series tasks, such as:
+        - Long Term Forecasting
+        - Short Term Forecasting
+
+    This interface with MomentFM focuses on the forecasting task, in which the
+    foundation model uses a user fine tuned 'forecasting head' to predict h steps ahead.
+    This model does NOT have zero shot capabilities and requires fine-tuning
+    to achieve performance on user inputted data.
+
+    For more information: see
+    https://github.com/moment-timeseries-foundation-model/moment
+
+    For information regarding licensing and use of the momentfm model please visit:
+    https://huggingface.co/AutonLab/MOMENT-1-large
+
+    pretrained_model_name_or_path : str
+        Path to the pretrained Momentfm model. Default is AutonLab/MOMENT-1-large
+
+    freeze_encoder : bool
+        Selection of whether or not to freeze the weights of the encoder
+        Default = True
+
+    freeze_embedder : bool
+        Selection whether or not to freeze the patch embedding layer
+        Default = True
+
+    freeze_head : bool
+        Selection whether or not to freeze the forecasting head.
+        Recommendation is that the linear forecasting head must be trained
+        Default = False
+
+    dropout : float
+        Dropout value of the model. Values range between [0.0, 1.0]
+        Default = 0.1
+
+    head_dropout : float
+        Dropout value of the forecasting head. Values range between [0.0, 1.0]
+        Default = 0.1
+
+    seq_len : int
+        length of sequences or length of historical values that are passed
+        to the model for training at each time point. the momentfm model requires
+        sequence lengths to be 512 exactly, so if less, padding will be used.
+        If the sequence length is > 512, it will be reduced to 512.
+        default = 512
+
+    batch_size : int
+        size of batches to train the model on
+        default = 32
+
+    eval_batch_size : int or "all"
+        size of batches to evaluate the model on. If the string "all" is
+        specified, then we process the entire validation set as a single batch
+        default = 32
+
+    epochs : int
+        Number of epochs to fit tune the model on
+        default = 1
+
+    max_lr : float
+        Maximum learning rate that the learning rate scheduler will use
+        default = 1e-4
+
+    device : str
+        torch device to use
+        default = "auto"
+        If set to auto, it will automatically use whatever device that
+        `accelerate` detects.
+
+    pct_start : float
+        percentage of total iterations where the learning rate rises during
+        one epoch
+        default = 0.3
+
+    max_norm : float
+        Float value used to clip gradients during training
+        default = 5.0
+
+    train_val_split : float
+        float value between 0 and 1 to determine portions of training
+        and validation splits
+        default = 0.2
+
+    transformer_backbone : str
+        d_model of a pre-trained transformer model to use. See
+        SUPPORTED_HUGGINGFACE_MODELS to specify valid models to use.
+        Default is 'google/flan-t5-large'.
+
+    config : dict, default = {}
+        If desired, user can pass in a config detailing all momentfm parameters
+        that they wish to set in dictionary form, so that parameters do not need
+        to be individually set. If a parameter inside a config is a
+        duplicate of one already passed in individually, it will be overwritten.
+
+    criterion : criterion, default = torch.nn.MSELoss
+        Criterion to use during training.
+
+    return_model_to_cpu : bool, default = False
+        After fitting and training, will return the `momentfm` model to the cpu.
+
+    References
+    ----------
+    Paper: https://arxiv.org/abs/2402.03885
+    Github: https://github.com/moment-timeseries-foundation-model/moment/tree/main
+
+    Examples
+    --------
+    >>> from sktime.forecasting.hf_momentfm_forecaster import MomentFMForecaster
+    >>> from sktime.datasets import load_airline
+    >>> y = load_airline()
+    >>> forecaster = MomentFMForecaster(seq_len = 2)
+    >>> forecaster.fit(y, fh=[1, 2, 3]) # doctest: +SKIP
+    >>> y_pred = forecaster.predict(y = y) # doctest: +SKIP
+    """
+
+    _tags = {
+        "scitype:y": "both",
+        "authors": ["julian-fong"],
+        "maintainers": ["julian-fong"],
+        "handles-missing-data": False,
+        "y_inner_mtype": [
+            "pd.DataFrame",
+            "pd-multiindex",
+            "pd_multiindex_hier",
+        ],
+        "ignores-exogeneous-X": True,
+        "requires-fh-in-fit": True,
+        "python_dependencies": [
+            "torch",
+            "tqdm",
+            "huggingface-hub",
+            # "momentfm",
+            "accelerate",
+            "transformers",
+        ],
+        "capability:global_forecasting": True,
+        "python_version": ">= 3.10",
+        "capability:insample": False,
+        "capability:pred_int:insample": False,
+        "capability:pred_int": False,
+    }
+
+    def __init__(
+        self,
+        pretrained_model_name_or_path="AutonLab/MOMENT-1-large",
+        freeze_encoder=True,
+        freeze_embedder=True,
+        freeze_head=False,
+        dropout=0.1,
+        head_dropout=0.1,
+        seq_len=512,
+        batch_size=32,
+        eval_batch_size=32,
+        epochs=1,
+        max_lr=1e-4,
+        device="auto",
+        pct_start=0.3,
+        max_norm=5.0,
+        train_val_split=0.2,
+        transformer_backbone="google/flan-t5-large",
+        criterion=None,
+        config=None,
+        return_model_to_cpu=False,
+    ):
+        super().__init__()
+        from torch.nn import MSELoss
+
+        self.pretrained_model_name_or_path = pretrained_model_name_or_path
+        self.freeze_encoder = freeze_encoder
+        self.freeze_embedder = freeze_embedder
+        self.freeze_head = freeze_head
+        self.dropout = dropout
+        self.head_dropout = head_dropout
+        self.seq_len = seq_len
+        self.batch_size = batch_size
+        self.eval_batch_size = eval_batch_size
+        self.epochs = epochs
+        self.max_lr = max_lr
+        self.device = device
+        self.pct_start = pct_start
+        self.max_norm = max_norm
+        self.train_val_split = train_val_split
+        self.transformer_backbone = transformer_backbone
+        self.config = config
+        self._config = config if config is not None else {}
+        self.criterion = criterion
+        self._criterion = self.criterion if self.criterion else MSELoss()
+        self._moment_seq_len = 512
+        self.return_model_to_cpu = return_model_to_cpu
+
+    def _fit(self, fh, y, X=None):
+        """Assumes y is a single or multivariate time series."""
+        from accelerate import Accelerator
+        from torch.optim import Adam
+        from torch.optim.lr_scheduler import OneCycleLR
+        from torch.utils.data import DataLoader
+
+        # keep a copy of y in case y is None in predict
+        self._y = y
+        self._y_index = self._y.index
+        self._y_cols = self._y.columns
+        self._y_shape = self._y.values.shape
+
+        self._pretrained_model_name_or_path = self._config.get(
+            "pretrained_model_name_or_path", self.pretrained_model_name_or_path
+        )
+        self._freeze_encoder = self._config.get("freeze_encoder", self.freeze_encoder)
+        self._pretrained_model_name_or_path = self._config.get(
+            "_pretrained_model_name_or_path", self._pretrained_model_name_or_path
+        )
+        self._freeze_embedder = self._config.get(
+            "freeze_embedder", self.freeze_embedder
+        )
+        self._freeze_head = self._config.get("freeze_head", self.freeze_head)
+        self._dropout = self._config.get("dropout", self.dropout)
+        self._head_dropout = self._config.get("head_dropout", self.head_dropout)
+        self._transformer_backbone = self._config.get(
+            "transformer_backbone", self.transformer_backbone
+        )
+        self._criterion = self._config.get("criterion", self._criterion)
+        self._seq_len = self._config.get("seq_len", self.seq_len)
+
+        if self._seq_len > self._moment_seq_len:
+            warnings.warn(
+                f"length of {self._seq_len} was found which is greater than 512. "
+                "The most recent 512"
+                " values will be used when fitting.",
+                stacklevel=2,
+            )
+            self._seq_len = 512
+        # in the case the config contains 'forecast_horizon', we'll set
+        # fh as that, otherwise we override it using the fh param
+
+        self._fh_config = self._config.get("forecast_horizon", None)
+
+        # device initialization
+        self._device = self._config.get("device", self.device)
+
+        # check availability of user specified device
+        self._device = _check_device(self._device)
+        # initialize accelerator
+        accelerator = Accelerator()
+        if self._device == "auto":
+            self._device = accelerator.device
+
+        cur_epoch = 0
+        max_epoch = self.epochs
+        if fh is not None:
+            self._fh_input = max(fh.to_relative(self.cutoff))
+        self._fh = self._fh_input if fh is not None else self._fh_config
+        # self._model_fh is guaranteed to be an int value
+        self._model_fh = self._fh
+        # revert self._fh back to fh to pass checks
+        self._fh = fh
+
+        self.model = MOMENTPipeline.from_pretrained(
+            self._pretrained_model_name_or_path,
+            model_kwargs={
+                "task_name": "forecasting",
+                "dropout": self._dropout,
+                "head_dropout": self._head_dropout,
+                "freeze_encoder": self._freeze_encoder,
+                "freeze_embedder": self._freeze_embedder,
+                "seq_len": 512,  # forced to be hard coded
+                "freeze_head": self._freeze_head,
+                "device": self._device,
+                "transformer_backbone": self._transformer_backbone,
+                "forecast_horizon": self._model_fh,
+            },
+        )
+        self.model.init()
+        # preparing the datasets
+        y_train, y_test = temporal_train_test_split(
+            y, train_size=1 - self.train_val_split, test_size=self.train_val_split
+        )
+
+        train_dataset = MomentPytorchDataset(
+            y=y_train,
+            fh=self._model_fh,
+            seq_len=self._seq_len,
+            device=self._device,
+        )
+
+        train_dataloader = DataLoader(
+            train_dataset, batch_size=self.batch_size, shuffle=True
+        )
+
+        if not y_test.empty:
+            val_dataset = MomentPytorchDataset(
+                y=y_test,
+                fh=self._model_fh,
+                seq_len=self._seq_len,
+                device=self._device,
+            )
+
+            if self.eval_batch_size == "all":
+                self._eval_batch_size = len(val_dataset)
+            else:
+                self._eval_batch_size = self.eval_batch_size
+
+            val_dataloader = DataLoader(
+                val_dataset, batch_size=self._eval_batch_size, shuffle=True
+            )
+        else:
+            val_dataloader = None
+
+        criterion = self._criterion
+        optimizer = Adam(self.model.parameters(), lr=self.max_lr)
+        # Enable mixed precision training
+
+        # Create a OneCycleLR scheduler
+        total_steps = len(train_dataloader) * max_epoch
+        scheduler = OneCycleLR(
+            optimizer,
+            max_lr=self.max_lr,
+            total_steps=total_steps,
+            pct_start=self.pct_start,
+        )
+
+        # Gradient clipping value
+        max_norm = self.max_norm
+
+        self.model, optimizer, train_dataloader, val_dataloader, scheduler = (
+            accelerator.prepare(
+                self.model, optimizer, train_dataloader, val_dataloader, scheduler
+            )
+        )
+
+        while cur_epoch < max_epoch:
+            cur_epoch = _run_epoch(
+                cur_epoch,
+                accelerator,
+                criterion,
+                optimizer,
+                scheduler,
+                self.model,
+                max_norm,
+                train_dataloader,
+                val_dataloader,
+            )
+
+        if self.return_model_to_cpu:
+            self.model.to("cpu")
+            empty_cache()
+
+        return self
+
+    def _predict(self, y, X=None, fh=None):
+        """Predict method to forecast timesteps into the future.
+
+        fh should not be passed here and
+        must be the same length as the one used to fit the model.
+        """
+        # use y values from fit if y is None in predict
+        if y is None:
+            y = self._y
+
+        index = self._fh.to_absolute_index(self.cutoff)
+        from torch import from_numpy
+
+        self.model = self.model.to(self._device)
+        self.model.eval()
+        y_index_names = list(y.index.names)
+        if isinstance(y.index, pd.MultiIndex):
+            y_ = _frame2numpy(y)
+        else:
+            y_ = np.expand_dims(y.values, axis=0)
+
+        num_instances, sequence_length, num_channels = (
+            y_.shape
+        )  # shape of our input to predict
+        # raise warning if sequence length of y is greater than the sequence
+        # length used to fit the model
+        if sequence_length > self._seq_len:
+            warnings.warn(
+                f"Sequence length of {sequence_length} was found which is greater "
+                "than sequence "
+                f"length {self._seq_len} used to fit the model. The most recent"
+                f" {self._seq_len} values will be used.",
+                stacklevel=2,
+            )
+            # truncate code
+            # only retain the most recent self._seq_len values if greater than
+            # self._seq_len
+            y_ = y_[:, -self._seq_len :, :]
+            sequence_length = self._seq_len
+        if sequence_length < self._moment_seq_len:
+            # if smaller, need to pad values
+            y_ = _create_padding(
+                y_,
+                (num_instances, self._moment_seq_len - sequence_length, num_channels),
+            )
+            input_mask = _create_mask(
+                sequence_length, self._moment_seq_len - sequence_length
+            )
+        else:  # this means sequence_length = self._seq_len == 512
+            input_mask = _create_mask(self._moment_seq_len)
+        if num_channels != self._y_shape[1]:
+            raise ValueError(
+                "The number of multivariate time series "
+                f"{num_channels} passed in predict does not match the "
+                f"number of multivariate time series {self._y_shape[1]} "
+                "used to train the model in fit."
+            )
+            pass
+        # transpose it to change it into (C, S) size
+        y_ = y_.transpose(0, 2, 1)
+        # returns a timeseriesoutput object
+        y_torch_input = from_numpy(y_).float().to(self._device)
+        input_mask = input_mask.to(self._device)
+        output = self.model(x_enc=y_torch_input, mask=input_mask)
+        forecast_output = output.forecast
+        # forecast_output = forecast_output.squeeze(0)
+
+        pred = forecast_output.detach().cpu().numpy()
+        # revert back to (B, S, C)
+        pred = np.transpose(pred, (0, 2, 1))
+
+        if isinstance(y.index, pd.MultiIndex):
+            ins = np.array(list(np.unique(y.index.droplevel(-1)).repeat(pred.shape[1])))
+            ins = [ins[..., i] for i in range(ins.shape[-1])] if ins.ndim > 1 else [ins]
+
+            idx = (
+                ForecastingHorizon(range(1, pred.shape[1] + 1), freq=self.fh.freq)
+                .to_absolute(self._cutoff)
+                ._values.tolist()
+                * pred.shape[0]
+            )
+            index = pd.MultiIndex.from_arrays(
+                ins + [idx],
+                names=y.index.names,
+            )
+        else:
+            index = (
+                ForecastingHorizon(range(1, pred.shape[1] + 1))
+                .to_absolute(self._cutoff)
+                ._values
+            )
+
+        df_pred = pd.DataFrame(
+            # batch_size * num_timestams, n_cols
+            pred.reshape(-1, pred.shape[-1]),
+            index=index,
+            columns=self._y_cols,
+        )
+
+        absolute_horizons = fh.to_absolute_index(self.cutoff)
+        dateindex = df_pred.index.get_level_values(-1).map(
+            lambda x: x in absolute_horizons
+        )
+        df_pred = df_pred.loc[dateindex]
+        df_pred.index.names = y_index_names
+
+        if self.return_model_to_cpu:
+            self.model.to("cpu")
+            empty_cache()
+
+        return df_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.
+
+        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_set = []
+        params1 = {"seq_len": 2, "return_model_to_cpu": True, "train_val_split": 0.0}
+        params_set.append(params1)
+        params2 = {
+            "batch_size": 16,
+            "seq_len": 2,
+            "return_model_to_cpu": True,
+            "train_val_split": 0.0,
+        }
+        params_set.append(params2)
+
+        return params_set
+
+
+def _create_padding(x, pad_shape):
+    """Return zero padded tensor of size seq_len, num_cols."""
+    # For example, if num_rows = 500 and seq_len = 512
+    # then x.shape[0] = 500 and pad_shape[0] = 12
+    # then cat(x, zero_pad) should return (num_cols,512)
+    from torch import cat, zeros
+
+    if isinstance(x, np.ndarray):
+        from torch import from_numpy
+
+        x_ = from_numpy(x)
+    else:
+        x_ = x
+
+    pad_shape_dim = x_.dim()
+    # create a padding of size 2 for fit()
+    if pad_shape_dim == 2:
+        zero_pad = zeros(pad_shape)
+        axis = 0
+    else:
+        # if its 3 dimensions, we are in predict()
+        zero_pad = zeros(pad_shape)
+        axis = 1
+
+    out = cat((x_, zero_pad), axis=axis)
+    if isinstance(x, np.ndarray):
+        out = np.array(out)
+    return out
+
+
+def _create_mask(ones_length, zeros_length=0):
+    from torch import cat, ones, zeros
+
+    zeros_tensor = zeros(zeros_length)
+    ones_tensor = ones(ones_length)
+
+    input_mask = cat((ones_tensor, zeros_tensor))
+    return input_mask
+
+
+def _run_epoch(
+    cur_epoch,
+    accelerator,
+    criterion,
+    optimizer,
+    scheduler,
+    model,
+    max_norm,
+    train_dataloader,
+    val_dataloader,
+):
+    import torch.cuda.amp
+    from tqdm import tqdm
+
+    from sktime.libs.momentfm.utils.forecasting_metrics import get_forecasting_metrics
+
+    losses = []
+    for data in tqdm(train_dataloader, total=len(train_dataloader)):
+        # Move the data to the GPU
+        timeseries = data["historical_y"]
+        input_mask = data["input_mask"]
+        forecast = data["future_y"]
+        with torch.cuda.amp.autocast():
+            output = model(x_enc=timeseries, mask=input_mask)
+        loss = criterion(output.forecast, forecast)
+
+        accelerator.backward(loss)
+
+        # Clip gradients
+        torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm)
+
+        optimizer.step()
+        optimizer.zero_grad(set_to_none=True)
+
+        losses.append(loss.item())
+
+    losses = np.array(losses)
+    average_loss = np.average(losses)
+    tqdm.write(f"Epoch {cur_epoch}: Train loss: {average_loss:.3f}")
+
+    # Step the learning rate scheduler
+    scheduler.step()
+
+    # Evaluate the model on the test split
+    if val_dataloader:
+        trues, preds, histories, losses = [], [], [], []
+        model.eval()
+        with torch.no_grad():
+            for data in tqdm(val_dataloader, total=len(val_dataloader)):
+                # Move the data to the specified device
+                timeseries = data["historical_y"]
+                input_mask = data["input_mask"]
+                forecast = data["future_y"]
+
+                with torch.cuda.amp.autocast():
+                    output = model(x_enc=timeseries, mask=input_mask)
+
+                loss = criterion(output.forecast, forecast)
+                losses.append(loss.item())
+
+                trues.append(forecast.detach().cpu().numpy())
+                preds.append(output.forecast.detach().cpu().numpy())
+                histories.append(timeseries.detach().cpu().numpy())
+
+        losses = np.array(losses)
+        average_loss = np.average(losses)
+        model.train()
+
+        trues = np.concatenate(trues, axis=0)
+        preds = np.concatenate(preds, axis=0)
+        histories = np.concatenate(histories, axis=0)
+
+        metrics = get_forecasting_metrics(y=trues, y_hat=preds, reduction="mean")
+        tqdm.write(
+            f"Epoch {cur_epoch}: Test MSE: {metrics.mse:.3f}Test MAE: {metrics.mae:.3f}"
+        )
+    cur_epoch += 1
+    return cur_epoch
+
+
+def _check_device(device):
+    if device == "auto":
+        return device
+    mps = False
+    cuda = False
+    if device == "mps":
+        from torch.backends.mps import is_available, is_built
+
+        if not is_available():
+            if not is_built():
+                print(
+                    "MPS not available because the current PyTorch install was not "
+                    "built with MPS enabled."
+                )
+            else:
+                print(
+                    "MPS not available because the current MacOS version is not 12.3+ "
+                    "and/or you do not have an MPS-enabled device on this machine."
+                )
+        else:
+            _device = "mps"
+            mps = True
+    elif device == "gpu" or device == "cuda":
+        from torch.cuda import is_available
+
+        if is_available():
+            _device = "cuda"
+            cuda = True
+    elif "cuda" in device:  # for specific cuda devices like cuda:0 etc
+        from torch.cuda import is_available
+
+        if is_available():
+            _device = device
+            cuda = True
+    if mps or cuda:
+        return _device
+    else:
+        _device = "cpu"
+
+    return _device
+
+
+def _same_index(data: pd.DataFrame):
+    data = data.groupby(level=list(range(len(data.index.levels) - 1))).apply(
+        lambda x: x.index.get_level_values(-1)
+    )
+    assert data.map(lambda x: x.equals(data.iloc[0])).all(), (
+        "All series must has the same index"
+    )
+    return data.iloc[0], len(data.iloc[0])
+
+
+def _frame2numpy(data: pd.DataFrame):
+    idx, length = _same_index(data)
+    arr = np.array(data.values, dtype=np.float32).reshape(
+        (-1, length, len(data.columns))
+    )
+    return arr
+
+
+class MomentPytorchDataset(Dataset):
+    """Customized Pytorch dataset for the momentfm model."""
+
+    def __init__(self, y, fh, seq_len, device):
+        self.y = y
+        self.moment_seq_len = 512  # forced seq_len by pre-trained model
+        self.seq_len = seq_len
+        self.fh = fh
+        self.shape = y.shape
+        self.device = device
+
+        # multi-index conversion
+        if isinstance(y.index, pd.MultiIndex):
+            self.y = _frame2numpy(y)
+        else:
+            self.y = np.expand_dims(y.values, axis=0)
+
+        # n_timestamps should be the seq length for a single series in both
+        # cases, multivariate dataframe
+        # or a panel/hierarchical dataset
+        # if hier/panel, self.n_sequences will be > 1, else will be = 1 if
+        # # a regular multivariate df
+        self.n_sequences, self.n_timestamps, self.n_columns = self.y.shape
+
+        # self.single_length is defined as the length of one series under
+        # one instance in the panel/hier case
+        # else it is just the seq_len of the multivariate data
+        self.single_length = self.n_timestamps - self.seq_len - self.fh + 1
+
+        # code block to figure out masking sizes in case seq_len < 512
+        if self.seq_len < self.moment_seq_len:
+            self._pad_shape = (self.moment_seq_len - self.seq_len, self.n_columns)
+            self.input_mask = _create_mask(
+                self.seq_len, self.moment_seq_len - self.seq_len
+            ).float()
+            # Concatenate the tensors
+        elif self.seq_len > self.moment_seq_len:
+            # for now if seq_len > 512 than we reduce it back to 512
+            self.seq_len = 512
+            self.input_mask = _create_mask(self.seq_len).float()
+        else:
+            self.input_mask = _create_mask(self.seq_len).float()
+
+    def __len__(self):
+        """Return length of dataset."""
+        # in the case of a regular multivariate df, we just return the trivial case
+        # self.n_timestamps - self.seq_len - self.fh + 1
+        # but in case the data is panel/hier, we need to count the total
+        # #number of instances
+        # i.e self.n_sequences
+        return self.single_length * self.n_sequences
+
+    def __getitem__(self, i):
+        """Return dataset items from index i."""
+        # batches must be returned in format (B, C, S)
+        # where B = batch_size, C = channels, S = sequence_length
+        from torch import from_numpy
+
+        # select the correct instance
+        n = i // self.single_length
+        # select the correct timepoint starting index based on the selected instance
+        m = i % self.single_length
+
+        hist_end = m + self.seq_len
+        pred_end = m + self.seq_len + self.fh
+
+        historical_y = from_numpy(self.y[n, m:hist_end, :]).float()
+        # historical_y = historical_y.reshape(historical_y.shape[1], -1)
+
+        if self.seq_len < self.moment_seq_len:
+            historical_y = _create_padding(historical_y, self._pad_shape)
+        historical_y = historical_y.float().T
+        future_y = from_numpy(self.y[n, hist_end:pred_end, :]).float().T
+        # future_y = future_y.reshape(future_y.shape[1], -1).T
+        return {
+            "future_y": future_y,
+            "historical_y": historical_y,
+            "input_mask": self.input_mask,
+        }
diff --git a/sktime/libs/README.md b/sktime/libs/README.md
index 890dd6abe..711e07467 100644
--- a/sktime/libs/README.md
+++ b/sktime/libs/README.md
@@ -26,6 +26,9 @@ This folder contains libraries directly distributed with, and maintained by, `sk
 * `granite_ttm` - a package implementing TinyTimeMixer.
   Unofficial fork of package which is not available on pypi.
 
+* `momentfm` - a package implementing the `momentfm` library, unofficial fork
+  maintained since April 2025.
+
 * `time_llm` - partial fork of the `time_llm` package, from [KimMeen/time-LLM](https://github.com/KimMeen/Time-LLM). Unofficial fork of partial code specific to the forecaster. An official package is not available.
 
 * `timemoe` - partial fork of `time-moe` package, from [Time-MoE/Time-MoE](https://github.com/Time-MoE/Time-MoE). Unofficial fork of partial code specific to the forecaster. An official package is not available.
diff --git a/sktime/libs/momentfm/README.md b/sktime/libs/momentfm/README.md
new file mode 100644
index 000000000..e2fdbd830
--- /dev/null
+++ b/sktime/libs/momentfm/README.md
@@ -0,0 +1,30 @@
+# MomentFM: A Family of Open Time-series Foundation Models
+
+Paper: https://arxiv.org/abs/2402.03885
+
+Github Repository: https://github.com/moment-timeseries-foundation-model/moment/tree/main
+
+Authors: Mononito Goswami, Konrad Szafer, Arjun Choudhry, Yifu Cai, Shuo Li, Artur Dubrawski
+
+## Citation
+```
+@inproceedings{goswami2024moment,
+  title={MOMENT: A Family of Open Time-series Foundation Models},
+  author={Mononito Goswami and Konrad Szafer and Arjun Choudhry and Yifu Cai and Shuo Li and Artur Dubrawski},
+  booktitle={International Conference on Machine Learning},
+  year={2024}
+}
+```
+
+## License
+```
+MIT License
+
+Copyright (c) 2024 Auton Lab, Carnegie Mellon University
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+```
diff --git a/sktime/libs/momentfm/__init__.py b/sktime/libs/momentfm/__init__.py
new file mode 100644
index 000000000..b3f958cb7
--- /dev/null
+++ b/sktime/libs/momentfm/__init__.py
@@ -0,0 +1,48 @@
+"""Temporary python implementation of momentfm forecaster.
+
+Unofficial fork of the ``momentfm`` package, maintained in ``sktime``.
+Once a proper package emerges, will be replaced by the package as a dependency
+For more details, please visit https://github.com/moment-timeseries-foundation-model/moment/tree/main
+
+sktime migration: 2024, August
+
+Original authors: Mononito Goswami and Konrad Szafer and Arjun Choudhry and
+Yifu Cai and Shuo Li and Artur Dubrawski
+
+2024 releases subject to following license:
+
+All momentfm code contained is released under the license below.
+
+MIT License
+
+Copyright (c) 2024 Auton Lab, Carnegie Mellon University
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+of the Software, and to permit persons to whom the Software is furnished to do
+so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
+FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
+COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
+IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+Citation:
+
+@inproceedings{goswami2024moment,
+  title={MOMENT: A Family of Open Time-series Foundation Models},
+  author={Mononito Goswami and Konrad Szafer and Arjun Choudhry and Yifu Cai
+  and Shuo Li and Artur Dubrawski},
+  booktitle={International Conference on Machine Learning},
+  year={2024}
+}
+"""
+
+from sktime.libs.momentfm.models.moment import MOMENT, MOMENTPipeline
diff --git a/sktime/libs/momentfm/common.py b/sktime/libs/momentfm/common.py
new file mode 100644
index 000000000..624e523ec
--- /dev/null
+++ b/sktime/libs/momentfm/common.py
@@ -0,0 +1,13 @@
+"""Common file for momentfm."""
+
+from dataclasses import dataclass
+
+
+@dataclass
+class TASKS:
+    """Tasks modules for momentfm models."""
+
+    RECONSTRUCTION: str = "reconstruction"
+    FORECASTING: str = "forecasting"
+    CLASSIFICATION: str = "classification"
+    EMBED: str = "embedding"
diff --git a/sktime/libs/momentfm/data/__init__.py b/sktime/libs/momentfm/data/__init__.py
new file mode 100644
index 000000000..f9a01e698
--- /dev/null
+++ b/sktime/libs/momentfm/data/__init__.py
@@ -0,0 +1,1 @@
+"""Implements the base module for momentfm data types."""
diff --git a/sktime/libs/momentfm/data/base.py b/sktime/libs/momentfm/data/base.py
new file mode 100644
index 000000000..a80e4df47
--- /dev/null
+++ b/sktime/libs/momentfm/data/base.py
@@ -0,0 +1,21 @@
+"""Base module for momentfm TimeSeriesOutputs."""
+
+from dataclasses import dataclass
+
+import numpy.typing as npt
+
+
+@dataclass
+class TimeseriesOutputs:
+    """TimeseriesOutput class."""
+
+    forecast: npt.NDArray = None
+    anomaly_scores: npt.NDArray = None
+    logits: npt.NDArray = None
+    labels: int = None
+    input_mask: npt.NDArray = None
+    pretrain_mask: npt.NDArray = None
+    reconstruction: npt.NDArray = None
+    embeddings: npt.NDArray = None
+    metadata: dict = None
+    illegal_output: bool = False
diff --git a/sktime/libs/momentfm/models/__init__.py b/sktime/libs/momentfm/models/__init__.py
new file mode 100644
index 000000000..00e79a581
--- /dev/null
+++ b/sktime/libs/momentfm/models/__init__.py
@@ -0,0 +1,1 @@
+"""Implement the base model for momentfm."""
diff --git a/sktime/libs/momentfm/models/layers/__init__.py b/sktime/libs/momentfm/models/layers/__init__.py
new file mode 100644
index 000000000..e8a1056a9
--- /dev/null
+++ b/sktime/libs/momentfm/models/layers/__init__.py
@@ -0,0 +1,1 @@
+"""Implement the layers for the momentfm model."""
diff --git a/sktime/libs/momentfm/models/layers/embed.py b/sktime/libs/momentfm/models/layers/embed.py
new file mode 100644
index 000000000..07e2d4611
--- /dev/null
+++ b/sktime/libs/momentfm/models/layers/embed.py
@@ -0,0 +1,326 @@
+"""Embeddings file for momentfm."""
+
+import math
+import warnings
+
+from skbase.utils.dependencies import _check_soft_dependencies
+
+from sktime.libs.momentfm.utils.masking import Masking
+
+if _check_soft_dependencies(["torch"], severity="none"):
+    import torch
+    import torch.nn as nn
+
+    class PositionalEmbedding(nn.Module):
+        """Positional Embedding."""
+
+        def __init__(self, d_model, max_len=5000, model_name="MOMENT"):
+            super().__init__()
+            self.model_name = model_name
+
+            # Compute the positional encodings once in log space.
+            pe = torch.zeros(max_len, d_model).float()
+            pe.require_grad = False
+
+            position = torch.arange(0, max_len).float().unsqueeze(1)
+            div_term = (
+                torch.arange(0, d_model, 2).float() * -(math.log(10000.0) / d_model)
+            ).exp()
+
+            pe[:, 0::2] = torch.sin(position * div_term)
+            pe[:, 1::2] = torch.cos(position * div_term)
+
+            pe = pe.unsqueeze(0)
+            self.register_buffer("pe", pe)
+
+        def forward(self, x):
+            """Forward function."""
+            if (
+                self.model_name == "MOMENT"
+                or self.model_name == "TimesNet"
+                or self.model_name == "GPT4TS"
+            ):
+                return self.pe[:, : x.size(2)]
+            else:
+                return self.pe[:, : x.size(1)]
+
+    class TokenEmbedding(nn.Module):
+        """Token Embeddings."""
+
+        def __init__(self, c_in, d_model):
+            super().__init__()
+            padding = 1 if torch.__version__ >= "1.5.0" else 2
+            self.tokenConv = nn.Conv1d(
+                in_channels=c_in,
+                out_channels=d_model,
+                kernel_size=3,
+                padding=padding,
+                padding_mode="circular",
+                bias=False,
+            )
+            for m in self.modules():
+                if isinstance(m, nn.Conv1d):
+                    nn.init.kaiming_normal_(
+                        m.weight, mode="fan_in", nonlinearity="leaky_relu"
+                    )
+
+        def forward(self, x):
+            """Forward Function."""
+            # x = x.permute(0, 2, 1)
+            x = self.tokenConv(x)
+            x = x.transpose(1, 2)
+            # batch_size x seq_len x d_model
+            return x
+
+    class FixedEmbedding(nn.Module):
+        """Fixed Embeddings."""
+
+        def __init__(self, c_in, d_model):
+            super().__init__()
+
+            w = torch.zeros(c_in, d_model).float()
+            w.require_grad = False
+
+            position = torch.arange(0, c_in).float().unsqueeze(1)
+            div_term = (
+                torch.arange(0, d_model, 2).float() * -(math.log(10000.0) / d_model)
+            ).exp()
+
+            w[:, 0::2] = torch.sin(position * div_term)
+            w[:, 1::2] = torch.cos(position * div_term)
+
+            self.emb = nn.Embedding(c_in, d_model)
+            self.emb.weight = nn.Parameter(w, requires_grad=False)
+
+        def forward(self, x):
+            """Forward Function."""
+            return self.emb(x).detach()
+
+    class TemporalEmbedding(nn.Module):
+        """Temporal Embedding."""
+
+        def __init__(self, d_model, embed_type="fixed", freq="h"):
+            super().__init__()
+
+            minute_size = 4
+            hour_size = 24
+            weekday_size = 7
+            day_size = 32
+            month_size = 13
+
+            Embed = FixedEmbedding if embed_type == "fixed" else nn.Embedding
+            if freq == "t":
+                self.minute_embed = Embed(minute_size, d_model)
+            self.hour_embed = Embed(hour_size, d_model)
+            self.weekday_embed = Embed(weekday_size, d_model)
+            self.day_embed = Embed(day_size, d_model)
+            self.month_embed = Embed(month_size, d_model)
+
+        def forward(self, x):
+            """Forward Function."""
+            x = x.long()
+            minute_x = (
+                self.minute_embed(x[:, :, 4]) if hasattr(self, "minute_embed") else 0.0
+            )
+            hour_x = self.hour_embed(x[:, :, 3])
+            weekday_x = self.weekday_embed(x[:, :, 2])
+            day_x = self.day_embed(x[:, :, 1])
+            month_x = self.month_embed(x[:, :, 0])
+
+            return hour_x + weekday_x + day_x + month_x + minute_x
+
+    class TimeFeatureEmbedding(nn.Module):
+        """Time Feature Embedding."""
+
+        def __init__(self, d_model, embed_type="timeF", freq="h"):
+            super().__init__()
+
+            freq_map = {"h": 4, "t": 5, "s": 6, "m": 1, "a": 1, "w": 2, "d": 3, "b": 3}
+            d_inp = freq_map[freq]
+            self.embed = nn.Linear(d_inp, d_model, bias=False)
+
+        def forward(self, x):
+            """Foward Function."""
+            return self.embed(x)
+
+    class DataEmbedding(nn.Module):
+        """Data Embeddings."""
+
+        def __init__(
+            self, c_in, d_model, model_name, embed_type="fixed", freq="h", dropout=0.1
+        ):
+            super().__init__()
+            self.value_embedding = TokenEmbedding(c_in=c_in, d_model=d_model)
+            self.position_embedding = PositionalEmbedding(
+                d_model=d_model, model_name=model_name
+            )
+            self.temporal_embedding = (
+                TemporalEmbedding(d_model=d_model, embed_type=embed_type, freq=freq)
+                if embed_type != "timeF"
+                else TimeFeatureEmbedding(
+                    d_model=d_model, embed_type=embed_type, freq=freq
+                )
+            )
+            self.dropout = nn.Dropout(p=dropout)
+
+        def forward(self, x, x_mark=None):
+            """Foward Function."""
+            if x_mark is None:
+                x = self.value_embedding(x) + self.position_embedding(x)
+            else:
+                x = (
+                    self.value_embedding(x)
+                    + self.temporal_embedding(x_mark)
+                    + self.position_embedding(x)
+                )
+            return self.dropout(x)
+
+    class DataEmbedding_wo_pos(nn.Module):
+        """Data Embedding without positioning."""
+
+        def __init__(self, c_in, d_model, embed_type="fixed", freq="h", dropout=0.1):
+            super().__init__()
+
+            self.value_embedding = TokenEmbedding(c_in=c_in, d_model=d_model)
+            self.position_embedding = PositionalEmbedding(d_model=d_model)
+            self.temporal_embedding = (
+                TemporalEmbedding(d_model=d_model, embed_type=embed_type, freq=freq)
+                if embed_type != "timeF"
+                else TimeFeatureEmbedding(
+                    d_model=d_model, embed_type=embed_type, freq=freq
+                )
+            )
+            self.dropout = nn.Dropout(p=dropout)
+
+        def forward(self, x, x_mark):
+            """Foward function."""
+            if x_mark is None:
+                x = self.value_embedding(x)
+            else:
+                x = self.value_embedding(x) + self.temporal_embedding(x_mark)
+            return self.dropout(x)
+
+    class PatchEmbedding(nn.Module):
+        """Patch Embedding."""
+
+        def __init__(
+            self,
+            d_model: int = 768,
+            seq_len: int = 512,
+            patch_len: int = 8,
+            stride: int = 8,
+            dropout: int = 0.1,
+            add_positional_embedding: bool = False,
+            value_embedding_bias: bool = False,
+            orth_gain: float = 1.41,
+        ):
+            super().__init__()
+            self.patch_len = patch_len
+            self.seq_len = seq_len
+            self.stride = stride
+            self.d_model = d_model
+            self.add_positional_embedding = add_positional_embedding
+
+            self.value_embedding = nn.Linear(
+                patch_len, d_model, bias=value_embedding_bias
+            )
+            self.mask_embedding = nn.Parameter(torch.zeros(d_model))
+
+            if orth_gain is not None:
+                torch.nn.init.orthogonal_(self.value_embedding.weight, gain=orth_gain)
+                if value_embedding_bias:
+                    self.value_embedding.bias.data.zero_()
+                # torch.nn.init.orthogonal_(self.mask_embedding, gain=orth_gain) # Fails
+
+            # Positional embedding
+            if self.add_positional_embedding:
+                self.position_embedding = PositionalEmbedding(d_model)
+
+            # Residual dropout
+            self.dropout = nn.Dropout(dropout)
+
+        def forward(self, x: torch.Tensor, mask: torch.Tensor = None) -> torch.Tensor:
+            """Foward Function."""
+            mask = Masking.convert_seq_to_patch_view(
+                mask, patch_len=self.patch_len
+            ).unsqueeze(-1)
+            # mask : [batch_size x n_patches x 1]
+            n_channels = x.shape[1]
+            mask = (
+                mask.repeat_interleave(self.d_model, dim=-1)
+                .unsqueeze(1)
+                .repeat(1, n_channels, 1, 1)
+            )
+            # mask : [batch_size x n_channels x n_patches x d_model]
+
+            # Input encoding
+            x = mask * self.value_embedding(x) + (1 - mask) * self.mask_embedding
+            if self.add_positional_embedding:
+                x = x + self.position_embedding(x)
+
+            return self.dropout(x)
+
+    class Patching(nn.Module):
+        """Patching."""
+
+        def __init__(self, patch_len: int, stride: int):
+            super().__init__()
+            self.patch_len = patch_len
+            self.stride = stride
+            if self.stride != self.patch_len:
+                warnings.warn(
+                    "Stride and patch length are not equal. "
+                    "This may lead to unexpected behavior."
+                )
+
+        def forward(self, x):
+            """Foward Function."""
+            x = x.unfold(dimension=-1, size=self.patch_len, step=self.stride)
+            # x : [batch_size x n_channels x num_patch x patch_len]
+            return x
+else:
+
+    class PositionalEmbedding:
+        """Dummy class if torch is unavailable."""
+
+        pass
+
+    class TokenEmbedding:
+        """Dummy class if torch is unavailable."""
+
+        pass
+
+    class FixedEmbedding:
+        """Dummy class if torch is unavailable."""
+
+        pass
+
+    class TimeFeatureEmbedding:
+        """Dummy class if torch is unavailable."""
+
+        pass
+
+    class TemporalEmbedding:
+        """Dummy class if torch is unavailable."""
+
+        pass
+
+    class DataEmbedding_wo_pos:
+        """Dummy class if torch is unavailable."""
+
+        pass
+
+    class DataEmbedding:
+        """Dummy class if torch is unavailable."""
+
+        pass
+
+    class PatchEmbedding:
+        """Dummy class if torch is unavailable."""
+
+        pass
+
+    class Patching:
+        """Dummy class if torch is unavailable."""
+
+        pass
diff --git a/sktime/libs/momentfm/models/layers/revin.py b/sktime/libs/momentfm/models/layers/revin.py
new file mode 100644
index 000000000..8e94707e6
--- /dev/null
+++ b/sktime/libs/momentfm/models/layers/revin.py
@@ -0,0 +1,103 @@
+"""Layer revin file."""
+
+from skbase.utils.dependencies import _check_soft_dependencies
+
+if _check_soft_dependencies(["torch"], severity="none"):
+    import torch
+    from torch import nn
+
+    def nanvar(tensor, dim=None, keepdim=False):
+        """Calculate the variance of all non-NaN elements."""
+        tensor_mean = tensor.nanmean(dim=dim, keepdim=True)
+        output = (tensor - tensor_mean).square().nanmean(dim=dim, keepdim=keepdim)
+        return output
+
+    def nanstd(tensor, dim=None, keepdim=False):
+        """Calculate the std of all non-NaN elements."""
+        output = nanvar(tensor, dim=dim, keepdim=keepdim)
+        output = output.sqrt()
+        return output
+
+    class RevIN(nn.Module):
+        """RevIN class."""
+
+        def __init__(self, num_features: int, eps: float = 1e-5, affine: bool = False):
+            """Initialize RevIN layer.
+
+            :param num_features: the number of features or channels
+            :param eps: a value added for numerical stability
+            :param affine: if True, RevIN has learnable affine parameters
+            """
+            super().__init__()
+            self.num_features = num_features
+            self.eps = eps
+            self.affine = affine
+
+            if self.affine:
+                self._init_params()
+
+        def forward(
+            self, x: torch.Tensor, mode: str = "norm", mask: torch.Tensor = None
+        ):
+            """Forward Function.
+
+            :param x: input tensor of shape (batch_size, n_channels, seq_len)
+            :param mode: 'norm' or 'denorm'
+            :param mask: input mask of shape (batch_size, seq_len)
+            :return: RevIN transformed tensor
+            """
+            if mode == "norm":
+                self._get_statistics(x, mask=mask)
+                x = self._normalize(x)
+            elif mode == "denorm":
+                x = self._denormalize(x)
+            else:
+                raise NotImplementedError
+            return x
+
+        def _init_params(self):
+            # initialize RevIN params: (C,)
+            self.affine_weight = nn.Parameter(torch.ones(1, self.num_features, 1))
+            self.affine_bias = nn.Parameter(torch.zeros(1, self.num_features, 1))
+
+        def _get_statistics(self, x, mask=None):
+            """Get statistics.
+
+            x    : batch_size x n_channels x seq_len
+            mask : batch_size x seq_len
+            """
+            if mask is None:
+                mask = torch.ones((x.shape[0], x.shape[-1]))
+            n_channels = x.shape[1]
+            mask = mask.unsqueeze(1).repeat(1, n_channels, 1).bool()
+            # Set masked positions to NaN, and unmasked positions are taken from x
+            masked_x = torch.where(mask, x, torch.nan)
+            self.mean = torch.nanmean(masked_x, dim=-1, keepdim=True).detach()
+            self.stdev = nanstd(masked_x, dim=-1, keepdim=True).detach() + self.eps
+            # self.stdev = torch.sqrt(
+            #     torch.var(masked_x, dim=-1, keepdim=True) + self.eps)
+            # .get_data().detach()
+            # NOTE: By default not bessel correction
+
+        def _normalize(self, x):
+            x = x - self.mean
+            x = x / self.stdev
+
+            if self.affine:
+                x = x * self.affine_weight
+                x = x + self.affine_bias
+            return x
+
+        def _denormalize(self, x):
+            if self.affine:
+                x = x - self.affine_bias
+                x = x / (self.affine_weight + self.eps * self.eps)
+            x = x * self.stdev
+            x = x + self.mean
+            return x
+else:
+
+    class RevIN:
+        """Dummy class if torch is unavailable."""
+
+    pass
diff --git a/sktime/libs/momentfm/models/moment.py b/sktime/libs/momentfm/models/moment.py
new file mode 100644
index 000000000..9b9bbf66d
--- /dev/null
+++ b/sktime/libs/momentfm/models/moment.py
@@ -0,0 +1,695 @@
+"""Moment model and heads file."""
+
+import logging
+import warnings
+from argparse import Namespace
+from copy import deepcopy
+from math import ceil
+
+from skbase.utils.dependencies import _check_soft_dependencies
+
+from sktime.libs.momentfm.common import TASKS
+from sktime.libs.momentfm.data.base import TimeseriesOutputs
+from sktime.libs.momentfm.models.layers.embed import PatchEmbedding, Patching
+from sktime.libs.momentfm.models.layers.revin import RevIN
+from sktime.libs.momentfm.utils.masking import Masking
+from sktime.libs.momentfm.utils.utils import (
+    NamespaceWithDefaults,
+    get_anomaly_criterion,
+    get_huggingface_model_dimensions,
+)
+
+SUPPORTED_HUGGINGFACE_MODELS = [
+    "google/flan-t5-small",
+    "google/flan-t5-base",
+    "google/flan-t5-large",
+    "google/flan-t5-xl",
+    "google/flan-t5-xxl",
+]
+
+if _check_soft_dependencies(
+    ["torch", "huggingface-hub", "transformers"], severity="none"
+):
+    import torch
+    from huggingface_hub import PyTorchModelHubMixin
+    from torch import nn
+    from transformers import T5Config, T5EncoderModel, T5Model
+
+    class PretrainHead(nn.Module):
+        """Pretrained Head."""
+
+        def __init__(
+            self,
+            d_model: int = 768,
+            patch_len: int = 8,
+            head_dropout: float = 0.1,
+            orth_gain: float = 1.41,
+        ):
+            super().__init__()
+            self.dropout = nn.Dropout(head_dropout)
+            self.linear = nn.Linear(d_model, patch_len)
+
+            if orth_gain is not None:
+                torch.nn.init.orthogonal_(self.linear.weight, gain=orth_gain)
+                self.linear.bias.data.zero_()
+
+        def forward(self, x):
+            """Foward Function."""
+            x = self.linear(self.dropout(x))
+            x = x.flatten(start_dim=2, end_dim=3)
+            return x
+
+    class ClassificationHead(nn.Module):
+        """Classification Head."""
+
+        def __init__(
+            self,
+            n_channels: int = 1,
+            d_model: int = 768,
+            n_classes: int = 2,
+            head_dropout: int = 0.1,
+            reduction: str = "concat",
+        ):
+            super().__init__()
+            self.dropout = nn.Dropout(head_dropout)
+            if reduction == "mean":
+                self.linear = nn.Linear(d_model, n_classes)
+            elif reduction == "concat":
+                self.linear = nn.Linear(n_channels * d_model, n_classes)
+            else:
+                raise ValueError(
+                    f"Reduction method {reduction} not implemented. "
+                    f"Only 'mean' and 'concat' are supported."
+                )
+
+        def forward(self, x, input_mask: torch.Tensor = None):
+            """Foward Function."""
+            x = torch.mean(x, dim=1)
+            x = self.dropout(x)
+            y = self.linear(x)
+            return y
+
+    class ForecastingHead(nn.Module):
+        """Forecasting Head."""
+
+        def __init__(
+            self,
+            head_nf: int = 768 * 64,
+            forecast_horizon: int = 96,
+            head_dropout: int = 0,
+        ):
+            super().__init__()
+            self.flatten = nn.Flatten(start_dim=-2)
+            self.dropout = nn.Dropout(head_dropout)
+            self.linear = nn.Linear(head_nf, forecast_horizon)
+
+        def forward(self, x, input_mask: torch.Tensor = None):
+            """Foward Function."""
+            x = self.flatten(x)
+            x = self.linear(x)
+            x = self.dropout(x)
+            return x
+
+    class MOMENT(nn.Module):
+        """Moment Class."""
+
+        def __init__(self, config, **kwargs: dict):
+            super().__init__()
+            config = self._update_inputs(config, **kwargs)
+            config = self._validate_inputs(config)
+            self.config = config
+            self.task_name = config.task_name
+            self.seq_len = config.seq_len
+            self.patch_len = config.patch_len
+
+            self.normalizer = RevIN(
+                num_features=1, affine=config.getattr("revin_affine", False)
+            )
+            self.tokenizer = Patching(
+                patch_len=config.patch_len, stride=config.patch_stride_len
+            )
+            self.patch_embedding = PatchEmbedding(
+                d_model=config.d_model,
+                seq_len=config.seq_len,
+                patch_len=config.patch_len,
+                stride=config.patch_stride_len,
+                dropout=config.getattr("dropout", 0.1),
+                add_positional_embedding=config.getattr(
+                    "add_positional_embedding", True
+                ),
+                value_embedding_bias=config.getattr("value_embedding_bias", False),
+                orth_gain=config.getattr("orth_gain", 1.41),
+            )
+            self.mask_generator = Masking(mask_ratio=config.getattr("mask_ratio", 0.0))
+            self.encoder = self._get_transformer_backbone(config)
+            self.head = self._get_head(self.task_name)
+
+            # Frozen parameters
+            self.freeze_embedder = config.getattr("freeze_embedder", True)
+            self.freeze_encoder = config.getattr("freeze_encoder", True)
+            self.freeze_head = config.getattr("freeze_head", False)
+
+            if self.freeze_embedder:
+                self.patch_embedding = freeze_parameters(self.patch_embedding)
+            if self.freeze_encoder:
+                self.encoder = freeze_parameters(self.encoder)
+            if self.freeze_head:
+                self.head = freeze_parameters(self.head)
+
+        def _update_inputs(self, config, **kwargs: dict) -> NamespaceWithDefaults:
+            """Update Inputs."""
+            if isinstance(config, dict) and "model_kwargs" in kwargs:
+                return NamespaceWithDefaults(**{**config, **kwargs["model_kwargs"]})
+            else:
+                return NamespaceWithDefaults.from_namespace(config)
+
+        def _validate_inputs(
+            self, config: NamespaceWithDefaults
+        ) -> NamespaceWithDefaults:
+            """Validate Inputs."""
+            if (
+                config.d_model is None
+                and config.transformer_backbone in SUPPORTED_HUGGINGFACE_MODELS
+            ):
+                config.d_model = get_huggingface_model_dimensions(
+                    config.transformer_backbone
+                )
+                logging.info(f"Setting d_model to {config.d_model}")
+            elif config.d_model is None:
+                raise ValueError(
+                    "d_model must be specified if transformer backbone "
+                    "unless transformer backbone is a Huggingface model."
+                )
+
+            if config.transformer_type not in [
+                "encoder_only",
+                "decoder_only",
+                "encoder_decoder",
+            ]:
+                raise ValueError(
+                    "transformer_type must be one of "
+                    "['encoder_only', 'decoder_only', 'encoder_decoder']"
+                )
+
+            if config.patch_stride_len != config.patch_len:
+                warnings.warn("Patch stride length is not equal to patch length.")
+            return config
+
+        def _get_head(self, task_name: str) -> nn.Module:
+            if task_name == TASKS.RECONSTRUCTION:
+                return PretrainHead(
+                    self.config.d_model,
+                    self.config.patch_len,
+                    self.config.getattr("dropout", 0.1),
+                    self.config.getattr("orth_gain", 1.41),
+                )
+            elif task_name == TASKS.CLASSIFICATION:
+                return ClassificationHead(
+                    self.config.n_channels,
+                    self.config.d_model,
+                    self.config.num_class,
+                    self.config.getattr("dropout", 0.1),
+                    reduction=self.config.getattr("reduction", "concat"),
+                )
+            elif task_name == TASKS.FORECASTING:
+                num_patches = (
+                    max(self.config.seq_len, self.config.patch_len)
+                    - self.config.patch_len
+                ) // self.config.patch_stride_len + 1
+                self.head_nf = self.config.d_model * num_patches
+                return ForecastingHead(
+                    self.head_nf,
+                    self.config.forecast_horizon,
+                    self.config.getattr("head_dropout", 0.1),
+                )
+            elif task_name == TASKS.EMBED:
+                return nn.Identity()
+            else:
+                raise NotImplementedError(f"Task {task_name} not implemented.")
+
+        def _get_transformer_backbone(self, config) -> nn.Module:
+            if config.getattr("randomly_initialize_backbone", False):
+                model_config = T5Config.from_pretrained(config.transformer_backbone)
+                transformer_backbone = T5Model(model_config)
+                logging.info(
+                    f"Initializing randomly initialized transformer from "
+                    f"{config.transformer_backbone}."
+                )
+            else:
+                transformer_backbone = T5EncoderModel.from_pretrained(
+                    config.transformer_backbone
+                )
+                logging.info(
+                    f"Initializing pre-trained transformer from "
+                    f"{config.transformer_backbone}."
+                )
+
+            transformer_backbone = transformer_backbone.get_encoder()
+
+            if config.getattr("enable_gradient_checkpointing", True):
+                transformer_backbone.gradient_checkpointing_enable()
+                logging.info("Enabling gradient checkpointing.")
+
+            return transformer_backbone
+
+        def __call__(self, *args, **kwargs) -> TimeseriesOutputs:
+            """__call__ function."""
+            return self.forward(*args, **kwargs)
+
+        def embed(
+            self,
+            x_enc: torch.Tensor,
+            input_mask: torch.Tensor = None,
+            reduction: str = "mean",
+            **kwargs,
+        ) -> TimeseriesOutputs:
+            """Embed function."""
+            batch_size, n_channels, seq_len = x_enc.shape
+
+            if input_mask is None:
+                input_mask = torch.ones((batch_size, seq_len)).to(x_enc.device)
+
+            x_enc = self.normalizer(x=x_enc, mask=input_mask, mode="norm")
+            x_enc = torch.nan_to_num(x_enc, nan=0, posinf=0, neginf=0)
+
+            input_mask_patch_view = Masking.convert_seq_to_patch_view(
+                input_mask, self.patch_len
+            )
+
+            x_enc = self.tokenizer(x=x_enc)
+            enc_in = self.patch_embedding(x_enc, mask=input_mask)
+
+            n_patches = enc_in.shape[2]
+            enc_in = enc_in.reshape(
+                (batch_size * n_channels, n_patches, self.config.d_model)
+            )
+
+            patch_view_mask = Masking.convert_seq_to_patch_view(
+                input_mask, self.patch_len
+            )
+            attention_mask = patch_view_mask.repeat_interleave(n_channels, dim=0)
+            outputs = self.encoder(inputs_embeds=enc_in, attention_mask=attention_mask)
+            enc_out = outputs.last_hidden_state
+
+            enc_out = enc_out.reshape((-1, n_channels, n_patches, self.config.d_model))
+            # [batch_size x n_channels x n_patches x d_model]
+
+            if reduction == "mean":
+                enc_out = enc_out.mean(dim=1, keepdim=False)  # Mean across channels
+                # [batch_size x n_patches x d_model]
+                input_mask_patch_view = input_mask_patch_view.unsqueeze(-1).repeat(
+                    1, 1, self.config.d_model
+                )
+                enc_out = (input_mask_patch_view * enc_out).sum(
+                    dim=1
+                ) / input_mask_patch_view.sum(dim=1)
+            else:
+                raise NotImplementedError(
+                    f"Reduction method {reduction} not implemented."
+                )
+
+            return TimeseriesOutputs(
+                embeddings=enc_out, input_mask=input_mask, metadata=reduction
+            )
+
+        def reconstruction(
+            self,
+            x_enc: torch.Tensor,
+            input_mask: torch.Tensor = None,
+            mask: torch.Tensor = None,
+            **kwargs,
+        ) -> TimeseriesOutputs:
+            """Reconstruction Function."""
+            batch_size, n_channels, _ = x_enc.shape
+
+            if mask is None:
+                mask = self.mask_generator.generate_mask(x=x_enc, input_mask=input_mask)
+                mask = mask.to(x_enc.device)  # mask: [batch_size x seq_len]
+
+            x_enc = self.normalizer(x=x_enc, mask=mask * input_mask, mode="norm")
+            # Prevent too short time-series from causing NaNs
+            x_enc = torch.nan_to_num(x_enc, nan=0, posinf=0, neginf=0)
+
+            x_enc = self.tokenizer(x=x_enc)
+            enc_in = self.patch_embedding(x_enc, mask=mask)
+
+            n_patches = enc_in.shape[2]
+            enc_in = enc_in.reshape(
+                (batch_size * n_channels, n_patches, self.config.d_model)
+            )
+
+            patch_view_mask = Masking.convert_seq_to_patch_view(
+                input_mask, self.patch_len
+            )
+            attention_mask = patch_view_mask.repeat_interleave(n_channels, dim=0)
+            if self.config.transformer_type == "encoder_decoder":
+                outputs = self.encoder(
+                    inputs_embeds=enc_in,
+                    decoder_inputs_embeds=enc_in,
+                    attention_mask=attention_mask,
+                )
+            else:
+                outputs = self.encoder(
+                    inputs_embeds=enc_in, attention_mask=attention_mask
+                )
+            enc_out = outputs.last_hidden_state
+
+            enc_out = enc_out.reshape((-1, n_channels, n_patches, self.config.d_model))
+
+            dec_out = self.head(enc_out)  # [batch_size x n_channels x seq_len]
+            dec_out = self.normalizer(x=dec_out, mode="denorm")
+
+            if self.config.getattr("debug", False):
+                illegal_output = self._check_model_weights_for_illegal_values()
+            else:
+                illegal_output = None
+
+            return TimeseriesOutputs(
+                input_mask=input_mask,
+                reconstruction=dec_out,
+                pretrain_mask=mask,
+                illegal_output=illegal_output,
+            )
+
+        def reconstruct(
+            self,
+            x_enc: torch.Tensor,
+            input_mask: torch.Tensor = None,
+            mask: torch.Tensor = None,
+            **kwargs,
+        ) -> TimeseriesOutputs:
+            """Reconstruct Function."""
+            if mask is None:
+                mask = torch.ones_like(input_mask)
+
+            batch_size, n_channels, _ = x_enc.shape
+            x_enc = self.normalizer(x=x_enc, mask=mask * input_mask, mode="norm")
+
+            x_enc = self.tokenizer(x=x_enc)
+            enc_in = self.patch_embedding(x_enc, mask=mask)
+
+            n_patches = enc_in.shape[2]
+            enc_in = enc_in.reshape(
+                (batch_size * n_channels, n_patches, self.config.d_model)
+            )
+            # [batch_size * n_channels x n_patches x d_model]
+
+            patch_view_mask = Masking.convert_seq_to_patch_view(
+                input_mask, self.patch_len
+            )
+            attention_mask = patch_view_mask.repeat_interleave(n_channels, dim=0).to(
+                x_enc.device
+            )
+
+            n_tokens = 0
+            if "prompt_embeds" in kwargs:
+                prompt_embeds = kwargs["prompt_embeds"].to(x_enc.device)
+
+                if isinstance(prompt_embeds, nn.Embedding):
+                    prompt_embeds = prompt_embeds.weight.data.unsqueeze(0)
+
+                n_tokens = prompt_embeds.shape[1]
+
+                enc_in = self._cat_learned_embedding_to_input(prompt_embeds, enc_in)
+                attention_mask = self._extend_attention_mask(attention_mask, n_tokens)
+
+            if self.config.transformer_type == "encoder_decoder":
+                outputs = self.encoder(
+                    inputs_embeds=enc_in,
+                    decoder_inputs_embeds=enc_in,
+                    attention_mask=attention_mask,
+                )
+            else:
+                outputs = self.encoder(
+                    inputs_embeds=enc_in, attention_mask=attention_mask
+                )
+            enc_out = outputs.last_hidden_state
+            enc_out = enc_out[:, n_tokens:, :]
+
+            enc_out = enc_out.reshape((-1, n_channels, n_patches, self.config.d_model))
+            # [batch_size x n_channels x n_patches x d_model]
+
+            dec_out = self.head(enc_out)  # [batch_size x n_channels x seq_len]
+            dec_out = self.normalizer(x=dec_out, mode="denorm")
+
+            return TimeseriesOutputs(input_mask=input_mask, reconstruction=dec_out)
+
+        def detect_anomalies(
+            self,
+            x_enc: torch.Tensor,
+            input_mask: torch.Tensor = None,
+            anomaly_criterion: str = "mse",
+            **kwargs,
+        ) -> TimeseriesOutputs:
+            """Detect Anomalies Function."""
+            outputs = self.reconstruct(x_enc=x_enc, input_mask=input_mask)
+            self.anomaly_criterion = get_anomaly_criterion(anomaly_criterion)
+
+            anomaly_scores = self.anomaly_criterion(x_enc, outputs.reconstruction)
+
+            return TimeseriesOutputs(
+                input_mask=input_mask,
+                reconstruction=outputs.reconstruction,
+                anomaly_scores=anomaly_scores,
+                metadata={"anomaly_criterion": anomaly_criterion},
+            )
+
+        def forecast(
+            self, x_enc: torch.Tensor, input_mask: torch.Tensor = None, **kwargs
+        ) -> TimeseriesOutputs:
+            """Forecast Function."""
+            batch_size, n_channels, seq_len = x_enc.shape
+
+            x_enc = self.normalizer(x=x_enc, mask=input_mask, mode="norm")
+            x_enc = torch.nan_to_num(x_enc, nan=0, posinf=0, neginf=0)
+
+            x_enc = self.tokenizer(x=x_enc)
+            enc_in = self.patch_embedding(x_enc, mask=torch.ones_like(input_mask))
+
+            n_patches = enc_in.shape[2]
+            enc_in = enc_in.reshape(
+                (batch_size * n_channels, n_patches, self.config.d_model)
+            )
+
+            patch_view_mask = Masking.convert_seq_to_patch_view(
+                input_mask, self.patch_len
+            )
+            attention_mask = patch_view_mask.repeat_interleave(n_channels, dim=0)
+            outputs = self.encoder(inputs_embeds=enc_in, attention_mask=attention_mask)
+            enc_out = outputs.last_hidden_state
+            enc_out = enc_out.reshape((-1, n_channels, n_patches, self.config.d_model))
+            # [batch_size x n_channels x n_patches x d_model]
+
+            dec_out = self.head(enc_out)  # [batch_size x n_channels x forecast_horizon]
+            dec_out = self.normalizer(x=dec_out, mode="denorm")
+
+            return TimeseriesOutputs(input_mask=input_mask, forecast=dec_out)
+
+        def short_forecast(
+            self,
+            x_enc: torch.Tensor,
+            input_mask: torch.Tensor = None,
+            forecast_horizon: int = 1,
+            **kwargs,
+        ) -> TimeseriesOutputs:
+            """Short_forecast Function."""
+            batch_size, n_channels, seq_len = x_enc.shape
+            num_masked_patches = ceil(forecast_horizon / self.patch_len)
+            num_masked_timesteps = num_masked_patches * self.patch_len
+
+            x_enc = self.normalizer(x=x_enc, mask=input_mask, mode="norm")
+            x_enc = torch.nan_to_num(x_enc, nan=0, posinf=0, neginf=0)
+
+            # Shift the time-series and mask the last few timesteps for forecasting
+            x_enc = torch.roll(x_enc, shifts=-num_masked_timesteps, dims=2)
+            input_mask = torch.roll(input_mask, shifts=-num_masked_timesteps, dims=1)
+
+            # Attending to mask tokens
+            input_mask[:, -num_masked_timesteps:] = 1
+            mask = torch.ones_like(input_mask)
+            mask[:, -num_masked_timesteps:] = 0
+
+            x_enc = self.tokenizer(x=x_enc)
+            enc_in = self.patch_embedding(x_enc, mask=mask)
+
+            n_patches = enc_in.shape[2]
+            enc_in = enc_in.reshape(
+                (batch_size * n_channels, n_patches, self.config.d_model)
+            )
+            # [batch_size * n_channels x n_patches x d_model]
+
+            patch_view_mask = Masking.convert_seq_to_patch_view(
+                input_mask, self.patch_len
+            )
+            attention_mask = patch_view_mask.repeat_interleave(n_channels, dim=0)
+            outputs = self.encoder(inputs_embeds=enc_in, attention_mask=attention_mask)
+            enc_out = outputs.last_hidden_state
+            enc_out = enc_out.reshape((-1, n_channels, n_patches, self.config.d_model))
+
+            dec_out = self.head(enc_out)  # [batch_size x n_channels x seq_len]
+
+            end = -num_masked_timesteps + forecast_horizon
+            end = None if end == 0 else end
+
+            dec_out = self.normalizer(x=dec_out, mode="denorm")
+            forecast = dec_out[:, :, -num_masked_timesteps:end]
+
+            return TimeseriesOutputs(
+                input_mask=input_mask,
+                reconstruction=dec_out,
+                forecast=forecast,
+                metadata={"forecast_horizon": forecast_horizon},
+            )
+
+        def classify(
+            self,
+            x_enc: torch.Tensor,
+            input_mask: torch.Tensor = None,
+            reduction: str = "concat",
+            **kwargs,
+        ) -> TimeseriesOutputs:
+            """Classify Function."""
+            batch_size, n_channels, seq_len = x_enc.shape
+
+            if input_mask is None:
+                input_mask = torch.ones((batch_size, seq_len)).to(x_enc.device)
+
+            x_enc = self.normalizer(x=x_enc, mask=input_mask, mode="norm")
+            x_enc = torch.nan_to_num(x_enc, nan=0, posinf=0, neginf=0)
+
+            # input_mask_patch_view = Masking.convert_seq_to_patch_view(
+            #     input_mask, self.patch_len
+            # )
+
+            x_enc = self.tokenizer(x=x_enc)
+            enc_in = self.patch_embedding(x_enc, mask=input_mask)
+
+            n_patches = enc_in.shape[2]
+            enc_in = enc_in.reshape(
+                (batch_size * n_channels, n_patches, self.config.d_model)
+            )
+
+            patch_view_mask = Masking.convert_seq_to_patch_view(
+                input_mask, self.patch_len
+            )
+            attention_mask = patch_view_mask.repeat_interleave(n_channels, dim=0)
+            outputs = self.encoder(inputs_embeds=enc_in, attention_mask=attention_mask)
+            enc_out = outputs.last_hidden_state
+
+            enc_out = enc_out.reshape((-1, n_channels, n_patches, self.config.d_model))
+            # [batch_size x n_channels x n_patches x d_model]
+
+            # Mean across channels
+            if reduction == "mean":
+                # [batch_size x n_patches x d_model]
+                enc_out = enc_out.mean(dim=1, keepdim=False)
+            # Concatenate across channels
+            elif reduction == "concat":
+                # [batch_size x n_patches x d_model * n_channels]
+                enc_out = enc_out.permute(0, 2, 3, 1).reshape(
+                    batch_size, n_patches, self.config.d_model * n_channels
+                )
+
+            else:
+                raise NotImplementedError(
+                    f"Reduction method {reduction} not implemented."
+                )
+
+            logits = self.head(enc_out, input_mask=input_mask)
+
+            return TimeseriesOutputs(
+                embeddings=enc_out, logits=logits, metadata=reduction
+            )
+
+        def forward(
+            self,
+            x_enc: torch.Tensor,
+            mask: torch.Tensor = None,
+            input_mask: torch.Tensor = None,
+            **kwargs,
+        ) -> TimeseriesOutputs:
+            """Foward Function."""
+            if input_mask is None:
+                input_mask = torch.ones_like(x_enc[:, 0, :])
+
+            if self.task_name == TASKS.RECONSTRUCTION:
+                return self.reconstruction(
+                    x_enc=x_enc, mask=mask, input_mask=input_mask, **kwargs
+                )
+            elif self.task_name == TASKS.EMBED:
+                return self.embed(x_enc=x_enc, input_mask=input_mask, **kwargs)
+            elif self.task_name == TASKS.FORECASTING:
+                return self.forecast(x_enc=x_enc, input_mask=input_mask, **kwargs)
+            elif self.task_name == TASKS.CLASSIFICATION:
+                return self.classify(x_enc=x_enc, input_mask=input_mask, **kwargs)
+            else:
+                raise NotImplementedError(f"Task {self.task_name} not implemented.")
+
+    class MOMENTPipeline(MOMENT, PyTorchModelHubMixin):
+        """Moment Pipeline."""
+
+        def __init__(self, config, **kwargs: dict):
+            self._validate_model_kwargs(**kwargs)
+            self.new_task_name = kwargs.get("model_kwargs", {}).pop(
+                "task_name", TASKS.RECONSTRUCTION
+            )
+            super().__init__(config, **kwargs)
+
+        def _validate_model_kwargs(self, **kwargs: dict) -> None:
+            kwargs = deepcopy(kwargs)
+            kwargs.setdefault("model_kwargs", {"task_name": TASKS.RECONSTRUCTION})
+            kwargs["model_kwargs"].setdefault("task_name", TASKS.RECONSTRUCTION)
+            config = Namespace(**kwargs["model_kwargs"])
+
+            if config.task_name == TASKS.FORECASTING:
+                if not hasattr(config, "forecast_horizon"):
+                    raise ValueError(
+                        "forecast_horizon must be specified for long-horizon "
+                        "forecasting."
+                    )
+
+            if config.task_name == TASKS.CLASSIFICATION:
+                if not hasattr(config, "n_channels"):
+                    raise ValueError("n_channels must be specified for classification.")
+                if not hasattr(config, "num_class"):
+                    raise ValueError("num_class must be specified for classification.")
+
+        def init(self) -> None:
+            """Init function."""
+            if self.new_task_name != TASKS.RECONSTRUCTION:
+                self.task_name = self.new_task_name
+                self.head = self._get_head(self.new_task_name)
+
+    def freeze_parameters(model):
+        """Freeze parameters of the model."""
+        # Freeze the parameters
+        for name, param in model.named_parameters():
+            param.requires_grad = False
+
+        return model
+else:
+
+    class PretrainHead:
+        """Dummy class if soft dependencies are not available."""
+
+        pass
+
+    class ClassificationHead:
+        """Dummy class if soft dependencies are not available."""
+
+        pass
+
+    class ForecastingHead:
+        """Dummy class if soft dependencies are not available."""
+
+        pass
+
+    class MOMENT:
+        """Dummy class if soft dependencies are not available."""
+
+        pass
+
+    class MOMENTPipeline:
+        """Dummy class if soft dependencies are not available."""
+
+        pass
diff --git a/sktime/libs/momentfm/utils/__init__.py b/sktime/libs/momentfm/utils/__init__.py
new file mode 100644
index 000000000..139c40ed0
--- /dev/null
+++ b/sktime/libs/momentfm/utils/__init__.py
@@ -0,0 +1,1 @@
+"""Implement the utilities for momentfm."""
diff --git a/sktime/libs/momentfm/utils/anomaly_detection_metrics.py b/sktime/libs/momentfm/utils/anomaly_detection_metrics.py
new file mode 100644
index 000000000..919cb4a99
--- /dev/null
+++ b/sktime/libs/momentfm/utils/anomaly_detection_metrics.py
@@ -0,0 +1,85 @@
+"""Anomaly detection metrics."""
+
+import numpy as np
+
+
+def f1_score(predict, actual):
+    """Calculate f1 score."""
+    TP = np.sum(predict * actual)
+    # TN = np.sum((1 - predict) * (1 - actual))
+    FP = np.sum(predict * (1 - actual))
+    FN = np.sum((1 - predict) * actual)
+    precision = TP / (TP + FP + 0.00001)
+    recall = TP / (TP + FN + 0.00001)
+    f1 = 2 * precision * recall / (precision + recall + 0.00001)
+    return f1
+
+
+def adjust_predicts(score, label, threshold=None, pred=None, calc_latency=False):
+    """Calculate Adjusted predict labels.
+
+    Calculate adjusted predict labels using given `score`,
+    `threshold` (or given `pred`) and `label`.
+    Args:
+        score (np.ndarray): The anomaly score
+        label (np.ndarray): The ground-truth label
+        threshold (float): The threshold of anomaly score.
+            A point is labeled as "anomaly" if its score is lower than the
+            threshold.
+        pred (np.ndarray or None): if not None, adjust `pred` and ignore
+        `score` and `threshold`,
+        calc_latency (bool):
+    Returns:
+        np.ndarray: predict labels
+    """
+    if len(score) != len(label):
+        raise ValueError("score and label must have the same length")
+    score = np.asarray(score)
+    label = np.asarray(label)
+    latency = 0
+    if pred is None:
+        predict = score < threshold
+    else:
+        predict = pred
+    actual = label > 0.1
+    anomaly_state = False
+    anomaly_count = 0
+    for i in range(len(score)):
+        if actual[i] and predict[i] and not anomaly_state:
+            anomaly_state = True
+            anomaly_count += 1
+            for j in range(i, 0, -1):
+                if not actual[j]:
+                    break
+                else:
+                    if not predict[j]:
+                        predict[j] = True
+                        latency += 1
+        elif not actual[i]:
+            anomaly_state = False
+        if anomaly_state:
+            predict[i] = True
+    if calc_latency:
+        return predict, latency / (anomaly_count + 1e-4)
+    else:
+        return predict
+
+
+def adjbestf1(y_true: np.array, y_scores: np.array, n_splits: int = 100):
+    """Calculate adjusted best f1."""
+    thresholds = np.linspace(y_scores.min(), y_scores.max(), n_splits)
+    adjusted_f1 = np.zeros(thresholds.shape)
+
+    for i, threshold in enumerate(thresholds):
+        y_pred = y_scores >= threshold
+        y_pred = adjust_predicts(
+            score=y_scores,
+            label=(y_true > 0),
+            pred=y_pred,
+            threshold=None,
+            calc_latency=False,
+        )
+        adjusted_f1[i] = f1_score(y_pred, y_true)
+
+    best_adjusted_f1 = np.max(adjusted_f1)
+    return best_adjusted_f1
diff --git a/sktime/libs/momentfm/utils/data.py b/sktime/libs/momentfm/utils/data.py
new file mode 100644
index 000000000..2bb6c9bef
--- /dev/null
+++ b/sktime/libs/momentfm/utils/data.py
@@ -0,0 +1,203 @@
+"""Data Utility file in momentfm."""
+
+import numpy as np
+
+
+def load_from_tsfile(
+    full_file_path_and_name,
+    replace_missing_vals_with="NaN",
+    return_meta_data=False,
+    return_type="auto",
+):
+    """Load time series .ts file into X and (optionally) y.
+
+    Parameters
+    ----------
+    full_file_path_and_name : string
+        full path of the file to load, .ts extension is assumed.
+    replace_missing_vals_with : string, default="NaN"
+        issing values in the file are replaces with this value
+    return_meta_data : boolean, default=False
+        return a dictionary with the meta data loaded from the file
+    return_type : string, default = "auto"
+        data type to convert to.
+        If "auto", returns numpy3D for equal length and list of numpy2D for unequal.
+        If "numpy2D", will squash a univariate equal length into a numpy2D (n_cases,
+        n_timepoints). Other options are available but not supported medium term.
+
+    Returns
+    -------
+    data: Union[np.ndarray,list]
+        time series data, np.ndarray (n_cases, n_channels, series_length) if equal
+        length time series, list of [n_cases] np.ndarray (n_channels, n_timepoints)
+        if unequal length series.
+    y : target variable, np.ndarray of string or int
+    meta_data : dict (optional).
+        dictionary of characteristics, with keys
+        "problemname" (string), booleans: "timestamps", "missing", "univariate",
+        "equallength", "classlabel", "targetlabel" and "class_values": [],
+
+    Raises
+    ------
+    IOError if the load fails.
+    """
+    # Check file ends in .ts, if not, insert
+    if not full_file_path_and_name.endswith(".ts"):
+        full_file_path_and_name = full_file_path_and_name + ".ts"
+    # Open file
+    with open(full_file_path_and_name, encoding="utf-8") as file:
+        # Read in headers
+        meta_data = _load_header_info(file)
+        # load into list of numpy
+        data, y, meta_data = _load_data(file, meta_data)
+
+    # if equal load to 3D numpy
+    if meta_data["equallength"]:
+        data = np.array(data)
+        if return_type == "numpy2D" and meta_data["univariate"]:
+            data = data.squeeze()
+    # If regression problem, convert y to float
+    if meta_data["targetlabel"]:
+        y = y.astype(float)
+
+    return (data, y, meta_data) if return_meta_data else (data, y)
+
+
+def _load_header_info(file):
+    """Load the meta data from a .ts file and advance file to the data.
+
+    Parameters
+    ----------
+    file : stream.
+        input file to read header from, assumed to be just opened
+
+    Returns
+    -------
+    meta_data : dict.
+        dictionary with the data characteristics stored in the header.
+    """
+    meta_data = {
+        "problemname": "none",
+        "timestamps": False,
+        "missing": False,
+        "univariate": True,
+        "equallength": True,
+        "classlabel": True,
+        "targetlabel": False,
+        "class_values": [],
+    }
+    boolean_keys = ["timestamps", "missing", "univariate", "equallength", "targetlabel"]
+    for line in file:
+        line = line.strip().lower()
+        if line and not line.startswith("#"):
+            tokens = line.split(" ")
+            token_len = len(tokens)
+            key = tokens[0][1:]
+            if key == "data":
+                if line != "@data":
+                    raise OSError("data tag should not have an associated value")
+                return meta_data
+            if key in meta_data.keys():
+                if key in boolean_keys:
+                    if token_len != 2:
+                        raise OSError(f"{tokens[0]} tag requires a boolean value")
+                    if tokens[1] == "true":
+                        meta_data[key] = True
+                    elif tokens[1] == "false":
+                        meta_data[key] = False
+                elif key == "problemname":
+                    meta_data[key] = tokens[1]
+                elif key == "classlabel":
+                    if tokens[1] == "true":
+                        meta_data["classlabel"] = True
+                        if token_len == 2:
+                            raise OSError(
+                                "if the classlabel tag is true then class values "
+                                "must be supplied"
+                            )
+                    elif tokens[1] == "false":
+                        meta_data["classlabel"] = False
+                    else:
+                        raise OSError("invalid class label value")
+                    meta_data["class_values"] = [token.strip() for token in tokens[2:]]
+        if meta_data["targetlabel"]:
+            meta_data["classlabel"] = False
+    return meta_data
+
+
+def _load_data(file, meta_data, replace_missing_vals_with="NaN"):
+    """Load data from a file with no header.
+
+    this assumes each time series has the same number of channels, but allows unequal
+    length series between cases.
+
+    Parameters
+    ----------
+    file : stream, input file to read data from, assume no comments or header info
+    meta_data : dict.
+        with meta data in the file header loaded with _load_header_info
+
+    Returns
+    -------
+    data: list[np.ndarray].
+        list of numpy arrays of floats: the time series
+    y_values : np.ndarray.
+        numpy array of strings: the class/target variable values
+    meta_data :  dict.
+        dictionary of characteristics enhanced with number of channels and series length
+        "problemname" (string), booleans: "timestamps", "missing", "univariate",
+        "equallength", "classlabel", "targetlabel" and "class_values": [],
+
+    """
+    data = []
+    n_cases = 0
+    n_channels = 0  # Assumed the same for all
+    current_channels = 0
+    series_length = 0
+    y_values = []
+    for line in file:
+        line = line.strip().lower()
+        line = line.replace("?", replace_missing_vals_with)
+        channels = line.split(":")
+        n_cases += 1
+        current_channels = len(channels)
+        if meta_data["classlabel"] or meta_data["targetlabel"]:
+            current_channels -= 1
+        if n_cases == 1:  # Find n_channels and length  from first if not unequal
+            n_channels = current_channels
+            if meta_data["equallength"]:
+                series_length = len(channels[0].split(","))
+        else:
+            if current_channels != n_channels:
+                raise OSError(
+                    f"Inconsistent number of dimensions in case {n_cases}. "
+                    f"Expecting {n_channels} but have read {current_channels}"
+                )
+            if meta_data["univariate"]:
+                if current_channels > 1:
+                    raise OSError(
+                        f"Seen {current_channels} in case {n_cases}."
+                        f"Expecting univariate from meta data"
+                    )
+        if meta_data["equallength"]:
+            current_length = series_length
+        else:
+            current_length = len(channels[0].split(","))
+        np_case = np.zeros(shape=(n_channels, current_length))
+        for i in range(0, n_channels):
+            single_channel = channels[i].strip()
+            data_series = single_channel.split(",")
+            data_series = [float(x) for x in data_series]
+            if len(data_series) != current_length:
+                raise OSError(
+                    f"Unequal length series, in case {n_cases} meta "
+                    f"data specifies all equal {series_length} but saw "
+                    f"{len(single_channel)}"
+                )
+            np_case[i] = np.array(data_series)
+        data.append(np_case)
+        if meta_data["classlabel"] or meta_data["targetlabel"]:
+            y_values.append(channels[n_channels])
+    if meta_data["equallength"]:
+        data = np.array(data)
+    return data, np.asarray(y_values), meta_data
diff --git a/sktime/libs/momentfm/utils/forecasting_metrics.py b/sktime/libs/momentfm/utils/forecasting_metrics.py
new file mode 100644
index 000000000..56fa07b79
--- /dev/null
+++ b/sktime/libs/momentfm/utils/forecasting_metrics.py
@@ -0,0 +1,364 @@
+"""Adapted from https://github.com/Nixtla/datasetsforecast/blob/main/datasetsforecast/losses.py."""
+
+import warnings
+from dataclasses import dataclass
+from typing import Optional, Union
+
+import numpy as np
+import numpy.typing as npt
+from skbase.utils.dependencies import _check_soft_dependencies
+
+from .utils import _reduce
+
+if _check_soft_dependencies(["torch"], severity="none"):
+    import torch
+    import torch.nn.functional as F
+    from torch import Tensor
+    from torch.nn.modules.loss import _Loss
+
+    class sMAPELoss(_Loss):
+        """sMAPELoss class."""
+
+        __constants__ = ["reduction"]
+
+        def __init__(
+            self, size_average=None, reduce=None, reduction: str = "mean"
+        ) -> None:
+            super().__init__(size_average, reduce, reduction)
+
+        def _abs(self, input):
+            return F.l1_loss(input, torch.zeros_like(input), reduction="none")
+
+        def _divide_no_nan(self, a: float, b: float) -> float:
+            """Auxiliary funtion to handle divide by 0."""
+            div = a / b
+            div[div != div] = 0.0
+            div[div == float("inf")] = 0.0
+            return div
+
+        def forward(self, input: Tensor, target: Tensor) -> Tensor:
+            """Forward function."""
+            delta_y = self._abs(input - target)
+            scale = self._abs(target) + self._abs(input)
+            error = self._divide_no_nan(delta_y, scale)
+            error = 200 * torch.nanmean(error)
+
+            return error
+
+    def _divide_no_nan(a: float, b: float) -> float:
+        """Auxiliary funtion to handle divide by 0."""
+        div = a / b
+        div[div != div] = 0.0
+        div[div == float("inf")] = 0.0
+        return div
+
+    @dataclass
+    class ForecastingMetrics:
+        """Forecasting Metrics."""
+
+        mae: Union[float, np.ndarray] = None
+        mse: Union[float, np.ndarray] = None
+        mape: Union[float, np.ndarray] = None
+        smape: Union[float, np.ndarray] = None
+        rmse: Union[float, np.ndarray] = None
+
+    def mae(
+        y: np.ndarray,
+        y_hat: np.ndarray,
+        reduction: str = "mean",
+        axis: Optional[int] = None,
+    ) -> Union[float, np.ndarray]:
+        r"""Calculate MAE.
+
+        Calculates Mean Absolute Error (MAE) between
+        y and y_hat. MAE measures the relative prediction
+        accuracy of a forecasting method by calculating the
+        deviation of the prediction and the true
+        value at a given time and averages these devations
+        over the length of the series.
+
+        $$ \\mathrm{MAE}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}_{\\tau}) =
+            \\frac{1}{H} \\sum^{t+H}_{\\tau=t+1}
+            |y_{\\tau} - \\hat{y}_{\\tau}| $$
+
+        Parameters
+        ----------
+            y: numpy array.
+                Observed values.
+            y_hat: numpy array
+                Predicted values.
+            reduction: str, optional.
+                Type of reduction to apply to the output: 'none' | 'mean' | 'sum'.
+                'none': no reduction will be applied,
+                'mean': the sum of the output will be divided by the number of
+                elements in the output,
+                'sum': the output will be summed.
+            axis: None or int, optional.
+                Axis or axes along which to average a.
+                The default, axis=None, will average over all of the elements of
+                the input array. If axis is negative it counts from last to first.
+
+        Returns
+        -------
+            mae: numpy array or double.
+                Return the MAE along the specified axis.
+        """
+        delta_y = np.abs(y - y_hat)
+        return _reduce(delta_y, reduction=reduction, axis=axis)
+
+    def mse(
+        y: np.ndarray,
+        y_hat: np.ndarray,
+        reduction: str = "mean",
+        axis: Optional[int] = None,
+    ) -> Union[float, np.ndarray]:
+        r"""Calculate MSE.
+
+        Calculates Mean Squared Error (MSE) between
+        y and y_hat. MSE measures the relative prediction
+        accuracy of a forecasting method by calculating the
+        squared deviation of the prediction and the true
+        value at a given time, and averages these devations
+        over the length of the series.
+
+        $$ \\mathrm{MSE}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}_{\\tau}) =
+            \\frac{1}{H} \\sum^{t+H}_{\\tau=t+1} (y_{\\tau} -
+            \\hat{y}_{\\tau})^{2} $$
+
+        Parameters
+        ----------
+            y: numpy array.
+                Actual test values.
+            y_hat: numpy array.
+                Predicted values.
+            reduction: str, optional.
+                Type of reduction to apply to the output: 'none' | 'mean' | 'sum'.
+                'none': no reduction will be applied,
+                'mean': the sum of the output will be divided by the number of
+                elements in the output,
+                'sum': the output will be summed.
+            axis: None or int, optional.
+                Axis or axes along which to average a.
+                The default, axis=None, will average over all of the
+                elements of the input array. If axis is negative it counts
+                from the last to the first axis.
+
+        Returns
+        -------
+            mse: numpy array or double.
+                Return the MSE along the specified axis.
+        """
+        delta_y = np.square(y - y_hat)
+        return _reduce(delta_y, reduction=reduction, axis=axis)
+
+    def rmse(
+        y: np.ndarray,
+        y_hat: np.ndarray,
+        reduction: str = "mean",
+        axis: Optional[int] = None,
+    ) -> Union[float, np.ndarray]:
+        r"""Calculate Rmse.
+
+        Calculates Root Mean Squared Error (RMSE) between
+        y and y_hat. RMSE measures the relative prediction
+        accuracy of a forecasting method by calculating the squared deviation
+        of the prediction and the observed value at a given time and
+        averages these devations over the length of the series.
+        Finally the RMSE will be in the same scale
+        as the original time series so its comparison with other
+        series is possible only if they share a common scale.
+        RMSE has a direct connection to the L2 norm.
+
+        $$ \\mathrm{RMSE}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}_{\\tau}) =
+            \\sqrt{\\frac{1}{H} \\sum^{t+H}_{\\tau=t+1} (y_{\\tau} -
+            \\hat{y}_{\\tau})^{2}} $$
+
+        Parameters
+        ----------
+            y: numpy array.
+                Observed values.
+            y_hat: numpy array.
+                Predicted values.
+            reduction: str, optional.
+                Type of reduction to apply to the output: 'none' | 'mean' | 'sum'.
+                'none': no reduction will be applied,
+                'mean': the sum of the output will be divided by the number of
+                elements in the output,
+                'sum': the output will be summed.
+            axis: None or int, optional.
+                Axis or axes along which to average a.
+                The default, axis=None, will average over all of the elements of
+                the input array. If axis is negative it counts from the last to
+                first.
+
+        Returns
+        -------
+            rmse: numpy array or double.
+                Return the RMSE along the specified axis.
+        """
+        return np.sqrt(mse(y, y_hat, reduction, axis))
+
+    def mape(
+        y: np.ndarray,
+        y_hat: np.ndarray,
+        reduction: str = "mean",
+        axis: Optional[int] = None,
+    ) -> Union[float, np.ndarray]:
+        r"""Calculate mape.
+
+        Calculates Mean Absolute Percentage Error (MAPE) between
+        y and y_hat. MAPE measures the relative prediction
+        accuracy of a forecasting method by calculating the percentual deviation
+        of the prediction and the observed value at a given time and
+        averages these devations over the length of the series.
+        The closer to zero an observed value is, the higher penalty MAPE loss
+        assigns to the corresponding error.
+
+        $$ \\mathrm{MAPE}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}_{\\tau}) =
+            \\frac{1}{H} \\sum^{t+H}_{\\tau=t+1}
+            \\frac{|y_{\\tau}-\\hat{y}_{\\tau}|}{|y_{\\tau}|} $$
+
+        Parameters
+        ----------
+            y: numpy array.
+                Observed values.
+            y_hat: numpy array.
+                Predicted values.
+            reduction: str, optional.
+                Type of reduction to apply to the output: 'none' | 'mean' | 'sum'.
+                'none': no reduction will be applied,
+                'mean': the sum of the output will be divided by the number of
+                elements
+                in the output,
+                'sum': the output will be summed.
+            axis: None or int, optional.
+                Axis or axes along which to average a.
+                The default, axis=None, will average over all of the elements of
+                the input array. If axis is negative it counts from the last to
+                first.
+
+        Returns
+        -------
+            mape: numpy array or double.
+                Return the MAPE along the specified axis.
+        """
+        delta_y = np.abs(y - y_hat)
+        scale = np.abs(y)
+        error = _divide_no_nan(delta_y, scale)
+        return 100 * _reduce(error, reduction=reduction, axis=axis)
+
+    def smape(
+        y: np.ndarray,
+        y_hat: np.ndarray,
+        reduction: str = "mean",
+        axis: Optional[int] = None,
+    ) -> Union[float, np.ndarray]:
+        r"""Calculate smape.
+
+        Calculates Symmetric Mean Absolute Percentage Error (SMAPE) between
+        y and y_hat. SMAPE measures the relative prediction
+        accuracy of a forecasting method by calculating the relative deviation
+        of the prediction and the observed value scaled by the sum of the
+        absolute values for the prediction and observed value at a
+        given time, then averages these devations over the length
+        of the series. This allows the SMAPE to have bounds between
+        0% and 200% which is desireble compared to normal MAPE that
+        may be undetermined when the target is zero.
+
+        $$ \\mathrm{SMAPE}_{2}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}_{\\tau}) =
+        \\frac{1}{H} \\sum^{t+H}_{\\tau=t+1}
+        \\frac{|y_{\\tau}-\\hat{y}_{\\tau}|}{|y_{\\tau}|+|\\hat{y}_{\\tau}|} $$
+
+        Parameters
+        ----------
+            y: numpy array.
+                Observed values.
+            y_hat: numpy array.
+                Predicted values.
+            reduction: str, optional.
+                Type of reduction to apply to the output: 'none' | 'mean' | 'sum'.
+                'none': no reduction will be applied,
+                'mean': the sum of the output will be divided by the number of
+                elements
+                in the output,
+                'sum': the output will be summed.
+            axis: None or int, optional.
+                Axis or axes along which to average a.
+                The default, axis=None, will average over all of the elements of
+                the input array. If axis is negative it counts from the last to
+                first.
+
+        Returns
+        -------
+            smape: numpy array or double.
+                Return the SMAPE along the specified axis.
+        """
+        delta_y = np.abs(y - y_hat)
+        scale = np.abs(y) + np.abs(y_hat)
+        error = _divide_no_nan(delta_y, scale)
+        error = 200 * _reduce(error, reduction=reduction, axis=axis)
+
+        if isinstance(error, float):
+            # assert error <= 200, "SMAPE should be lower than 200"
+            if error > 200:
+                warnings.warn(
+                    f"SMAPE should be lower than 200 but was found to be {error}"
+                )
+        else:
+            # assert all(error <= 200), "SMAPE should be lower than 200"
+            if all(error > 200):
+                warnings.warn(
+                    f"SMAPE should be lower than 200 but was found to be {error}"
+                )
+
+        return error
+
+    def get_forecasting_metrics(
+        y: npt.NDArray,
+        y_hat: npt.NDArray,
+        reduction: str = "mean",
+        axis: Optional[int] = None,
+    ) -> Union[float, np.ndarray]:
+        """Get Forecasting Metrics."""
+        return ForecastingMetrics(
+            mae=mae(y=y, y_hat=y_hat, axis=axis, reduction=reduction),
+            mse=mse(y=y, y_hat=y_hat, axis=axis, reduction=reduction),
+            mape=mape(y=y, y_hat=y_hat, axis=axis, reduction=reduction),
+            smape=smape(y=y, y_hat=y_hat, axis=axis, reduction=reduction),
+            rmse=rmse(y=y, y_hat=y_hat, axis=axis, reduction=reduction),
+        )
+else:
+
+    class sMAPELoss:
+        """Create dummy class if torch is unavailable."""
+
+        pass
+
+    class ForecastingMetrics:
+        """Create dummy class torch is unavailable."""
+
+        pass
+
+    def mae():
+        """Create dummy function if torch is not available."""
+        pass
+
+    def mse():
+        """Create dummy function if torch is not available."""
+        pass
+
+    def rmse():
+        """Create dummy function if torch is not available."""
+        pass
+
+    def mape():
+        """Create dummy function if torch is not available."""
+        pass
+
+    def smape():
+        """Create dummy function if torch is not available."""
+        pass
+
+    def get_forecasting_metrics():
+        """Create dummy function if torch is not available."""
+        pass
diff --git a/sktime/libs/momentfm/utils/masking.py b/sktime/libs/momentfm/utils/masking.py
new file mode 100644
index 000000000..8e840a534
--- /dev/null
+++ b/sktime/libs/momentfm/utils/masking.py
@@ -0,0 +1,139 @@
+"""Masking module for momentfm."""
+
+from typing import Optional
+
+from skbase.utils.dependencies import _check_soft_dependencies
+
+if _check_soft_dependencies(["torch"], severity="none"):
+    import torch
+
+    class Masking:
+        """Masking Module."""
+
+        def __init__(
+            self,
+            mask_ratio: float = 0.3,
+            patch_len: int = 8,
+            stride: Optional[int] = None,
+        ):
+            """Masking class for momentfm.
+
+            Indices with 0 mask are hidden, and with 1 are observed.
+            """
+            self.mask_ratio = mask_ratio
+            self.patch_len = patch_len
+            self.stride = patch_len if stride is None else stride
+
+        @staticmethod
+        def convert_seq_to_patch_view(
+            mask: torch.Tensor, patch_len: int = 8, stride: Optional[int] = None
+        ):
+            """Convert sequence to patch function.
+
+            Input:
+                mask : torch.Tensor of shape [batch_size x seq_len]
+            Output
+                mask : torch.Tensor of shape [batch_size x n_patches]
+            """
+            stride = patch_len if stride is None else stride
+            mask = mask.unfold(dimension=-1, size=patch_len, step=stride)
+            # mask : [batch_size x n_patches x patch_len]
+            return (mask.sum(dim=-1) == patch_len).long()
+
+        @staticmethod
+        def convert_patch_to_seq_view(
+            mask: torch.Tensor,
+            patch_len: int = 8,
+        ):
+            """Convert patch to sequence function.
+
+            Input:
+                mask : torch.Tensor of shape [batch_size x n_patches]
+            Output:
+                mask : torch.Tensor of shape [batch_size x seq_len]
+            """
+            return mask.repeat_interleave(patch_len, dim=-1)
+
+        def generate_mask(
+            self, x: torch.Tensor, input_mask: Optional[torch.Tensor] = None
+        ):
+            """Generate Mask Function.
+
+            Input:
+                x : torch.Tensor of shape
+                [batch_size x n_channels x n_patches x patch_len] or
+                [batch_size x n_channels x seq_len]
+                input_mask: torch.Tensor of shape [batch_size x seq_len] or
+                [batch_size x n_patches]
+            Output:
+                mask : torch.Tensor of shape [batch_size x seq_len]
+            """
+            if x.ndim == 4:
+                return self._mask_patch_view(x, input_mask=input_mask)
+            elif x.ndim == 3:
+                return self._mask_seq_view(x, input_mask=input_mask)
+
+        def _mask_patch_view(self, x, input_mask=None):
+            """Return mask for inputted patch view.
+
+            Input:
+                x : torch.Tensor of shape
+                [batch_size x n_channels x n_patches x patch_len]
+                input_mask: torch.Tensor of shape [batch_size x seq_len]
+            Output:
+                mask : torch.Tensor of shape [batch_size x n_patches]
+            """
+            input_mask = self.convert_seq_to_patch_view(
+                input_mask, self.patch_len, self.stride
+            )
+            n_observed_patches = input_mask.sum(dim=-1, keepdim=True)  # batch_size x 1
+
+            batch_size, _, n_patches, _ = x.shape
+            len_keep = torch.ceil(n_observed_patches * (1 - self.mask_ratio)).long()
+            noise = torch.rand(
+                batch_size, n_patches, device=x.device
+            )  # noise in [0, 1], batch_size x n_channels x n_patches
+            noise = torch.where(
+                input_mask == 1, noise, torch.ones_like(noise)
+            )  # only keep the noise of observed patches
+
+            # Sort noise for each sample
+            ids_shuffle = torch.argsort(
+                noise, dim=1
+            )  # Ascend: small is keep, large is remove
+            ids_restore = torch.argsort(
+                ids_shuffle, dim=1
+            )  # ids_restore: [batch_size x n_patches]
+
+            # Generate the binary mask: 0 is keep, 1 is remove
+            mask = torch.zeros(
+                [batch_size, n_patches], device=x.device
+            )  # mask: [batch_size x n_patches]
+            for i in range(batch_size):
+                mask[i, : len_keep[i]] = 1
+
+            # Unshuffle to get the binary mask
+            mask = torch.gather(mask, dim=1, index=ids_restore)
+
+            return mask.long()
+
+        def _mask_seq_view(self, x, input_mask=None):
+            """Return mask for inputted tensor.
+
+            Input:
+                x : torch.Tensor of shape
+                [batch_size x n_channels x seq_len]
+                input_mask: torch.Tensor of shape [batch_size x seq_len]
+            Output:
+                mask : torch.Tensor of shape [batch_size x seq_len]
+            """
+            x = x.unfold(dimension=-1, size=self.patch_len, step=self.stride)
+            mask = self._mask_patch_view(x, input_mask=input_mask)
+            return self.convert_patch_to_seq_view(mask, self.patch_len).long()
+
+else:
+
+    class Masking:
+        """Dummy class if torch is unavailable."""
+
+        pass
diff --git a/sktime/libs/momentfm/utils/utils.py b/sktime/libs/momentfm/utils/utils.py
new file mode 100644
index 000000000..2ca94ba17
--- /dev/null
+++ b/sktime/libs/momentfm/utils/utils.py
@@ -0,0 +1,162 @@
+"""Utility Function for momentfm."""
+
+import os
+import random
+from argparse import Namespace
+
+import numpy as np
+from skbase.utils.dependencies import _check_soft_dependencies
+
+if _check_soft_dependencies(["torch", "transformers"], severity="none"):
+    import torch
+
+    def control_randomness(seed: int = 13):
+        """Control randomness seed function."""
+        random.seed(seed)
+        os.environ["PYTHONHASHSEED"] = str(seed)
+        np.random.seed(seed)
+        torch.manual_seed(seed)
+        torch.cuda.manual_seed(seed)
+        torch.backends.cudnn.deterministic = True
+        torch.backends.cudnn.benchmark = False
+
+    def dtype_map(dtype: str):
+        """Check dtype map."""
+        map = {
+            "float16": torch.float16,
+            "float32": torch.float32,
+            "float64": torch.float64,
+            "bfloat16": torch.bfloat16,
+            "uint8": torch.uint8,
+            "int8": torch.int8,
+            "int16": torch.int16,
+            "int32": torch.int32,
+            "int64": torch.int64,
+            "bool": torch.bool,
+        }
+        return map[dtype]
+
+    def get_huggingface_model_dimensions(model_name: str = "flan-t5-base"):
+        """Get HF Model dimensions."""
+        from transformers import T5Config
+
+        config = T5Config.from_pretrained(model_name)
+
+        return config.d_model
+
+    def get_anomaly_criterion(anomaly_criterion: str = "mse"):
+        """Get anomaly criterion."""
+        if anomaly_criterion == "mse":
+            return torch.nn.MSELoss(reduction="none")
+        elif anomaly_criterion == "mae":
+            return torch.nn.L1Loss(reduction="none")
+        else:
+            raise ValueError(f"Anomaly criterion {anomaly_criterion} not supported.")
+
+    def _reduce(metric, reduction="mean", axis=None):
+        if reduction == "mean":
+            return np.nanmean(metric, axis=axis)
+        elif reduction == "sum":
+            return np.nansum(metric, axis=axis)
+        elif reduction == "none":
+            return metric
+
+    class EarlyStopping:
+        """Class to implement early stopping."""
+
+        def __init__(self, patience: int = 3, verbose: bool = False, delta: float = 0):
+            self.patience = patience
+            self.verbose = verbose
+            self.counter = 0
+            self.best_score = None
+            self.early_stop = False
+            self.val_loss_min = np.Inf
+            self.delta = delta
+
+        def __call__(self, validation_loss):
+            """Call method for early stopping."""
+            score = -validation_loss
+            if self.best_score is None:
+                self.best_score = score
+
+            elif score < self.best_score + self.delta:
+                self.counter += 1
+                if self.verbose:
+                    print(
+                        f"EarlyStopping counter: {self.counter} out of {self.patience}"
+                    )
+                if self.counter >= self.patience:
+                    self.early_stop = True
+            else:
+                self.best_score = score
+                self.counter = 0
+
+else:
+
+    class EarlyStopping:
+        """Dummy class if torch is unavailable."""
+
+        pass
+
+    def get_anomaly_criterion(anomaly_criterion: str = "mse"):
+        """Get anomaly criterion."""
+        pass
+
+    def control_randomness(seed: int = 13):
+        """Control randomness seed function."""
+        pass
+
+    def dtype_map(dtype: str):
+        """Check dtype map."""
+        pass
+
+    def get_huggingface_model_dimensions(model_name: str = "flan-t5-base"):
+        """Get HF Model dimensions."""
+        pass
+
+    def _reduce(metric, reduction="mean", axis=None):
+        pass
+
+
+class NamespaceWithDefaults(Namespace):
+    """Namespace class."""
+
+    @classmethod
+    def from_namespace(cls, namespace):
+        """Set attribute function."""
+        new_instance = cls()
+        for attr in dir(namespace):
+            if not attr.startswith("__"):
+                setattr(new_instance, attr, getattr(namespace, attr))
+        return new_instance
+
+    def getattr(self, key, default=None):
+        """Get attr function."""
+        return getattr(self, key, default)
+
+
+def parse_config(config: dict) -> NamespaceWithDefaults:
+    """Parse Config function."""
+    args = NamespaceWithDefaults(**config)
+    return args
+
+
+def make_dir_if_not_exists(path, verbose=True):
+    """Make Directory if does not exist."""
+    if not is_directory(path):
+        path = path.split(".")[0]
+    if not os.path.exists(path=path):
+        os.makedirs(path)
+        if verbose:
+            print(f"Making directory: {path}...")
+    return True
+
+
+def is_directory(path):
+    """Check is_directory."""
+    extensions = [".pth", ".txt", ".json", ".yaml"]
+
+    for ext in extensions:
+        if ext in path:
+            return False
+    return True


Test patch:
diff --git a/sktime/forecasting/tests/test_known_values.py b/sktime/forecasting/tests/test_known_values.py
new file mode 100644
index 000000000..b8c10d7d6
--- /dev/null
+++ b/sktime/forecasting/tests/test_known_values.py
@@ -0,0 +1,78 @@
+#!/usr/bin/env python3 -u
+"""Tests for ForecastKnownValues."""
+# copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
+
+__author__ = ["RobKuebler"]
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from sktime.forecasting.dummy import ForecastKnownValues
+from sktime.forecasting.tests._config import TEST_OOS_FHS
+from sktime.tests.test_switch import run_test_for_class
+
+
+@pytest.fixture
+def y_known():
+    index = pd.MultiIndex.from_product(
+        [["A", "B"], ["X"], [0, 1, 2]],
+        names=["Level1", "Level2", "Date"],
+    )
+    data = range(len(index))
+
+    return pd.DataFrame(data, index=index, columns=["Value"])
+
+
+@pytest.mark.skipif(
+    not run_test_for_class(ForecastKnownValues),
+    reason="run test only if softdeps are present and incrementally (if requested)",
+)
+@pytest.mark.parametrize("fh", TEST_OOS_FHS)
+def test_multiindex(fh, y_known) -> None:
+    """Test multiindex y_known."""
+    f = ForecastKnownValues(y_known=y_known)
+    f.fit(y_known)
+    y_pred = f.predict(fh)
+
+    # Create expected data
+    if not isinstance(fh, np.ndarray):
+        fh = np.array([fh])
+
+    index = pd.MultiIndex.from_product(
+        [["A", "B"], ["X"], fh + 2],
+        names=["Level1", "Level2", "Date"],
+    )
+    expected = pd.DataFrame(None, index=index, columns=["Value"])
+
+    pd.testing.assert_frame_equal(
+        y_pred,
+        expected,
+        check_dtype=False,
+        check_index_type=False,
+    )
+
+
+@pytest.mark.skipif(
+    not run_test_for_class(ForecastKnownValues),
+    reason="run test only if softdeps are present and incrementally (if requested)",
+)
+@pytest.mark.parametrize("fh", TEST_OOS_FHS)
+def test_singleindex(fh, y_known) -> None:
+    """Test singleindex y_known."""
+    f = ForecastKnownValues(y_known=y_known.loc["A", "X"])
+    f.fit(y_known.loc["A", "X"])
+    y_pred = f.predict(fh)
+
+    # Create expected data
+    if not isinstance(fh, np.ndarray):
+        fh = np.array([fh])
+
+    expected = pd.DataFrame(None, index=fh + 2, columns=["Value"])
+
+    pd.testing.assert_frame_equal(
+        y_pred,
+        expected,
+        check_dtype=False,
+        check_index_type=False,
+    )
diff --git a/sktime/tests/_config.py b/sktime/tests/_config.py
index df62d7bb6..5aecbe942 100644
--- a/sktime/tests/_config.py
+++ b/sktime/tests/_config.py
@@ -49,6 +49,8 @@ EXCLUDE_ESTIMATORS = [
     # sporadic timeouts, see #6344
     "ShapeletLearningClassifierTslearn",
     "DartsXGBModel",
+    # models with large weights
+    "MomentFMForecaster",
     # Large datasets
     "M5Dataset",
     # Test estimators