From 2be83da27d20691cfe0e6dd5157f5411504ef864 Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Tue, 9 Dec 2025 09:24:04 +0100
Subject: [PATCH 01/13] wip

---
 .../models/forecasting/median_forecaster.py   | 431 ++++++++++++++++++
 1 file changed, 431 insertions(+)
 create mode 100644 packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py

diff --git a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
new file mode 100644
index 000000000..5ba6554c5
--- /dev/null
+++ b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
@@ -0,0 +1,431 @@
+# SPDX-FileCopyrightText: 2025 Contributors to the OpenSTEF project <short.term.energy.forecasts@alliander.com>
+#
+# SPDX-License-Identifier: MPL-2.0
+
+"""Median regressor based forecasting models for energy forecasting.
+
+Provides median regression models for multi-quantile energy forecasting.
+Optimized for time series data with specialized loss functions and
+comprehensive hyperparameter control for production forecasting workflows.
+
+Note that this is a autoregressive model, meaning that it uses the previous
+    predictions to predict the next value.
+
+    This regressor is good for predicting two types of signals:
+    - Signals with very slow dynamics compared to the sampling rate, possibly
+      with a lot of noise.
+    - Signals that switch between two or more states, which random in nature or
+    depend on unknown features, but tend to be stable in each state. An example of
+    this may be waste heat delivered from an industrial process. Using a median
+    over the last few timesteps adds some hysterisis to avoid triggering on noise.
+
+    Tips for using this regressor:
+    - Set the lags to be evenly spaced and at a frequency mathching the
+      frequency of the input data. For example, if the input data is at 15
+      minute intervals, set the lags to be at 15 minute intervals as well.
+    - Use a small training dataset, since there are no actual parameters to train.
+    - Set the frequency of the input data index to avoid inferring it. Inference might be
+    a problem if we get very small chunks of data in training or validation sets.
+    - Use only one training horizon, since the regressor will use the same lags for all
+      training horizons.
+    - Allow for missing data by setting completeness_threshold to 0. If the prediction horizon
+      is larger than the context window there will be a lot of nans in the input data, but
+      the autoregression solves that.
+"""
+
+from datetime import timedelta
+from typing import Self, override
+
+import numpy as np
+import pandas as pd
+from pydantic import Field
+
+from openstef_core.datasets.validated_datasets import ForecastDataset, ForecastInputDataset
+from openstef_core.exceptions import ModelLoadingError
+from openstef_core.mixins import State
+from openstef_core.mixins.predictor import HyperParams
+from openstef_core.types import LeadTime, Quantile
+from openstef_models.explainability.mixins import ExplainableForecaster
+from openstef_models.models.forecasting.forecaster import Forecaster, ForecasterConfig
+
+
+class MedianForecasterHyperParams(HyperParams):
+    """Hyperparameter configuration for base case forecaster."""
+
+    primary_lag: timedelta = Field(
+        default=timedelta(days=7),
+        description="Primary lag to use for predictions (default: 7 days for weekly patterns)",
+    )
+    fallback_lag: timedelta = Field(
+        default=timedelta(days=14),
+        description="Fallback lag to use when primary lag data is unavailable (default: 14 days)",
+    )
+
+
+class MedianForecasterConfig(ForecasterConfig):
+    """Configuration for base case forecaster."""
+
+    quantiles: list[Quantile] = Field(
+        default=[Quantile(0.5)],
+        description=(
+            "Probability levels for uncertainty estimation. Each quantile represents a confidence level "
+            "(e.g., 0.1 = 10th percentile, 0.5 = median, 0.9 = 90th percentile). "
+            "Models must generate predictions for all specified quantiles."
+        ),
+        min_length=1,
+        max_length=1,
+    )
+    horizons: list[LeadTime] = Field(
+        default=...,
+        description=(
+            "Lead times for predictions, accounting for data availability and versioning cutoffs. "
+            "Each horizon defines how far ahead the model should predict."
+        ),
+        min_length=1,
+        max_length=1,
+    )
+
+    hyperparams: MedianForecasterHyperParams = Field(
+        default_factory=MedianForecasterHyperParams,
+    )
+
+class MedianForecaster(Forecaster, ExplainableForecaster):
+    """Median forecaster using lag features for predictions.
+
+    This forecaster predicts the median value based on specified lag features.
+    It is particularly useful for signals with slow dynamics or state-switching behavior.
+
+    Hyperparameters:
+        lags: List of time deltas representing the lag features to use for prediction.
+              These should be aligned with the data sampling frequency.
+        context_window: Time delta representing the context window size for input data.
+                        This defines how much historical data is considered for making predictions.
+
+                _config: BaseCaseForecasterConfig
+    """
+
+    def __init__(
+        self,
+        config: MedianForecasterConfig,
+    ) -> None:
+        """Initialize the base case forecaster.
+
+        Args:
+            config: Configuration specifying quantiles, horizons, and lag hyperparameters.
+                   If None, uses default configuration with 7-day primary and 14-day fallback lags.
+        """
+        self._config = config
+
+    @property
+    @override
+    def config(self) -> MedianForecasterConfig:
+        return self._config
+
+    @property
+    @override
+    def hyperparams(self) -> MedianForecasterHyperParams:
+        return self._config.hyperparams
+
+    @override
+    def to_state(self) -> State:
+        return {
+            "version": MODEL_CODE_VERSION,
+            "config": self.config.model_dump(mode="json"),
+        }
+
+    @override
+    def from_state(self, state: State) -> Self:
+        if not isinstance(state, dict) or "version" not in state or state["version"] > MODEL_CODE_VERSION:
+            raise ModelLoadingError("Invalid state for BaseCaseForecaster")
+
+        return self.__class__(config=MedianForecasterConfig.model_validate(state["config"]))
+
+    @property
+    @override
+    def is_fitted(self) -> bool:
+        return self._xgboost_model.__sklearn_is_fitted__() 
+
+
+
+    @property
+    @override
+    def feature_importances(self) -> pd.DataFrame:
+        booster = self._xgboost_model.get_booster()
+        weights_df = pd.DataFrame(
+            data=booster.get_score(importance_type="gain"),
+            index=[quantile.format() for quantile in self.config.quantiles],
+        ).transpose()
+        weights_df.index.name = "feature_name"
+        weights_df.columns.name = "quantiles"
+
+        weights_abs = weights_df.abs()
+        total = weights_abs.sum(axis=0).replace(to_replace=0, value=1.0)  # pyright: ignore[reportUnknownMemberType]
+
+        return weights_abs / total
+
+
+
+    @staticmethod
+    def _infer_frequency(index: pd.DatetimeIndex) -> pd.Timedelta:
+        """
+        Infer the frequency of a pandas DatetimeIndex if the freq attribute is not set.
+        This method calculates the most common time difference between consecutive timestamps,
+        which is more permissive of missing chunks of data than the pandas infer_freq method.
+
+        Args:
+            index (pd.DatetimeIndex): The datetime index to infer the frequency from.
+
+        Returns:
+            pd.Timedelta: The inferred frequency as a pandas Timedelta.
+        """
+        if len(index) < 2:
+            raise ValueError(
+                "Cannot infer frequency from an index with fewer than 2 timestamps."
+            )
+
+        # Calculate the differences between consecutive timestamps
+        deltas = index.to_series().diff().dropna()
+
+        # Find the most common difference
+        inferred_freq = deltas.mode().iloc[0]
+        return inferred_freq
+
+    def _frequency_matches(self, index: pd.DatetimeIndex) -> bool:
+        """
+        Check if the frequency of the input data matches the model frequency.
+
+        Args:
+            x (pd.DataFrame): The input data to check.
+
+        Returns:
+            bool: True if the frequencies match, False otherwise.
+        """
+        if not isinstance(index, pd.DatetimeIndex):
+            raise ValueError(
+                "The index of the input data must be a pandas DatetimeIndex."
+            )
+
+        if index.freq is None:
+            input_frequency = self._infer_frequency(index)
+        else:
+            input_frequency = index.freq
+
+        return input_frequency == pd.Timedelta(minutes=self.frequency)
+
+    @staticmethod
+    def _extract_and_validate_lags(
+        x: pd.DataFrame,
+    ) -> tuple[tuple[str], int, list[tuple[str, int]]]:
+        """Extract and validate the lag features from the input data.
+
+        This method checks that the lag features are evenly spaced and match the frequency of the input data.
+        It also extracts the lag features and their corresponding time deltas.
+        Args:
+            x (pd.DataFrame): The input data containing lag features.
+        Returns:
+            tuple: A tuple containing:
+                - A list of feature names, sorted by their lag in minutes.
+                - The frequency of the lag features in minutes.
+                - A list of tuples containing the lag feature names and their corresponding time deltas in minutes.
+        """
+        # Check that the input data contains the required lag features
+        feature_names = list(x.columns[x.columns.str.startswith("T-")])
+        if len(feature_names) == 0:
+            raise ValueError("No lag features found in the input data.")
+
+        # Convert all lags to minutes to make comparable
+        feature_to_lags_in_min = []
+        for feature in feature_names:
+            if feature.endswith("min"):
+                lag_in_min = int(feature.split("-")[1].split("min")[0])
+            elif feature.endswith("d"):
+                lag_in_min = int(feature.split("-")[1].split("d")[0]) * 60 * 24
+            else:
+                raise ValueError(
+                    f"Feature name '{feature}' does not follow the expected format."
+                    " Expected format is 'T-<lag_in_minutes>' or 'T-<lag_in_days>d'."
+                )
+            feature_to_lags_in_min.append((feature, lag_in_min))
+
+        # Sort the features by lag in minutes
+        feature_to_lags_in_min.sort(key=lambda x: x[1])
+        sorted_features, sorted_lags_in_min = zip(*feature_to_lags_in_min)
+
+        # Check that the lags are evenly spaced
+        diffs = np.diff(sorted_lags_in_min)
+        unique_diffs = np.unique(diffs)
+        if len(unique_diffs) > 1:
+            raise ValueError(
+                "Lag features are not evenly spaced. "
+                f"Got lags with differences: {unique_diffs} min. "
+                "Please ensure that the lag features are generated correctly."
+            )
+        frequency = unique_diffs[0]
+
+        return sorted_features, frequency, feature_to_lags_in_min
+
+    @staticmethod
+    def _fill_diagonal_with_median(
+        lag_array: np.ndarray, start: int, end: int, median: float
+    ) -> np.ndarray | None:
+        # Use the calculated median to fill in future lag values where this prediction would be used as input.
+
+        # If the start index is beyond the array bounds, no future updates are needed from this step.
+        if start >= lag_array.shape[0]:
+            return lag_array
+
+        # Ensure the end index does not exceed the array bounds.
+        end = min(end, lag_array.shape[0])
+
+        # Get a view of the sub-array where the diagonal needs to be filled.
+        # The slice represents future time steps (rows) and corresponding lag features (columns).
+        # Rows: from 'start' up to (but not including) 'end'
+        # Columns: from 0 up to (but not including) 'end - start'
+        # This selects the part of the array where lag_array[start + k, k] resides for k in range(end - start).
+        view = lag_array[start:end, 0 : (end - start)]
+
+        # Create a mask for NaNs on the diagonal
+        diagonal_nan_mask = np.isnan(np.diag(view))
+
+        # Only update if there are NaNs on the diagonal
+        if np.any(diagonal_nan_mask):
+            # Create a temporary array to hold the new diagonal
+            updated_diagonal = np.diag(view).copy()
+            updated_diagonal[diagonal_nan_mask] = median
+            np.fill_diagonal(view, updated_diagonal)
+
+    @override
+    def predict(self, data: ForecastInputDataset) -> ForecastDataset:
+        """
+        Predict the median of the lag features for each time step in the context window.
+
+        Args:
+            x (pd.DataFrame): The input data for prediction. This should be a pandas dataframe with lag features.
+
+        Returns:
+            np.array: The predicted median for each time step in the context window.
+            If any lag feature is NaN, this will be ignored.
+            If all lag features are NaN, the regressor will return NaN.
+        """
+
+        input_data: pd.DataFrame = data.input_data(start=data.forecast_start)
+
+        if not self._frequency_matches(input_data.index):
+            raise ValueError(
+                f"The input data frequency ({input_data.index.freq}) does not match the model frequency ({self.frequency})."
+            )
+
+        # Check that the input data contains the required lag features
+        missing_features = set(self.feature_names) - set(input_data.columns)
+        if missing_features:
+            raise ValueError(
+                f"The input data is missing the following lag features: {missing_features}"
+            )
+
+        # Reindex the input data to ensure there are no gaps in the time series.
+        # This is important for the autoregressive logic that follows.
+        # Store the original index to return predictions aligned with the input.
+        original_index = input_data.index.copy()
+        first_index = input_data.index[0]
+        last_index = input_data.index[-1]
+        freq = pd.Timedelta(minutes=self.frequency)
+        # Create a new date range with the expected frequency.
+        new_index = pd.date_range(first_index, last_index, freq=freq)
+        # Reindex the input DataFrame, filling any new timestamps with NaN.
+        input_data = input_data.reindex(new_index, fill_value=np.nan)
+
+        # Select only the lag feature columns in the specified order.
+        lag_df = input_data[self.feature_names]
+
+        # Convert the lag DataFrame and its index to NumPy arrays for faster processing.
+        lag_array = lag_df.to_numpy()
+        time_index = lag_df.index.to_numpy()
+        # Initialize the prediction array with NaNs.
+        prediction = np.full(lag_array.shape[0], np.nan)
+
+        # Calculate the time step size based on the model frequency.
+        step_size = pd.Timedelta(minutes=self.frequency)
+        # Determine the number of steps corresponding to the smallest and largest lags.
+        smallest_lag_steps = int(
+            self.lags_to_time_deltas_[self.feature_names[0]] / step_size
+        )
+        largest_lag_steps = int(
+            self.lags_to_time_deltas_[self.feature_names[-1]] / step_size
+        )
+
+        # Iterate through each time step in the reindexed data.
+        for time_step in range(lag_array.shape[0]):
+            # Get the lag features for the current time step.
+            current_lags = lag_array[time_step]
+            # Calculate the median of the available lag features, ignoring NaNs.
+            median = np.nanmedian(current_lags)
+            # Store the calculated median in the prediction array.
+            prediction[time_step] = median
+
+            # If the median calculation resulted in NaN (e.g., all lags were NaN), skip the autoregression step.
+            if np.isnan(median):
+                continue
+
+            # Auto-regressive step: update the lag array for future time steps.
+            # Calculate the start and end indices in the future time steps that will be affected.
+            start, end = (
+                time_step + smallest_lag_steps,
+                time_step + largest_lag_steps + 1,
+            )
+            self._fill_diagonal_with_median(lag_array, start, end, median)
+
+        # Convert the prediction array back to a pandas DataFrame using the reindexed time index.
+        prediction_df = pd.DataFrame(prediction, index=time_index, columns=["median"])
+        # Select only the predictions corresponding to the original input index.
+        prediction = prediction_df.loc[original_index].to_numpy().flatten()
+
+        # Return the final predictions as a ForecastDataset.
+        predictions = pd.DataFrame(
+            data=prediction,
+            index=input_data.index,
+            columns=[quantile.format() for quantile in self.config.quantiles],
+        )
+
+        return ForecastDataset(
+            data=predictions,
+            sample_interval=data.sample_interval,
+        )
+
+    @override
+    def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None = None) -> None:
+                """This model does not have any hyperparameters to fit,
+        but it does need to know the feature names of the lag features and the order of these.
+
+        Lag features are expected to be evently spaced and match the frequency of the input data.
+        The lag features are expected to be named in the format T-<lag_in_minutes> or T-<lag_in_days>d.
+        For example, T-1min, T-2min, T-3min or T-1d, T-2d.
+
+        Which lag features are used is determined by the feature engineering step.
+        """                
+        input_data_predictors = data.input_data()
+    
+        target_series = data.target_series
+    
+        (
+            feature_names,
+            frequency,
+            feature_to_lags_in_min,
+        ) = self._extract_and_validate_lags(input_data_predictors)
+
+        self.feature_names_ = list(feature_names)
+        self.frequency_ = frequency
+        self.lags_to_time_deltas_ = {
+            key: pd.Timedelta(minutes=val) for key, val in feature_to_lags_in_min
+        }
+
+        # Check that the frequency of the input data matches frequency of the lags
+        if not self._frequency_matches(
+            input_data_predictors.index.drop_duplicates()
+        ):  # Several training horizons give duplicates
+            raise ValueError(
+                f"The input data frequency ({input_data_predictors.index.freq}) does not match the model frequency ({self.frequency})."
+            )
+
+        self.feature_importances_ = np.ones(len(self.feature_names_)) / (
+            len(self.feature_names_) or 1.0
+        )

From f3d95e9fd58463d3b5065b95eb96026726beac54 Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Tue, 9 Dec 2025 09:24:12 +0100
Subject: [PATCH 02/13] wip

---
 .../openstef_models/models/forecasting/median_forecaster.py   | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
index 5ba6554c5..3df4fbfe6 100644
--- a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
+++ b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
@@ -426,6 +426,6 @@ def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None
                 f"The input data frequency ({input_data_predictors.index.freq}) does not match the model frequency ({self.frequency})."
             )
 
-        self.feature_importances_ = np.ones(len(self.feature_names_)) / (
-            len(self.feature_names_) or 1.0
+        self.feature_importances = np.ones(len(self.feature_names)) / (
+            len(self.feature_names) or 1.0
         )

From e33a01efa7bf2aeaf79cb13fda4cc25c36cbf8a3 Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Fri, 19 Dec 2025 10:52:42 +0100
Subject: [PATCH 03/13] Add median model

---
 .gitignore                                    |   4 +
 .../models/forecasting/median_forecaster.py   | 227 +++---------------
 .../presets/forecasting_workflow.py           |  19 +-
 3 files changed, 61 insertions(+), 189 deletions(-)

diff --git a/.gitignore b/.gitignore
index c66abc3dc..47f061709 100644
--- a/.gitignore
+++ b/.gitignore
@@ -125,3 +125,7 @@ certificates/
 
 # Benchmark outputs
 benchmark_results*/
+
+# Mlflow
+mlflow
+mlflow_artifacts_local
\ No newline at end of file
diff --git a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
index 3df4fbfe6..a55e4f871 100644
--- a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
+++ b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
@@ -1,3 +1,4 @@
+# ruff: noqa
 # SPDX-FileCopyrightText: 2025 Contributors to the OpenSTEF project <short.term.energy.forecasts@alliander.com>
 #
 # SPDX-License-Identifier: MPL-2.0
@@ -34,18 +35,18 @@
 """
 
 from datetime import timedelta
-from typing import Self, override
+from typing import override
 
 import numpy as np
 import pandas as pd
 from pydantic import Field
 
 from openstef_core.datasets.validated_datasets import ForecastDataset, ForecastInputDataset
-from openstef_core.exceptions import ModelLoadingError
-from openstef_core.mixins import State
+
 from openstef_core.mixins.predictor import HyperParams
 from openstef_core.types import LeadTime, Quantile
 from openstef_models.explainability.mixins import ExplainableForecaster
+from openstef_core.utils.pydantic import timedelta_from_isoformat
 from openstef_models.models.forecasting.forecaster import Forecaster, ForecasterConfig
 
 
@@ -89,6 +90,7 @@ class MedianForecasterConfig(ForecasterConfig):
         default_factory=MedianForecasterHyperParams,
     )
 
+
 class MedianForecaster(Forecaster, ExplainableForecaster):
     """Median forecaster using lag features for predictions.
 
@@ -101,9 +103,12 @@ class MedianForecaster(Forecaster, ExplainableForecaster):
         context_window: Time delta representing the context window size for input data.
                         This defines how much historical data is considered for making predictions.
 
-                _config: BaseCaseForecasterConfig
+                _config: MedianForecasterConfig
     """
 
+    Config = MedianForecasterConfig
+    HyperParams = MedianForecasterHyperParams
+
     def __init__(
         self,
         config: MedianForecasterConfig,
@@ -115,6 +120,8 @@ def __init__(
                    If None, uses default configuration with 7-day primary and 14-day fallback lags.
         """
         self._config = config
+        self.is_fitted_ = False
+        self.feature_names_: list[str] = []
 
     @property
     @override
@@ -126,148 +133,21 @@ def config(self) -> MedianForecasterConfig:
     def hyperparams(self) -> MedianForecasterHyperParams:
         return self._config.hyperparams
 
-    @override
-    def to_state(self) -> State:
-        return {
-            "version": MODEL_CODE_VERSION,
-            "config": self.config.model_dump(mode="json"),
-        }
-
-    @override
-    def from_state(self, state: State) -> Self:
-        if not isinstance(state, dict) or "version" not in state or state["version"] > MODEL_CODE_VERSION:
-            raise ModelLoadingError("Invalid state for BaseCaseForecaster")
-
-        return self.__class__(config=MedianForecasterConfig.model_validate(state["config"]))
-
     @property
     @override
     def is_fitted(self) -> bool:
-        return self._xgboost_model.__sklearn_is_fitted__() 
-
-
+        return self.is_fitted_
 
     @property
     @override
     def feature_importances(self) -> pd.DataFrame:
-        booster = self._xgboost_model.get_booster()
-        weights_df = pd.DataFrame(
-            data=booster.get_score(importance_type="gain"),
-            index=[quantile.format() for quantile in self.config.quantiles],
-        ).transpose()
-        weights_df.index.name = "feature_name"
-        weights_df.columns.name = "quantiles"
-
-        weights_abs = weights_df.abs()
-        total = weights_abs.sum(axis=0).replace(to_replace=0, value=1.0)  # pyright: ignore[reportUnknownMemberType]
-
-        return weights_abs / total
-
-
 
-    @staticmethod
-    def _infer_frequency(index: pd.DatetimeIndex) -> pd.Timedelta:
-        """
-        Infer the frequency of a pandas DatetimeIndex if the freq attribute is not set.
-        This method calculates the most common time difference between consecutive timestamps,
-        which is more permissive of missing chunks of data than the pandas infer_freq method.
-
-        Args:
-            index (pd.DatetimeIndex): The datetime index to infer the frequency from.
-
-        Returns:
-            pd.Timedelta: The inferred frequency as a pandas Timedelta.
-        """
-        if len(index) < 2:
-            raise ValueError(
-                "Cannot infer frequency from an index with fewer than 2 timestamps."
-            )
-
-        # Calculate the differences between consecutive timestamps
-        deltas = index.to_series().diff().dropna()
-
-        # Find the most common difference
-        inferred_freq = deltas.mode().iloc[0]
-        return inferred_freq
-
-    def _frequency_matches(self, index: pd.DatetimeIndex) -> bool:
-        """
-        Check if the frequency of the input data matches the model frequency.
-
-        Args:
-            x (pd.DataFrame): The input data to check.
-
-        Returns:
-            bool: True if the frequencies match, False otherwise.
-        """
-        if not isinstance(index, pd.DatetimeIndex):
-            raise ValueError(
-                "The index of the input data must be a pandas DatetimeIndex."
-            )
-
-        if index.freq is None:
-            input_frequency = self._infer_frequency(index)
-        else:
-            input_frequency = index.freq
-
-        return input_frequency == pd.Timedelta(minutes=self.frequency)
-
-    @staticmethod
-    def _extract_and_validate_lags(
-        x: pd.DataFrame,
-    ) -> tuple[tuple[str], int, list[tuple[str, int]]]:
-        """Extract and validate the lag features from the input data.
-
-        This method checks that the lag features are evenly spaced and match the frequency of the input data.
-        It also extracts the lag features and their corresponding time deltas.
-        Args:
-            x (pd.DataFrame): The input data containing lag features.
-        Returns:
-            tuple: A tuple containing:
-                - A list of feature names, sorted by their lag in minutes.
-                - The frequency of the lag features in minutes.
-                - A list of tuples containing the lag feature names and their corresponding time deltas in minutes.
-        """
-        # Check that the input data contains the required lag features
-        feature_names = list(x.columns[x.columns.str.startswith("T-")])
-        if len(feature_names) == 0:
-            raise ValueError("No lag features found in the input data.")
-
-        # Convert all lags to minutes to make comparable
-        feature_to_lags_in_min = []
-        for feature in feature_names:
-            if feature.endswith("min"):
-                lag_in_min = int(feature.split("-")[1].split("min")[0])
-            elif feature.endswith("d"):
-                lag_in_min = int(feature.split("-")[1].split("d")[0]) * 60 * 24
-            else:
-                raise ValueError(
-                    f"Feature name '{feature}' does not follow the expected format."
-                    " Expected format is 'T-<lag_in_minutes>' or 'T-<lag_in_days>d'."
-                )
-            feature_to_lags_in_min.append((feature, lag_in_min))
-
-        # Sort the features by lag in minutes
-        feature_to_lags_in_min.sort(key=lambda x: x[1])
-        sorted_features, sorted_lags_in_min = zip(*feature_to_lags_in_min)
-
-        # Check that the lags are evenly spaced
-        diffs = np.diff(sorted_lags_in_min)
-        unique_diffs = np.unique(diffs)
-        if len(unique_diffs) > 1:
-            raise ValueError(
-                "Lag features are not evenly spaced. "
-                f"Got lags with differences: {unique_diffs} min. "
-                "Please ensure that the lag features are generated correctly."
-            )
-        frequency = unique_diffs[0]
-
-        return sorted_features, frequency, feature_to_lags_in_min
+        return pd.DataFrame(
+            data=self.feature_importances_, columns=[self.config.quantiles[0].format()], index=self.feature_names_
+        )
 
     @staticmethod
-    def _fill_diagonal_with_median(
-        lag_array: np.ndarray, start: int, end: int, median: float
-    ) -> np.ndarray | None:
+    def _fill_diagonal_with_median(lag_array: np.ndarray, start: int, end: int, median: float) -> np.ndarray | None:
         # Use the calculated median to fill in future lag values where this prediction would be used as input.
 
         # If the start index is beyond the array bounds, no future updates are needed from this step.
@@ -310,17 +190,10 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
 
         input_data: pd.DataFrame = data.input_data(start=data.forecast_start)
 
-        if not self._frequency_matches(input_data.index):
-            raise ValueError(
-                f"The input data frequency ({input_data.index.freq}) does not match the model frequency ({self.frequency})."
-            )
-
         # Check that the input data contains the required lag features
-        missing_features = set(self.feature_names) - set(input_data.columns)
+        missing_features = set(self.feature_names_) - set(data.feature_names)
         if missing_features:
-            raise ValueError(
-                f"The input data is missing the following lag features: {missing_features}"
-            )
+            raise ValueError(f"The input data is missing the following lag features: {missing_features}")
 
         # Reindex the input data to ensure there are no gaps in the time series.
         # This is important for the autoregressive logic that follows.
@@ -328,14 +201,14 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         original_index = input_data.index.copy()
         first_index = input_data.index[0]
         last_index = input_data.index[-1]
-        freq = pd.Timedelta(minutes=self.frequency)
+        freq = self.frequency_
         # Create a new date range with the expected frequency.
         new_index = pd.date_range(first_index, last_index, freq=freq)
         # Reindex the input DataFrame, filling any new timestamps with NaN.
         input_data = input_data.reindex(new_index, fill_value=np.nan)
 
         # Select only the lag feature columns in the specified order.
-        lag_df = input_data[self.feature_names]
+        lag_df = input_data[self.feature_names_]
 
         # Convert the lag DataFrame and its index to NumPy arrays for faster processing.
         lag_array = lag_df.to_numpy()
@@ -344,14 +217,10 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         prediction = np.full(lag_array.shape[0], np.nan)
 
         # Calculate the time step size based on the model frequency.
-        step_size = pd.Timedelta(minutes=self.frequency)
+        step_size = self.frequency_
         # Determine the number of steps corresponding to the smallest and largest lags.
-        smallest_lag_steps = int(
-            self.lags_to_time_deltas_[self.feature_names[0]] / step_size
-        )
-        largest_lag_steps = int(
-            self.lags_to_time_deltas_[self.feature_names[-1]] / step_size
-        )
+        smallest_lag_steps = int(self.lags_to_time_deltas_[self.feature_names_[0]] / step_size)
+        largest_lag_steps = int(self.lags_to_time_deltas_[self.feature_names_[-1]] / step_size)
 
         # Iterate through each time step in the reindexed data.
         for time_step in range(lag_array.shape[0]):
@@ -377,23 +246,16 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         # Convert the prediction array back to a pandas DataFrame using the reindexed time index.
         prediction_df = pd.DataFrame(prediction, index=time_index, columns=["median"])
         # Select only the predictions corresponding to the original input index.
-        prediction = prediction_df.loc[original_index].to_numpy().flatten()
-
-        # Return the final predictions as a ForecastDataset.
-        predictions = pd.DataFrame(
-            data=prediction,
-            index=input_data.index,
-            columns=[quantile.format() for quantile in self.config.quantiles],
-        )
+        # prediction = prediction_df.loc[original_index].to_numpy().flatten()
 
         return ForecastDataset(
-            data=predictions,
+            data=prediction_df.dropna().rename(columns={"median": self.config.quantiles[0].format()}),
             sample_interval=data.sample_interval,
         )
 
     @override
     def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None = None) -> None:
-                """This model does not have any hyperparameters to fit,
+        """This model does not have any hyperparameters to fit,
         but it does need to know the feature names of the lag features and the order of these.
 
         Lag features are expected to be evently spaced and match the frequency of the input data.
@@ -401,31 +263,20 @@ def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None
         For example, T-1min, T-2min, T-3min or T-1d, T-2d.
 
         Which lag features are used is determined by the feature engineering step.
-        """                
-        input_data_predictors = data.input_data()
-    
-        target_series = data.target_series
-    
-        (
-            feature_names,
-            frequency,
-            feature_to_lags_in_min,
-        ) = self._extract_and_validate_lags(input_data_predictors)
-
-        self.feature_names_ = list(feature_names)
-        self.frequency_ = frequency
+        """
+        lag_perfix = f"{data.target_column}_lag_"
+        self.feature_names_ = [
+            feature_name for feature_name in data.feature_names if feature_name.startswith(lag_perfix)
+        ]
+
         self.lags_to_time_deltas_ = {
-            key: pd.Timedelta(minutes=val) for key, val in feature_to_lags_in_min
+            feature_name: timedelta_from_isoformat(feature_name.replace(lag_perfix, ""))
+            for feature_name in self.feature_names_
         }
 
-        # Check that the frequency of the input data matches frequency of the lags
-        if not self._frequency_matches(
-            input_data_predictors.index.drop_duplicates()
-        ):  # Several training horizons give duplicates
-            raise ValueError(
-                f"The input data frequency ({input_data_predictors.index.freq}) does not match the model frequency ({self.frequency})."
-            )
+        self.frequency_ = data.sample_interval
 
-        self.feature_importances = np.ones(len(self.feature_names)) / (
-            len(self.feature_names) or 1.0
-        )
+        self.feature_names_ = sorted(self.feature_names_, key=lambda f: self.lags_to_time_deltas_[f])
+
+        self.feature_importances_ = np.ones(len(self.feature_names_)) / (len(self.feature_names_) or 1.0)
+        self.is_fitted_ = True
diff --git a/packages/openstef-models/src/openstef_models/presets/forecasting_workflow.py b/packages/openstef-models/src/openstef_models/presets/forecasting_workflow.py
index a433ca692..2de522083 100644
--- a/packages/openstef-models/src/openstef_models/presets/forecasting_workflow.py
+++ b/packages/openstef-models/src/openstef_models/presets/forecasting_workflow.py
@@ -30,6 +30,7 @@
 from openstef_models.models import ForecastingModel
 from openstef_models.models.forecasting.flatliner_forecaster import FlatlinerForecaster
 from openstef_models.models.forecasting.gblinear_forecaster import GBLinearForecaster
+from openstef_models.models.forecasting.median_forecaster import MedianForecaster
 from openstef_models.models.forecasting.xgboost_forecaster import XGBoostForecaster
 from openstef_models.transforms.energy_domain import WindPowerFeatureAdder
 from openstef_models.transforms.general import (
@@ -99,7 +100,7 @@ class ForecastingWorkflowConfig(BaseConfig):  # PredictionJob
     run_name: str | None = Field(default=None, description="Optional name for this workflow run.")
 
     # Model configuration
-    model: Literal["xgboost", "gblinear", "flatliner"] = Field(
+    model: Literal["xgboost", "gblinear", "flatliner", "median"] = Field(
         description="Type of forecasting model to use."
     )  # TODO(#652): Implement median forecaster
     quantiles: list[Quantile] = Field(
@@ -373,6 +374,22 @@ def create_forecasting_workflow(config: ForecastingWorkflowConfig) -> CustomFore
                 add_quantiles_from_std=False,
             ),
         ]
+    elif config.model == "median":
+        preprocessing = [
+            LagsAdder(
+                history_available=config.predict_history,
+                horizons=config.horizons,
+                add_trivial_lags=True,
+                target_column=config.target_column,
+            )
+        ]
+        forecaster = MedianForecaster(
+            config=MedianForecaster.Config(
+                quantiles=config.quantiles,
+                horizons=config.horizons,
+            ),
+        )
+        postprocessing = []
     elif config.model == "flatliner":
         preprocessing = []
         forecaster = FlatlinerForecaster(

From 081eb4fbaf704939bc252c579892d062a76ffca2 Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Fri, 19 Dec 2025 13:23:06 +0100
Subject: [PATCH 04/13] Fix linting

---
 .../models/forecasting/median_forecaster.py   | 33 +++++++++----------
 1 file changed, 15 insertions(+), 18 deletions(-)

diff --git a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
index a55e4f871..19feeeb7f 100644
--- a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
+++ b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
@@ -1,4 +1,3 @@
-# ruff: noqa
 # SPDX-FileCopyrightText: 2025 Contributors to the OpenSTEF project <short.term.energy.forecasts@alliander.com>
 #
 # SPDX-License-Identifier: MPL-2.0
@@ -42,11 +41,10 @@
 from pydantic import Field
 
 from openstef_core.datasets.validated_datasets import ForecastDataset, ForecastInputDataset
-
 from openstef_core.mixins.predictor import HyperParams
 from openstef_core.types import LeadTime, Quantile
-from openstef_models.explainability.mixins import ExplainableForecaster
 from openstef_core.utils.pydantic import timedelta_from_isoformat
+from openstef_models.explainability.mixins import ExplainableForecaster
 from openstef_models.models.forecasting.forecaster import Forecaster, ForecasterConfig
 
 
@@ -173,37 +171,37 @@ def _fill_diagonal_with_median(lag_array: np.ndarray, start: int, end: int, medi
             updated_diagonal = np.diag(view).copy()
             updated_diagonal[diagonal_nan_mask] = median
             np.fill_diagonal(view, updated_diagonal)
+        return None
 
     @override
     def predict(self, data: ForecastInputDataset) -> ForecastDataset:
-        """
-        Predict the median of the lag features for each time step in the context window.
+        """Predict the median of the lag features for each time step in the context window.
 
         Args:
-            x (pd.DataFrame): The input data for prediction. This should be a pandas dataframe with lag features.
+            data (ForecastInputDataset): The input data for prediction.
+            This should be a pandas dataframe with lag features.
 
         Returns:
             np.array: The predicted median for each time step in the context window.
             If any lag feature is NaN, this will be ignored.
             If all lag features are NaN, the regressor will return NaN.
-        """
 
+        Raises:
+            ValueError: If the input data is missing any of the required lag features.
+        """
         input_data: pd.DataFrame = data.input_data(start=data.forecast_start)
 
         # Check that the input data contains the required lag features
         missing_features = set(self.feature_names_) - set(data.feature_names)
         if missing_features:
-            raise ValueError(f"The input data is missing the following lag features: {missing_features}")
+            msg = f"The input data is missing the following lag features: {missing_features}"
+            raise ValueError(msg)
 
         # Reindex the input data to ensure there are no gaps in the time series.
         # This is important for the autoregressive logic that follows.
         # Store the original index to return predictions aligned with the input.
-        original_index = input_data.index.copy()
-        first_index = input_data.index[0]
-        last_index = input_data.index[-1]
-        freq = self.frequency_
         # Create a new date range with the expected frequency.
-        new_index = pd.date_range(first_index, last_index, freq=freq)
+        new_index = pd.date_range(input_data.index[0], input_data.index[-1], freq=self.frequency_)
         # Reindex the input DataFrame, filling any new timestamps with NaN.
         input_data = input_data.reindex(new_index, fill_value=np.nan)
 
@@ -212,7 +210,6 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
 
         # Convert the lag DataFrame and its index to NumPy arrays for faster processing.
         lag_array = lag_df.to_numpy()
-        time_index = lag_df.index.to_numpy()
         # Initialize the prediction array with NaNs.
         prediction = np.full(lag_array.shape[0], np.nan)
 
@@ -244,9 +241,7 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
             self._fill_diagonal_with_median(lag_array, start, end, median)
 
         # Convert the prediction array back to a pandas DataFrame using the reindexed time index.
-        prediction_df = pd.DataFrame(prediction, index=time_index, columns=["median"])
-        # Select only the predictions corresponding to the original input index.
-        # prediction = prediction_df.loc[original_index].to_numpy().flatten()
+        prediction_df = pd.DataFrame(prediction, index=lag_df.index.to_numpy(), columns=["median"])
 
         return ForecastDataset(
             data=prediction_df.dropna().rename(columns={"median": self.config.quantiles[0].format()}),
@@ -255,7 +250,9 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
 
     @override
     def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None = None) -> None:
-        """This model does not have any hyperparameters to fit,
+        """Take car of fitting the median forecaster.
+
+        This regressor does not need any fitting,
         but it does need to know the feature names of the lag features and the order of these.
 
         Lag features are expected to be evently spaced and match the frequency of the input data.

From 4dff91a136d6d82a5c854bf153d6cd30be04c168 Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Wed, 7 Jan 2026 16:40:25 +0100
Subject: [PATCH 05/13] added and fixed unit tests

---
 .../models/forecasting/median_forecaster.py   | 108 ++++-
 .../forecasting/test_median_forecaster.py     | 453 ++++++++++++++++++
 2 files changed, 556 insertions(+), 5 deletions(-)
 create mode 100644 packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py

diff --git a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
index 19feeeb7f..2801b4639 100644
--- a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
+++ b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
@@ -143,6 +143,16 @@ def feature_importances(self) -> pd.DataFrame:
         return pd.DataFrame(
             data=self.feature_importances_, columns=[self.config.quantiles[0].format()], index=self.feature_names_
         )
+    
+    @property
+    def frequency(self) -> int:
+        """Retrieve the model input frequency.
+
+        Returns:
+            The frequency of the model input
+
+        """
+        return self.frequency_
 
     @staticmethod
     def _fill_diagonal_with_median(lag_array: np.ndarray, start: int, end: int, median: float) -> np.ndarray | None:
@@ -172,6 +182,53 @@ def _fill_diagonal_with_median(lag_array: np.ndarray, start: int, end: int, medi
             updated_diagonal[diagonal_nan_mask] = median
             np.fill_diagonal(view, updated_diagonal)
         return None
+    
+    @staticmethod
+    def _infer_frequency(index: pd.DatetimeIndex) -> pd.Timedelta:
+        """
+        Infer the frequency of a pandas DatetimeIndex if the freq attribute is not set.
+        This method calculates the most common time difference between consecutive timestamps,
+        which is more permissive of missing chunks of data than the pandas infer_freq method.
+
+        Args:
+            index (pd.DatetimeIndex): The datetime index to infer the frequency from.
+
+        Returns:
+            pd.Timedelta: The inferred frequency as a pandas Timedelta.
+        """
+        if len(index) < 2:
+            raise ValueError(
+                "Cannot infer frequency from an index with fewer than 2 timestamps."
+            )
+
+        # Calculate the differences between consecutive timestamps
+        deltas = index.to_series().diff().dropna()
+
+        # Find the most common difference
+        inferred_freq = deltas.mode().iloc[0]
+        return inferred_freq
+    
+    def _frequency_matches(self, index: pd.DatetimeIndex) -> bool:
+        """
+        Check if the frequency of the input data matches the model frequency.
+
+        Args:
+            index (pd.DatetimeIndex): The input data to check.
+
+        Returns:
+            bool: True if the frequencies match, False otherwise.
+        """
+        if not isinstance(index, pd.DatetimeIndex):
+            raise ValueError(
+                "The index of the input data must be a pandas DatetimeIndex."
+            )
+
+        if index.freq is None:
+            input_frequency = self._infer_frequency(index)
+        else:
+            input_frequency = index.freq
+
+        return input_frequency == self.frequency
 
     @override
     def predict(self, data: ForecastInputDataset) -> ForecastDataset:
@@ -189,6 +246,11 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         Raises:
             ValueError: If the input data is missing any of the required lag features.
         """
+
+        if not self.is_fitted:
+            msg = "This MedianForecaster instance is not fitted yet"
+            raise AttributeError(msg)
+
         input_data: pd.DataFrame = data.input_data(start=data.forecast_start)
 
         # Check that the input data contains the required lag features
@@ -196,10 +258,20 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         if missing_features:
             msg = f"The input data is missing the following lag features: {missing_features}"
             raise ValueError(msg)
+        
+        if not self._frequency_matches(data.input_data().index):
+            msg = (
+                f"The frequency of the input data does not match the model frequency. "
+                f"Input data frequency: {data.input_data().index.freq}, "
+                f"Model frequency: {pd.Timedelta(minutes=self.frequency_)}"
+            )
+            raise ValueError(msg)
+        
 
         # Reindex the input data to ensure there are no gaps in the time series.
         # This is important for the autoregressive logic that follows.
         # Store the original index to return predictions aligned with the input.
+        old_index = input_data.index
         # Create a new date range with the expected frequency.
         new_index = pd.date_range(input_data.index[0], input_data.index[-1], freq=self.frequency_)
         # Reindex the input DataFrame, filling any new timestamps with NaN.
@@ -225,13 +297,15 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
             current_lags = lag_array[time_step]
             # Calculate the median of the available lag features, ignoring NaNs.
             median = np.nanmedian(current_lags)
-            # Store the calculated median in the prediction array.
-            prediction[time_step] = median
-
             # If the median calculation resulted in NaN (e.g., all lags were NaN), skip the autoregression step.
-            if np.isnan(median):
+            if not np.isnan(median):
+                median = float(median)
+            else:
                 continue
 
+            # Store the calculated median in the prediction array.
+            prediction[time_step] = median
+            
             # Auto-regressive step: update the lag array for future time steps.
             # Calculate the start and end indices in the future time steps that will be affected.
             start, end = (
@@ -243,9 +317,13 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         # Convert the prediction array back to a pandas DataFrame using the reindexed time index.
         prediction_df = pd.DataFrame(prediction, index=lag_df.index.to_numpy(), columns=["median"])
 
+        # Reindex the prediction DataFrame back to the original input data index.
+        prediction_df = prediction_df.reindex(old_index)
+
         return ForecastDataset(
             data=prediction_df.dropna().rename(columns={"median": self.config.quantiles[0].format()}),
             sample_interval=data.sample_interval,
+            forecast_start=data.forecast_start,
         )
 
     @override
@@ -261,17 +339,37 @@ def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None
 
         Which lag features are used is determined by the feature engineering step.
         """
+        self.frequency_ = data.sample_interval
+        # Check that the frequency of the input data matches frequency of the lags
+        if not self._frequency_matches(
+            data.data.index.drop_duplicates()
+        ):  # Several training horizons give duplicates
+            raise ValueError(
+                f"The input data frequency ({data.data.index.freq}) does not match the model frequency ({self.frequency_})."
+            )
+
         lag_perfix = f"{data.target_column}_lag_"
         self.feature_names_ = [
             feature_name for feature_name in data.feature_names if feature_name.startswith(lag_perfix)
         ]
 
+        if not self.feature_names_:
+            msg = f"No lag features found in the input data with prefix '{lag_perfix}'."
+            raise ValueError(msg)
+
         self.lags_to_time_deltas_ = {
             feature_name: timedelta_from_isoformat(feature_name.replace(lag_perfix, ""))
             for feature_name in self.feature_names_
         }
 
-        self.frequency_ = data.sample_interval
+        # Check if lags are evenly spaced
+        lag_deltas = sorted(self.lags_to_time_deltas_.values())
+        lag_intervals = [
+            (lag_deltas[i] - lag_deltas[i - 1]).total_seconds() for i in range(1, len(lag_deltas))
+        ]
+        if not all(interval == lag_intervals[0] for interval in lag_intervals):
+            msg = "Lag features are not evenly spaced. Please ensure lag features are evenly spaced and match the data frequency."
+            raise ValueError(msg)
 
         self.feature_names_ = sorted(self.feature_names_, key=lambda f: self.lags_to_time_deltas_[f])
 
diff --git a/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py b/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
new file mode 100644
index 000000000..a9d66caaa
--- /dev/null
+++ b/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
@@ -0,0 +1,453 @@
+# ruff: noqa
+# SPDX-FileCopyrightText: 2017-2025 Contributors to the OpenSTEF project <openstef@lfenergy.org>
+#
+# SPDX-License-Identifier: MPL-2.0
+
+import numpy as np
+from openstef_models.models.forecasting.forecaster import ForecastDataset
+from openstef_models.utils.data_split import DataSplitter
+import pandas as pd
+import pytest
+from openstef_models.models.forecasting.median_forecaster import MedianForecaster, MedianForecasterConfig
+from openstef_core.datasets import ForecastInputDataset
+from openstef_core.types import LeadTime, Q
+from openstef_core.testing import create_timeseries_dataset
+from openstef_models.presets import ForecastingWorkflowConfig, create_forecasting_workflow
+
+
+def test_median_returns_median():
+    # Arrange     
+    index = pd.date_range("2020-01-01T00:00", periods=3, freq="h")
+    training_data = create_timeseries_dataset(
+        index=index,
+        load=[1.0,4.0,7.0],
+        load_lag_PT1H=[1.0,np.nan, np.nan],
+        load_lag_PT2H=[4.0,1.0, np.nan],
+        load_lag_PT3H=[7.0,4.0,1.0],
+        available_at=index,
+    )   
+
+    training_input_data = ForecastInputDataset.from_timeseries(
+            dataset=training_data,
+            target_column="load",
+            forecast_start=index[0],
+        )
+
+    expected_result = ForecastDataset(
+            data=pd.DataFrame(
+                {
+                    "quantile_P50": [4.0, 4.0, 4.0],
+                },
+                index=index,
+            ),
+            sample_interval=training_input_data.sample_interval,
+        )
+    expected_result.index.freq = None
+
+    config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT36H")])
+    model = MedianForecaster(config=config)
+
+
+    # Act
+    model.fit(training_input_data)
+    result = model.predict(training_input_data)
+    
+    # Assert
+    assert result.sample_interval == expected_result.sample_interval
+    pd.testing.assert_frame_equal(result.data, expected_result.data)
+
+
+def test_median_handles_some_missing_data():
+    # Arrange
+    index = pd.date_range("2023-01-01", periods=3, freq="h")
+    training_data = create_timeseries_dataset(
+        index=index,
+        load=[1.0, 2.0, 3.0],
+        load_lag_PT1H=[1.0, np.nan, np.nan],
+        load_lag_PT2H=[np.nan, 1.0, np.nan],
+        load_lag_PT3H=[3.0, np.nan, 1.0],
+        available_at=index,
+    )
+
+    training_input_data = ForecastInputDataset.from_timeseries(
+        dataset=training_data,
+        target_column="load",
+        forecast_start=index[0],
+    )
+
+    expected_result = ForecastDataset(
+        data=pd.DataFrame(
+            {
+                "quantile_P50": [2.0, 1.5, 1.5],
+            },
+            index=index,
+        ),
+        sample_interval=training_input_data.sample_interval,
+    )
+    expected_result.index.freq = None
+
+    config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT3M")])
+    model = MedianForecaster(config=config)
+
+    # Act
+    model.fit(training_input_data)
+    result = model.predict(training_input_data)
+
+    # Assert
+    assert result.sample_interval == expected_result.sample_interval
+    pd.testing.assert_frame_equal(result.data, expected_result.data)
+
+def test_median_handles_missing_data_for_some_horizons():
+    # Arrange
+    index = pd.date_range("2023-01-01", periods=3, freq="h")
+    training_data = create_timeseries_dataset(
+        index=index,
+        load=[5.0, 5.0, 5.0],
+        load_lag_PT1H=[5.0, np.nan, np.nan],
+        load_lag_PT2H=[np.nan, 5.0, np.nan],
+        load_lag_PT3H=[np.nan, np.nan, 5.0],
+        available_at=index,
+    )
+
+    training_input_data = ForecastInputDataset.from_timeseries(
+        dataset=training_data,
+        target_column="load",
+        forecast_start=index[0],
+    )
+
+    expected_result = ForecastDataset(
+        data=pd.DataFrame(
+            {
+                "quantile_P50": [5.0, 5.0, 5.0],
+            },
+            index=index,
+        ),
+        sample_interval=training_input_data.sample_interval,
+    )
+    expected_result.index.freq = None
+
+    config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT3M")])
+    model = MedianForecaster(config=config)
+
+    # Act
+    model.fit(training_input_data)
+    result = model.predict(training_input_data)
+
+    # Assert
+    assert result.sample_interval == expected_result.sample_interval
+    pd.testing.assert_frame_equal(result.data, expected_result.data)
+
+def test_median_handles_all_missing_data():
+    # Arrange
+    index = pd.date_range("2023-01-01", periods=3, freq="h")
+    training_data = create_timeseries_dataset(
+        index=index,
+        load=[np.nan, np.nan, np.nan],
+        load_lag_PT1H=[np.nan, np.nan, np.nan],
+        load_lag_PT2H=[np.nan, np.nan, np.nan],
+        load_lag_PT3H=[np.nan, np.nan, np.nan],
+        available_at=index,
+    )
+
+    training_input_data = ForecastInputDataset.from_timeseries(
+        dataset=training_data,
+        target_column="load",
+        forecast_start=index[0],
+    )
+
+    expected_result = ForecastDataset(
+        data=pd.DataFrame(
+            {
+                "quantile_P50": [],
+            },
+            index=pd.DatetimeIndex([], freq="h"),
+        ),
+        sample_interval=training_input_data.sample_interval,
+        forecast_start=training_input_data.forecast_start,
+    )
+
+
+    config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT3M")])
+    model = MedianForecaster(config=config)
+
+    # Act
+    model.fit(training_input_data)
+    result = model.predict(training_input_data)
+
+    # Assert
+    assert result.sample_interval == expected_result.sample_interval
+    pd.testing.assert_frame_equal(result.data, expected_result.data)
+
+def test_median_uses_lag_features_if_available():
+    # Arrange
+    index = pd.date_range("2023-01-01T00:00", periods=3, freq="h")
+    training_data = create_timeseries_dataset(
+        index=index,
+        load=[4.0, 5.0, 6.0],
+        load_lag_PT1H=[1.0, 2.0, 3.0],
+        load_lag_PT2H=[4.0, 5.0, 6.0],
+        load_lag_PT3H=[7.0, 8.0, 9.0],
+        available_at=index,
+    )
+
+    training_input_data = ForecastInputDataset.from_timeseries(
+        dataset=training_data,
+        target_column="load",
+        forecast_start=index[0],
+    )
+
+    expected_result = ForecastDataset(
+        data=pd.DataFrame(
+            {
+                "quantile_P50": [4.0, 5.0, 6.0],
+            },
+            index=index,
+        ),
+        sample_interval=training_input_data.sample_interval,
+    )
+    expected_result.index.freq = None
+
+    config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT3M")])
+    model = MedianForecaster(config=config)
+
+    # Act
+    model.fit(training_input_data)
+    result = model.predict(training_input_data)
+
+    # Assert
+    assert result.sample_interval == expected_result.sample_interval
+    pd.testing.assert_frame_equal(result.data, expected_result.data)
+
+def test_median_handles_small_gap():
+    # Arrange
+    index = pd.date_range("2023-01-01T00:00", periods=5, freq="h")
+    training_data = create_timeseries_dataset(
+        index=index,
+        load=[5.0, 4.0,3.0,2.0,1.0,],
+        load_lag_PT1H=[1.0,np.nan, np.nan, np.nan, np.nan],
+        load_lag_PT2H=[2.0,1.0, np.nan, np.nan, np.nan],
+        load_lag_PT3H=[3.0,2.0,1.0, np.nan, np.nan],
+        load_lag_PT4H=[4.0,3.0,2.0,1.0, np.nan],
+        load_lag_PT5H=[5.0,4.0,3.0,2.0,1.0],
+        available_at=index,
+    )
+
+    training_input_data = ForecastInputDataset.from_timeseries(
+        dataset=training_data,
+        target_column="load",
+        forecast_start=index[0],
+    )
+
+    # Remove the second row to create a small gap
+    training_input_data.data = training_input_data.data[training_input_data.data.index != "2023-01-01T01:00"]
+
+
+    expected_result = ForecastDataset(
+        data=pd.DataFrame(
+            {
+                "quantile_P50": [3.0,  3.0, 3.0, 3.0],
+            },
+            index=pd.date_range("2023-01-01T00:00", periods=5, freq="h").delete(1),
+        ),
+        sample_interval=training_input_data.sample_interval,
+    )
+    expected_result.index.freq = None
+
+    config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT3M")])
+    model = MedianForecaster(config=config)
+
+    # Act
+    model.fit(training_input_data)
+    result = model.predict(training_input_data)
+
+    # Assert
+    assert result.sample_interval == expected_result.sample_interval
+    pd.testing.assert_frame_equal(result.data, expected_result.data)
+
+def test_median_handles_large_gap():
+    # Arrange
+    index_1 = pd.date_range("2023-01-01T00:00", periods=3, freq="h")
+    training_data_1 = create_timeseries_dataset(
+        index=index_1,
+        load=[1.0, 2.0, 3.0],
+        load_lag_PT1H=[4.0, 5.0, 6.0],
+        load_lag_PT2H=[7.0, 8.0, 9.0],
+        available_at=index_1,
+    )
+    
+    index_2 = pd.date_range("2023-01-02T01:00", periods=3, freq="h")
+    training_data_2 = create_timeseries_dataset(
+        index=index_2,
+        load=[10.0, 11.0, 12.0],
+        load_lag_PT1H=[13.0, 14.0, 15.0],
+        load_lag_PT2H=[16.0, 17.0, 18.0],
+        available_at=index_2,
+    )
+    
+    training_data = training_data_1
+
+    training_data.data = pd.concat([training_data_1.data, training_data_2.data])
+
+
+    training_input_data = ForecastInputDataset.from_timeseries(
+        dataset=training_data,
+        target_column="load",
+        forecast_start=index_1[0],
+    )
+
+    expected_result = ForecastDataset(
+        data=pd.DataFrame(
+            {
+                "quantile_P50": [5.5, 6.5, 7.5, 14.5, 15.5, 16.5],
+            },
+            index=pd.DatetimeIndex(pd.concat([index_1.to_series(), index_2.to_series()])),
+        ),
+        sample_interval=training_input_data.sample_interval,
+    )
+    expected_result.index.freq = None
+
+    config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT3M")])
+    model = MedianForecaster(config=config)
+
+    # Act
+    model.fit(training_input_data)
+    result = model.predict(training_input_data)
+
+    # Assert
+    assert result.sample_interval == expected_result.sample_interval
+    pd.testing.assert_frame_equal(result.data, expected_result.data)
+
+def test_median_fit_with_missing_features_raises():
+    # Arrange
+    index = pd.date_range("2023-01-01", periods=3, freq="h")
+    training_data = create_timeseries_dataset(
+        index=index,
+        load=[1.0, 2.0, 3.0],
+        load_lag_PT1H=[1.0, 2.0, 3.0],
+        load_lag_PT2H=[4.0, 5.0, 6.0],
+        load_lag_PT3H=[7.0, 8.0, 9.0],
+        available_at=index,
+    )
+
+    training_input_data = ForecastInputDataset.from_timeseries(
+        dataset=training_data,
+        target_column="load",
+        forecast_start=index[0],
+    )
+
+    config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT3H")])
+    model = MedianForecaster(config=config)
+    model.fit(training_input_data)
+
+    # Create prediction data missing one lag feature
+    prediction_data = training_data.data.copy()
+    prediction_data = prediction_data.drop(columns=["load_lag_PT3H"])
+    prediction_input_data = ForecastInputDataset(
+        data=prediction_data,
+        target_column="load",
+        forecast_start=index[0],
+        sample_interval=training_input_data.sample_interval,
+    )
+
+    # Act & Assert
+    with pytest.raises(ValueError, match="The input data is missing the following lag features"):
+        model.predict(prediction_input_data)
+
+def test_median_fit_with_no_lag_features_raises():
+    # Arrange
+    index = pd.date_range("2023-01-01", periods=3, freq="h")
+    training_data = create_timeseries_dataset(
+        index=index,
+        load=[1.0, 2.0, 3.0],
+        unrelated_feature=[1.0, 2.0, 3.0],
+        available_at=index,
+    )
+
+    training_input_data = ForecastInputDataset.from_timeseries(
+        dataset=training_data,
+        target_column="load",
+        forecast_start=index[0],
+    )
+
+    config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT3H")])
+    model = MedianForecaster(config=config)
+
+    # Act & Assert
+    with pytest.raises(ValueError, match="No lag features found in the input data."):
+        model.fit(training_input_data)
+
+def test_median_fit_with_inconsistent_lag_features_raises():
+    # Arrange
+    index = pd.date_range("2023-01-01", periods=3, freq="h")
+    training_data = create_timeseries_dataset(
+        index=index,
+        load=[1.0, 2.0, 3.0],
+        load_lag_PT1H=[1.0, 2.0, 3.0],
+        load_lag_PT5H=[4.0, 5.0, 6.0],
+        load_lag_PT60H=[7.0, 8.0, 9.0],
+        load_lag_P4D=[10.0, 11.0, 12.0],
+        available_at=index,
+    )
+
+    training_input_data = ForecastInputDataset.from_timeseries(
+        dataset=training_data,
+        target_column="load",
+        forecast_start=index[0],
+    )
+
+    config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT3H")])
+    model = MedianForecaster(config=config)
+
+    # Act & Assert
+    with pytest.raises(ValueError, match="Lag features are not evenly spaced"):
+        model.fit(training_input_data)
+
+def test_median_fit_with_inconsistent_frequency_raises():
+    # Arrange
+    index = pd.date_range("2023-01-01", periods=3, freq="min")
+    training_data = create_timeseries_dataset(
+        index=index,
+        load=[1.0, 2.0, 3.0],
+        load_lag_PT1H=[1.0, 2.0, 3.0],
+        load_lag_PT2H=[4.0, 5.0, 6.0],
+        load_lag_PT3H=[7.0, 8.0, 9.0],
+        available_at=index,
+    )
+
+    training_input_data = ForecastInputDataset.from_timeseries(
+        dataset=training_data,
+        target_column="load",
+        forecast_start=index[0],
+    )
+
+    config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT3H")])
+    model = MedianForecaster(config=config)
+
+    # Act & Assert
+    with pytest.raises(ValueError, match="does not match the model frequency."):
+        model.fit(training_input_data)
+
+def test_predicting_without_fitting_raises():
+    # Arrange
+    index = pd.date_range("2023-01-01", periods=3, freq="min")
+    training_data = create_timeseries_dataset(
+        index=index,
+        load=[1.0, 2.0, 3.0],
+        load_lag_PT1M=[1.0, 2.0, 3.0],
+        load_lag_PT2M=[4.0, 5.0, 6.0],
+        load_lag_PT3M=[7.0, 8.0, 9.0],
+        available_at=index,
+    )
+
+    training_input_data = ForecastInputDataset.from_timeseries(
+        dataset=training_data,
+        target_column="load",
+        forecast_start=index[0],
+    )
+
+    config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT3M")])
+    model = MedianForecaster(config=config)
+
+    # Act & Assert
+    with pytest.raises(AttributeError, match="This MedianForecaster instance is not fitted yet"):
+        model.predict(training_input_data)

From b5ba22c48875346826187240feb77337dc00e38d Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Wed, 7 Jan 2026 16:41:33 +0100
Subject: [PATCH 06/13] ran linting

---
 .../models/forecasting/median_forecaster.py    | 18 +++++++-----------
 1 file changed, 7 insertions(+), 11 deletions(-)

diff --git a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
index 2801b4639..e5422a2d9 100644
--- a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
+++ b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
@@ -143,7 +143,7 @@ def feature_importances(self) -> pd.DataFrame:
         return pd.DataFrame(
             data=self.feature_importances_, columns=[self.config.quantiles[0].format()], index=self.feature_names_
         )
-    
+
     @property
     def frequency(self) -> int:
         """Retrieve the model input frequency.
@@ -182,11 +182,10 @@ def _fill_diagonal_with_median(lag_array: np.ndarray, start: int, end: int, medi
             updated_diagonal[diagonal_nan_mask] = median
             np.fill_diagonal(view, updated_diagonal)
         return None
-    
+
     @staticmethod
     def _infer_frequency(index: pd.DatetimeIndex) -> pd.Timedelta:
-        """
-        Infer the frequency of a pandas DatetimeIndex if the freq attribute is not set.
+        """Infer the frequency of a pandas DatetimeIndex if the freq attribute is not set.
         This method calculates the most common time difference between consecutive timestamps,
         which is more permissive of missing chunks of data than the pandas infer_freq method.
 
@@ -207,10 +206,9 @@ def _infer_frequency(index: pd.DatetimeIndex) -> pd.Timedelta:
         # Find the most common difference
         inferred_freq = deltas.mode().iloc[0]
         return inferred_freq
-    
+
     def _frequency_matches(self, index: pd.DatetimeIndex) -> bool:
-        """
-        Check if the frequency of the input data matches the model frequency.
+        """Check if the frequency of the input data matches the model frequency.
 
         Args:
             index (pd.DatetimeIndex): The input data to check.
@@ -246,7 +244,6 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         Raises:
             ValueError: If the input data is missing any of the required lag features.
         """
-
         if not self.is_fitted:
             msg = "This MedianForecaster instance is not fitted yet"
             raise AttributeError(msg)
@@ -258,7 +255,7 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         if missing_features:
             msg = f"The input data is missing the following lag features: {missing_features}"
             raise ValueError(msg)
-        
+
         if not self._frequency_matches(data.input_data().index):
             msg = (
                 f"The frequency of the input data does not match the model frequency. "
@@ -266,7 +263,6 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
                 f"Model frequency: {pd.Timedelta(minutes=self.frequency_)}"
             )
             raise ValueError(msg)
-        
 
         # Reindex the input data to ensure there are no gaps in the time series.
         # This is important for the autoregressive logic that follows.
@@ -305,7 +301,7 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
 
             # Store the calculated median in the prediction array.
             prediction[time_step] = median
-            
+
             # Auto-regressive step: update the lag array for future time steps.
             # Calculate the start and end indices in the future time steps that will be affected.
             start, end = (

From 8aead577eeabb8a7ccd232f4824d5efc99232d06 Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Wed, 7 Jan 2026 16:42:38 +0100
Subject: [PATCH 07/13] ran format

---
 .../models/forecasting/median_forecaster.py   | 16 +---
 .../forecasting/test_median_forecaster.py     | 76 +++++++++++--------
 2 files changed, 48 insertions(+), 44 deletions(-)

diff --git a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
index e5422a2d9..d314fa1f9 100644
--- a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
+++ b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
@@ -196,9 +196,7 @@ def _infer_frequency(index: pd.DatetimeIndex) -> pd.Timedelta:
             pd.Timedelta: The inferred frequency as a pandas Timedelta.
         """
         if len(index) < 2:
-            raise ValueError(
-                "Cannot infer frequency from an index with fewer than 2 timestamps."
-            )
+            raise ValueError("Cannot infer frequency from an index with fewer than 2 timestamps.")
 
         # Calculate the differences between consecutive timestamps
         deltas = index.to_series().diff().dropna()
@@ -217,9 +215,7 @@ def _frequency_matches(self, index: pd.DatetimeIndex) -> bool:
             bool: True if the frequencies match, False otherwise.
         """
         if not isinstance(index, pd.DatetimeIndex):
-            raise ValueError(
-                "The index of the input data must be a pandas DatetimeIndex."
-            )
+            raise ValueError("The index of the input data must be a pandas DatetimeIndex.")
 
         if index.freq is None:
             input_frequency = self._infer_frequency(index)
@@ -337,9 +333,7 @@ def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None
         """
         self.frequency_ = data.sample_interval
         # Check that the frequency of the input data matches frequency of the lags
-        if not self._frequency_matches(
-            data.data.index.drop_duplicates()
-        ):  # Several training horizons give duplicates
+        if not self._frequency_matches(data.data.index.drop_duplicates()):  # Several training horizons give duplicates
             raise ValueError(
                 f"The input data frequency ({data.data.index.freq}) does not match the model frequency ({self.frequency_})."
             )
@@ -360,9 +354,7 @@ def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None
 
         # Check if lags are evenly spaced
         lag_deltas = sorted(self.lags_to_time_deltas_.values())
-        lag_intervals = [
-            (lag_deltas[i] - lag_deltas[i - 1]).total_seconds() for i in range(1, len(lag_deltas))
-        ]
+        lag_intervals = [(lag_deltas[i] - lag_deltas[i - 1]).total_seconds() for i in range(1, len(lag_deltas))]
         if not all(interval == lag_intervals[0] for interval in lag_intervals):
             msg = "Lag features are not evenly spaced. Please ensure lag features are evenly spaced and match the data frequency."
             raise ValueError(msg)
diff --git a/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py b/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
index a9d66caaa..22eb5bb9b 100644
--- a/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
+++ b/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
@@ -16,42 +16,41 @@
 
 
 def test_median_returns_median():
-    # Arrange     
+    # Arrange
     index = pd.date_range("2020-01-01T00:00", periods=3, freq="h")
     training_data = create_timeseries_dataset(
         index=index,
-        load=[1.0,4.0,7.0],
-        load_lag_PT1H=[1.0,np.nan, np.nan],
-        load_lag_PT2H=[4.0,1.0, np.nan],
-        load_lag_PT3H=[7.0,4.0,1.0],
+        load=[1.0, 4.0, 7.0],
+        load_lag_PT1H=[1.0, np.nan, np.nan],
+        load_lag_PT2H=[4.0, 1.0, np.nan],
+        load_lag_PT3H=[7.0, 4.0, 1.0],
         available_at=index,
-    )   
+    )
 
     training_input_data = ForecastInputDataset.from_timeseries(
-            dataset=training_data,
-            target_column="load",
-            forecast_start=index[0],
-        )
+        dataset=training_data,
+        target_column="load",
+        forecast_start=index[0],
+    )
 
     expected_result = ForecastDataset(
-            data=pd.DataFrame(
-                {
-                    "quantile_P50": [4.0, 4.0, 4.0],
-                },
-                index=index,
-            ),
-            sample_interval=training_input_data.sample_interval,
-        )
+        data=pd.DataFrame(
+            {
+                "quantile_P50": [4.0, 4.0, 4.0],
+            },
+            index=index,
+        ),
+        sample_interval=training_input_data.sample_interval,
+    )
     expected_result.index.freq = None
 
     config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT36H")])
     model = MedianForecaster(config=config)
 
-
     # Act
     model.fit(training_input_data)
     result = model.predict(training_input_data)
-    
+
     # Assert
     assert result.sample_interval == expected_result.sample_interval
     pd.testing.assert_frame_equal(result.data, expected_result.data)
@@ -97,6 +96,7 @@ def test_median_handles_some_missing_data():
     assert result.sample_interval == expected_result.sample_interval
     pd.testing.assert_frame_equal(result.data, expected_result.data)
 
+
 def test_median_handles_missing_data_for_some_horizons():
     # Arrange
     index = pd.date_range("2023-01-01", periods=3, freq="h")
@@ -137,6 +137,7 @@ def test_median_handles_missing_data_for_some_horizons():
     assert result.sample_interval == expected_result.sample_interval
     pd.testing.assert_frame_equal(result.data, expected_result.data)
 
+
 def test_median_handles_all_missing_data():
     # Arrange
     index = pd.date_range("2023-01-01", periods=3, freq="h")
@@ -166,7 +167,6 @@ def test_median_handles_all_missing_data():
         forecast_start=training_input_data.forecast_start,
     )
 
-
     config = MedianForecasterConfig(quantiles=[Q(0.5)], horizons=[LeadTime.from_string("PT3M")])
     model = MedianForecaster(config=config)
 
@@ -178,6 +178,7 @@ def test_median_handles_all_missing_data():
     assert result.sample_interval == expected_result.sample_interval
     pd.testing.assert_frame_equal(result.data, expected_result.data)
 
+
 def test_median_uses_lag_features_if_available():
     # Arrange
     index = pd.date_range("2023-01-01T00:00", periods=3, freq="h")
@@ -218,17 +219,24 @@ def test_median_uses_lag_features_if_available():
     assert result.sample_interval == expected_result.sample_interval
     pd.testing.assert_frame_equal(result.data, expected_result.data)
 
+
 def test_median_handles_small_gap():
     # Arrange
     index = pd.date_range("2023-01-01T00:00", periods=5, freq="h")
     training_data = create_timeseries_dataset(
         index=index,
-        load=[5.0, 4.0,3.0,2.0,1.0,],
-        load_lag_PT1H=[1.0,np.nan, np.nan, np.nan, np.nan],
-        load_lag_PT2H=[2.0,1.0, np.nan, np.nan, np.nan],
-        load_lag_PT3H=[3.0,2.0,1.0, np.nan, np.nan],
-        load_lag_PT4H=[4.0,3.0,2.0,1.0, np.nan],
-        load_lag_PT5H=[5.0,4.0,3.0,2.0,1.0],
+        load=[
+            5.0,
+            4.0,
+            3.0,
+            2.0,
+            1.0,
+        ],
+        load_lag_PT1H=[1.0, np.nan, np.nan, np.nan, np.nan],
+        load_lag_PT2H=[2.0, 1.0, np.nan, np.nan, np.nan],
+        load_lag_PT3H=[3.0, 2.0, 1.0, np.nan, np.nan],
+        load_lag_PT4H=[4.0, 3.0, 2.0, 1.0, np.nan],
+        load_lag_PT5H=[5.0, 4.0, 3.0, 2.0, 1.0],
         available_at=index,
     )
 
@@ -241,11 +249,10 @@ def test_median_handles_small_gap():
     # Remove the second row to create a small gap
     training_input_data.data = training_input_data.data[training_input_data.data.index != "2023-01-01T01:00"]
 
-
     expected_result = ForecastDataset(
         data=pd.DataFrame(
             {
-                "quantile_P50": [3.0,  3.0, 3.0, 3.0],
+                "quantile_P50": [3.0, 3.0, 3.0, 3.0],
             },
             index=pd.date_range("2023-01-01T00:00", periods=5, freq="h").delete(1),
         ),
@@ -264,6 +271,7 @@ def test_median_handles_small_gap():
     assert result.sample_interval == expected_result.sample_interval
     pd.testing.assert_frame_equal(result.data, expected_result.data)
 
+
 def test_median_handles_large_gap():
     # Arrange
     index_1 = pd.date_range("2023-01-01T00:00", periods=3, freq="h")
@@ -274,7 +282,7 @@ def test_median_handles_large_gap():
         load_lag_PT2H=[7.0, 8.0, 9.0],
         available_at=index_1,
     )
-    
+
     index_2 = pd.date_range("2023-01-02T01:00", periods=3, freq="h")
     training_data_2 = create_timeseries_dataset(
         index=index_2,
@@ -283,12 +291,11 @@ def test_median_handles_large_gap():
         load_lag_PT2H=[16.0, 17.0, 18.0],
         available_at=index_2,
     )
-    
+
     training_data = training_data_1
 
     training_data.data = pd.concat([training_data_1.data, training_data_2.data])
 
-
     training_input_data = ForecastInputDataset.from_timeseries(
         dataset=training_data,
         target_column="load",
@@ -317,6 +324,7 @@ def test_median_handles_large_gap():
     assert result.sample_interval == expected_result.sample_interval
     pd.testing.assert_frame_equal(result.data, expected_result.data)
 
+
 def test_median_fit_with_missing_features_raises():
     # Arrange
     index = pd.date_range("2023-01-01", periods=3, freq="h")
@@ -353,6 +361,7 @@ def test_median_fit_with_missing_features_raises():
     with pytest.raises(ValueError, match="The input data is missing the following lag features"):
         model.predict(prediction_input_data)
 
+
 def test_median_fit_with_no_lag_features_raises():
     # Arrange
     index = pd.date_range("2023-01-01", periods=3, freq="h")
@@ -376,6 +385,7 @@ def test_median_fit_with_no_lag_features_raises():
     with pytest.raises(ValueError, match="No lag features found in the input data."):
         model.fit(training_input_data)
 
+
 def test_median_fit_with_inconsistent_lag_features_raises():
     # Arrange
     index = pd.date_range("2023-01-01", periods=3, freq="h")
@@ -402,6 +412,7 @@ def test_median_fit_with_inconsistent_lag_features_raises():
     with pytest.raises(ValueError, match="Lag features are not evenly spaced"):
         model.fit(training_input_data)
 
+
 def test_median_fit_with_inconsistent_frequency_raises():
     # Arrange
     index = pd.date_range("2023-01-01", periods=3, freq="min")
@@ -427,6 +438,7 @@ def test_median_fit_with_inconsistent_frequency_raises():
     with pytest.raises(ValueError, match="does not match the model frequency."):
         model.fit(training_input_data)
 
+
 def test_predicting_without_fitting_raises():
     # Arrange
     index = pd.date_range("2023-01-01", periods=3, freq="min")

From be94422ced7bf0c6aebbca0ca8feddd26a72b523 Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Wed, 7 Jan 2026 16:49:13 +0100
Subject: [PATCH 08/13] fix linting and formatting

---
 .../models/forecasting/median_forecaster.py   | 31 ++++++++++++-------
 1 file changed, 20 insertions(+), 11 deletions(-)

diff --git a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
index d314fa1f9..3d42983a6 100644
--- a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
+++ b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
@@ -186,6 +186,7 @@ def _fill_diagonal_with_median(lag_array: np.ndarray, start: int, end: int, medi
     @staticmethod
     def _infer_frequency(index: pd.DatetimeIndex) -> pd.Timedelta:
         """Infer the frequency of a pandas DatetimeIndex if the freq attribute is not set.
+
         This method calculates the most common time difference between consecutive timestamps,
         which is more permissive of missing chunks of data than the pandas infer_freq method.
 
@@ -194,16 +195,19 @@ def _infer_frequency(index: pd.DatetimeIndex) -> pd.Timedelta:
 
         Returns:
             pd.Timedelta: The inferred frequency as a pandas Timedelta.
+
+        Raises:
+            ValueError: If the index has fewer than 2 timestamps.
         """
-        if len(index) < 2:
+        minimum_required_length = 2
+        if len(index) < minimum_required_length:
             raise ValueError("Cannot infer frequency from an index with fewer than 2 timestamps.")
 
         # Calculate the differences between consecutive timestamps
         deltas = index.to_series().diff().dropna()
 
         # Find the most common difference
-        inferred_freq = deltas.mode().iloc[0]
-        return inferred_freq
+        return deltas.mode().iloc[0]
 
     def _frequency_matches(self, index: pd.DatetimeIndex) -> bool:
         """Check if the frequency of the input data matches the model frequency.
@@ -214,9 +218,6 @@ def _frequency_matches(self, index: pd.DatetimeIndex) -> bool:
         Returns:
             bool: True if the frequencies match, False otherwise.
         """
-        if not isinstance(index, pd.DatetimeIndex):
-            raise ValueError("The index of the input data must be a pandas DatetimeIndex.")
-
         if index.freq is None:
             input_frequency = self._infer_frequency(index)
         else:
@@ -229,8 +230,7 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         """Predict the median of the lag features for each time step in the context window.
 
         Args:
-            data (ForecastInputDataset): The input data for prediction.
-            This should be a pandas dataframe with lag features.
+            data (ForecastInputDataset): The input data for prediction, this should be a pandas dataframe with lag features.
 
         Returns:
             np.array: The predicted median for each time step in the context window.
@@ -239,6 +239,7 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
 
         Raises:
             ValueError: If the input data is missing any of the required lag features.
+            AttributeError: If the model is not fitted yet.
         """
         if not self.is_fitted:
             msg = "This MedianForecaster instance is not fitted yet"
@@ -330,13 +331,21 @@ def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None
         For example, T-1min, T-2min, T-3min or T-1d, T-2d.
 
         Which lag features are used is determined by the feature engineering step.
+
+        Args:
+            data (ForecastInputDataset): The training data containing lag features.
+            data_val (ForecastInputDataset | None): Optional validation data, not used in this regressor.
+
+        Raises:
+            ValueError: If the input data frequency does not match the model frequency.
+            ValueError: If no lag features are found in the input data.
+
         """
         self.frequency_ = data.sample_interval
         # Check that the frequency of the input data matches frequency of the lags
         if not self._frequency_matches(data.data.index.drop_duplicates()):  # Several training horizons give duplicates
-            raise ValueError(
-                f"The input data frequency ({data.data.index.freq}) does not match the model frequency ({self.frequency_})."
-            )
+            msg = f"The input data frequency ({data.data.index.freq}) does not match the model frequency ({self.frequency_})."
+            raise ValueError(msg)
 
         lag_perfix = f"{data.target_column}_lag_"
         self.feature_names_ = [

From 31ee823d5ee38b142f084fd268c86f2f93aeea19 Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Wed, 7 Jan 2026 16:57:45 +0100
Subject: [PATCH 09/13] fix more linting and formatting

---
 .../models/forecasting/median_forecaster.py   | 38 +++++++++----------
 .../forecasting/test_median_forecaster.py     | 14 +++----
 2 files changed, 25 insertions(+), 27 deletions(-)

diff --git a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
index 3d42983a6..585c5bde2 100644
--- a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
+++ b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
@@ -218,11 +218,7 @@ def _frequency_matches(self, index: pd.DatetimeIndex) -> bool:
         Returns:
             bool: True if the frequencies match, False otherwise.
         """
-        if index.freq is None:
-            input_frequency = self._infer_frequency(index)
-        else:
-            input_frequency = index.freq
-
+        input_frequency = self._infer_frequency(index) if index.freq is None else index.freq
         return input_frequency == self.frequency
 
     @override
@@ -230,7 +226,8 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         """Predict the median of the lag features for each time step in the context window.
 
         Args:
-            data (ForecastInputDataset): The input data for prediction, this should be a pandas dataframe with lag features.
+            data (ForecastInputDataset): The input data for prediction,
+            this should be a pandas dataframe with lag features.
 
         Returns:
             np.array: The predicted median for each time step in the context window.
@@ -263,15 +260,11 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
 
         # Reindex the input data to ensure there are no gaps in the time series.
         # This is important for the autoregressive logic that follows.
-        # Store the original index to return predictions aligned with the input.
-        old_index = input_data.index
         # Create a new date range with the expected frequency.
         new_index = pd.date_range(input_data.index[0], input_data.index[-1], freq=self.frequency_)
         # Reindex the input DataFrame, filling any new timestamps with NaN.
-        input_data = input_data.reindex(new_index, fill_value=np.nan)
-
         # Select only the lag feature columns in the specified order.
-        lag_df = input_data[self.feature_names_]
+        lag_df = input_data.reindex(new_index, fill_value=np.nan)[self.feature_names_]
 
         # Convert the lag DataFrame and its index to NumPy arrays for faster processing.
         lag_array = lag_df.to_numpy()
@@ -311,7 +304,7 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         prediction_df = pd.DataFrame(prediction, index=lag_df.index.to_numpy(), columns=["median"])
 
         # Reindex the prediction DataFrame back to the original input data index.
-        prediction_df = prediction_df.reindex(old_index)
+        prediction_df = prediction_df.reindex(input_data.index)
 
         return ForecastDataset(
             data=prediction_df.dropna().rename(columns={"median": self.config.quantiles[0].format()}),
@@ -343,21 +336,25 @@ def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None
         """
         self.frequency_ = data.sample_interval
         # Check that the frequency of the input data matches frequency of the lags
-        if not self._frequency_matches(data.data.index.drop_duplicates()):  # Several training horizons give duplicates
-            msg = f"The input data frequency ({data.data.index.freq}) does not match the model frequency ({self.frequency_})."
+        # Several training horizons give duplicates
+        if not self._frequency_matches(data.data.index.drop_duplicates()):
+            msg = (
+                f"The input data frequency ({data.data.index.freq}) "
+                f"does not match the model frequency ({self.frequency_})."
+            )
             raise ValueError(msg)
 
-        lag_perfix = f"{data.target_column}_lag_"
+        lag_prefix = f"{data.target_column}_lag_"
         self.feature_names_ = [
-            feature_name for feature_name in data.feature_names if feature_name.startswith(lag_perfix)
+            feature_name for feature_name in data.feature_names if feature_name.startswith(lag_prefix)
         ]
 
         if not self.feature_names_:
-            msg = f"No lag features found in the input data with prefix '{lag_perfix}'."
+            msg = f"No lag features found in the input data with prefix '{lag_prefix}'."
             raise ValueError(msg)
 
         self.lags_to_time_deltas_ = {
-            feature_name: timedelta_from_isoformat(feature_name.replace(lag_perfix, ""))
+            feature_name: timedelta_from_isoformat(feature_name.replace(lag_prefix, ""))
             for feature_name in self.feature_names_
         }
 
@@ -365,7 +362,10 @@ def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None
         lag_deltas = sorted(self.lags_to_time_deltas_.values())
         lag_intervals = [(lag_deltas[i] - lag_deltas[i - 1]).total_seconds() for i in range(1, len(lag_deltas))]
         if not all(interval == lag_intervals[0] for interval in lag_intervals):
-            msg = "Lag features are not evenly spaced. Please ensure lag features are evenly spaced and match the data frequency."
+            msg = (
+                "Lag features are not evenly spaced. "
+                "Please ensure lag features are evenly spaced and match the data frequency."
+            )
             raise ValueError(msg)
 
         self.feature_names_ = sorted(self.feature_names_, key=lambda f: self.lags_to_time_deltas_[f])
diff --git a/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py b/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
index 22eb5bb9b..293461251 100644
--- a/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
+++ b/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
@@ -1,18 +1,16 @@
-# ruff: noqa
 # SPDX-FileCopyrightText: 2017-2025 Contributors to the OpenSTEF project <openstef@lfenergy.org>
 #
 # SPDX-License-Identifier: MPL-2.0
 
 import numpy as np
-from openstef_models.models.forecasting.forecaster import ForecastDataset
-from openstef_models.utils.data_split import DataSplitter
 import pandas as pd
 import pytest
-from openstef_models.models.forecasting.median_forecaster import MedianForecaster, MedianForecasterConfig
+
 from openstef_core.datasets import ForecastInputDataset
-from openstef_core.types import LeadTime, Q
 from openstef_core.testing import create_timeseries_dataset
-from openstef_models.presets import ForecastingWorkflowConfig, create_forecasting_workflow
+from openstef_core.types import LeadTime, Q
+from openstef_models.models.forecasting.forecaster import ForecastDataset
+from openstef_models.models.forecasting.median_forecaster import MedianForecaster, MedianForecasterConfig
 
 
 def test_median_returns_median():
@@ -382,7 +380,7 @@ def test_median_fit_with_no_lag_features_raises():
     model = MedianForecaster(config=config)
 
     # Act & Assert
-    with pytest.raises(ValueError, match="No lag features found in the input data."):
+    with pytest.raises(ValueError, match=r"No lag features found in the input data."):
         model.fit(training_input_data)
 
 
@@ -435,7 +433,7 @@ def test_median_fit_with_inconsistent_frequency_raises():
     model = MedianForecaster(config=config)
 
     # Act & Assert
-    with pytest.raises(ValueError, match="does not match the model frequency."):
+    with pytest.raises(ValueError, match=r"does not match the model frequency."):
         model.fit(training_input_data)
 
 

From 486c43cb8a3ca5c045e59c6bccee8938cfc50ec6 Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Wed, 7 Jan 2026 17:05:08 +0100
Subject: [PATCH 10/13] fix type check

---
 .../models/forecasting/median_forecaster.py   | 20 +++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
index 585c5bde2..d5cd216ab 100644
--- a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
+++ b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
@@ -145,7 +145,7 @@ def feature_importances(self) -> pd.DataFrame:
         )
 
     @property
-    def frequency(self) -> int:
+    def frequency(self) -> timedelta:
         """Retrieve the model input frequency.
 
         Returns:
@@ -250,11 +250,11 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
             msg = f"The input data is missing the following lag features: {missing_features}"
             raise ValueError(msg)
 
-        if not self._frequency_matches(data.input_data().index):
+        if not self._frequency_matches(data.input_data().index):  # type: ignore
             msg = (
                 f"The frequency of the input data does not match the model frequency. "
-                f"Input data frequency: {data.input_data().index.freq}, "
-                f"Model frequency: {pd.Timedelta(minutes=self.frequency_)}"
+                f"Input data frequency: {data.input_data().index.freq}, "  # type: ignore
+                f"Model frequency: {self.frequency_}"
             )
             raise ValueError(msg)
 
@@ -282,10 +282,10 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
             # Get the lag features for the current time step.
             current_lags = lag_array[time_step]
             # Calculate the median of the available lag features, ignoring NaNs.
-            median = np.nanmedian(current_lags)
+            median = np.nanmedian(current_lags)  # type: ignore
             # If the median calculation resulted in NaN (e.g., all lags were NaN), skip the autoregression step.
-            if not np.isnan(median):
-                median = float(median)
+            if not np.isnan(median):  # type: ignore
+                median = float(median)  # type: ignore
             else:
                 continue
 
@@ -307,7 +307,7 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         prediction_df = prediction_df.reindex(input_data.index)
 
         return ForecastDataset(
-            data=prediction_df.dropna().rename(columns={"median": self.config.quantiles[0].format()}),
+            data=prediction_df.dropna().rename(columns={"median": self.config.quantiles[0].format()}),  # type: ignore
             sample_interval=data.sample_interval,
             forecast_start=data.forecast_start,
         )
@@ -337,9 +337,9 @@ def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None
         self.frequency_ = data.sample_interval
         # Check that the frequency of the input data matches frequency of the lags
         # Several training horizons give duplicates
-        if not self._frequency_matches(data.data.index.drop_duplicates()):
+        if not self._frequency_matches(data.data.index.drop_duplicates()):  # type: ignore
             msg = (
-                f"The input data frequency ({data.data.index.freq}) "
+                f"The input data frequency ({data.data.index.freq}) "  # type: ignore
                 f"does not match the model frequency ({self.frequency_})."
             )
             raise ValueError(msg)

From 031c780bd71129ae7eab402fa76be866909f1ecd Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Fri, 9 Jan 2026 13:08:04 +0100
Subject: [PATCH 11/13] Implemented comments

---
 .../models/forecasting/median_forecaster.py         |  9 ++-------
 .../openstef_models/presets/forecasting_workflow.py |  2 +-
 .../models/forecasting/test_median_forecaster.py    | 13 -------------
 3 files changed, 3 insertions(+), 21 deletions(-)

diff --git a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
index d5cd216ab..29a1ff7c0 100644
--- a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
+++ b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
@@ -5,8 +5,6 @@
 """Median regressor based forecasting models for energy forecasting.
 
 Provides median regression models for multi-quantile energy forecasting.
-Optimized for time series data with specialized loss functions and
-comprehensive hyperparameter control for production forecasting workflows.
 
 Note that this is a autoregressive model, meaning that it uses the previous
     predictions to predict the next value.
@@ -20,7 +18,7 @@
     over the last few timesteps adds some hysterisis to avoid triggering on noise.
 
     Tips for using this regressor:
-    - Set the lags to be evenly spaced and at a frequency mathching the
+    - Set the lags to be evenly spaced and at a frequency matching the
       frequency of the input data. For example, if the input data is at 15
       minute intervals, set the lags to be at 15 minute intervals as well.
     - Use a small training dataset, since there are no actual parameters to train.
@@ -28,9 +26,6 @@
     a problem if we get very small chunks of data in training or validation sets.
     - Use only one training horizon, since the regressor will use the same lags for all
       training horizons.
-    - Allow for missing data by setting completeness_threshold to 0. If the prediction horizon
-      is larger than the context window there will be a lot of nans in the input data, but
-      the autoregression solves that.
 """
 
 from datetime import timedelta
@@ -314,7 +309,7 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
 
     @override
     def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None = None) -> None:
-        """Take car of fitting the median forecaster.
+        """Take care of fitting the median forecaster.
 
         This regressor does not need any fitting,
         but it does need to know the feature names of the lag features and the order of these.
diff --git a/packages/openstef-models/src/openstef_models/presets/forecasting_workflow.py b/packages/openstef-models/src/openstef_models/presets/forecasting_workflow.py
index 2de522083..bed8c2c2e 100644
--- a/packages/openstef-models/src/openstef_models/presets/forecasting_workflow.py
+++ b/packages/openstef-models/src/openstef_models/presets/forecasting_workflow.py
@@ -102,7 +102,7 @@ class ForecastingWorkflowConfig(BaseConfig):  # PredictionJob
     # Model configuration
     model: Literal["xgboost", "gblinear", "flatliner", "median"] = Field(
         description="Type of forecasting model to use."
-    )  # TODO(#652): Implement median forecaster
+    )
     quantiles: list[Quantile] = Field(
         default=[Q(0.5)], description="List of quantiles to predict for probabilistic forecasting."
     )
diff --git a/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py b/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
index 293461251..35408de93 100644
--- a/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
+++ b/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
@@ -22,7 +22,6 @@ def test_median_returns_median():
         load_lag_PT1H=[1.0, np.nan, np.nan],
         load_lag_PT2H=[4.0, 1.0, np.nan],
         load_lag_PT3H=[7.0, 4.0, 1.0],
-        available_at=index,
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(
@@ -63,7 +62,6 @@ def test_median_handles_some_missing_data():
         load_lag_PT1H=[1.0, np.nan, np.nan],
         load_lag_PT2H=[np.nan, 1.0, np.nan],
         load_lag_PT3H=[3.0, np.nan, 1.0],
-        available_at=index,
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(
@@ -104,7 +102,6 @@ def test_median_handles_missing_data_for_some_horizons():
         load_lag_PT1H=[5.0, np.nan, np.nan],
         load_lag_PT2H=[np.nan, 5.0, np.nan],
         load_lag_PT3H=[np.nan, np.nan, 5.0],
-        available_at=index,
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(
@@ -145,7 +142,6 @@ def test_median_handles_all_missing_data():
         load_lag_PT1H=[np.nan, np.nan, np.nan],
         load_lag_PT2H=[np.nan, np.nan, np.nan],
         load_lag_PT3H=[np.nan, np.nan, np.nan],
-        available_at=index,
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(
@@ -186,7 +182,6 @@ def test_median_uses_lag_features_if_available():
         load_lag_PT1H=[1.0, 2.0, 3.0],
         load_lag_PT2H=[4.0, 5.0, 6.0],
         load_lag_PT3H=[7.0, 8.0, 9.0],
-        available_at=index,
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(
@@ -235,7 +230,6 @@ def test_median_handles_small_gap():
         load_lag_PT3H=[3.0, 2.0, 1.0, np.nan, np.nan],
         load_lag_PT4H=[4.0, 3.0, 2.0, 1.0, np.nan],
         load_lag_PT5H=[5.0, 4.0, 3.0, 2.0, 1.0],
-        available_at=index,
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(
@@ -278,7 +272,6 @@ def test_median_handles_large_gap():
         load=[1.0, 2.0, 3.0],
         load_lag_PT1H=[4.0, 5.0, 6.0],
         load_lag_PT2H=[7.0, 8.0, 9.0],
-        available_at=index_1,
     )
 
     index_2 = pd.date_range("2023-01-02T01:00", periods=3, freq="h")
@@ -287,7 +280,6 @@ def test_median_handles_large_gap():
         load=[10.0, 11.0, 12.0],
         load_lag_PT1H=[13.0, 14.0, 15.0],
         load_lag_PT2H=[16.0, 17.0, 18.0],
-        available_at=index_2,
     )
 
     training_data = training_data_1
@@ -332,7 +324,6 @@ def test_median_fit_with_missing_features_raises():
         load_lag_PT1H=[1.0, 2.0, 3.0],
         load_lag_PT2H=[4.0, 5.0, 6.0],
         load_lag_PT3H=[7.0, 8.0, 9.0],
-        available_at=index,
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(
@@ -367,7 +358,6 @@ def test_median_fit_with_no_lag_features_raises():
         index=index,
         load=[1.0, 2.0, 3.0],
         unrelated_feature=[1.0, 2.0, 3.0],
-        available_at=index,
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(
@@ -394,7 +384,6 @@ def test_median_fit_with_inconsistent_lag_features_raises():
         load_lag_PT5H=[4.0, 5.0, 6.0],
         load_lag_PT60H=[7.0, 8.0, 9.0],
         load_lag_P4D=[10.0, 11.0, 12.0],
-        available_at=index,
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(
@@ -420,7 +409,6 @@ def test_median_fit_with_inconsistent_frequency_raises():
         load_lag_PT1H=[1.0, 2.0, 3.0],
         load_lag_PT2H=[4.0, 5.0, 6.0],
         load_lag_PT3H=[7.0, 8.0, 9.0],
-        available_at=index,
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(
@@ -446,7 +434,6 @@ def test_predicting_without_fitting_raises():
         load_lag_PT1M=[1.0, 2.0, 3.0],
         load_lag_PT2M=[4.0, 5.0, 6.0],
         load_lag_PT3M=[7.0, 8.0, 9.0],
-        available_at=index,
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(

From 7266d8ea12096949b9f072a4a8c3559d9c074308 Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Mon, 12 Jan 2026 10:30:16 +0100
Subject: [PATCH 12/13] moved sample interval check to TimeseriesDataset and
 updated tests accordingly

---
 .../backtesting/backtest_pipeline.py          |  1 -
 .../backtesting/test_backtest_pipeline.py     | 11 ++--
 .../unit/backtesting/test_batch_prediction.py |  2 -
 .../unit/benchmarking/test_target_provider.py |  6 +-
 .../datasets/timeseries_dataset.py            | 57 ++++++++++++++++-
 .../models/forecasting/median_forecaster.py   | 64 +++----------------
 .../forecasting/test_median_forecaster.py     | 11 +++-
 .../test_confidence_interval_applicator.py    |  4 +-
 .../test_datetime_features_adder.py           |  6 +-
 .../test_rolling_aggregates_adder.py          |  2 +-
 .../test_daylight_feature_adder.py            |  4 +-
 11 files changed, 93 insertions(+), 75 deletions(-)

diff --git a/packages/openstef-beam/src/openstef_beam/backtesting/backtest_pipeline.py b/packages/openstef-beam/src/openstef_beam/backtesting/backtest_pipeline.py
index 65a8393af..597f6ac5f 100644
--- a/packages/openstef-beam/src/openstef_beam/backtesting/backtest_pipeline.py
+++ b/packages/openstef-beam/src/openstef_beam/backtesting/backtest_pipeline.py
@@ -180,7 +180,6 @@ def run(
 
         return TimeSeriesDataset(
             data=pd.concat([pred.data for pred in prediction_list], axis=0),
-            sample_interval=self.config.prediction_sample_interval,
         )
 
     def _process_train_event(self, event: BacktestEvent, dataset: VersionedTimeSeriesDataset) -> None:
diff --git a/packages/openstef-beam/tests/unit/backtesting/test_backtest_pipeline.py b/packages/openstef-beam/tests/unit/backtesting/test_backtest_pipeline.py
index 9e68fd5a1..cfb1b6778 100644
--- a/packages/openstef-beam/tests/unit/backtesting/test_backtest_pipeline.py
+++ b/packages/openstef-beam/tests/unit/backtesting/test_backtest_pipeline.py
@@ -139,8 +139,7 @@ def test_run_training_scenarios(
 
     # Assert
     assert isinstance(result, TimeSeriesDataset)
-    assert result.sample_interval == forecaster_config.predict_sample_interval
-
+    assert result.sample_interval == timedelta(hours=6)
     # Validate call counts
     if expected_train_calls == ">0":
         assert mock_forecaster.train_call_count > 0
@@ -245,7 +244,7 @@ def create_prediction(data: RestrictedHorizonVersionedTimeSeries) -> TimeSeriesD
     )
 
     # Assert - Basic structure
-    assert result.sample_interval == mock_forecaster.config.predict_sample_interval
+    assert result.sample_interval == timedelta(hours=6)
     assert mock_forecaster.predict_call_count >= 2
 
     # Assert - Output validation
@@ -352,18 +351,20 @@ def test_run_edge_cases(
         timestamps = pd.date_range("2025-01-01T12:00:00", "2025-01-01T15:00:00", freq="1h")
         start_time = "2025-01-01T12:00:00"
         end_time = "2025-01-01T15:00:00"
+        sample_interval = timedelta(hours=1)
     else:  # sparse
         timestamps = pd.DatetimeIndex(["2025-01-01T12:00:00", "2025-01-01T18:00:00"])
         start_time = "2025-01-01T18:00:00"
         end_time = "2025-01-01T20:00:00"
+        sample_interval = timedelta(hours=6)
 
     ground_truth = VersionedTimeSeriesDataset.from_dataframe(
         data=pd.DataFrame({"available_at": timestamps, "target": range(len(timestamps))}, index=timestamps),
-        sample_interval=timedelta(hours=1),
+        sample_interval=sample_interval,
     )
     predictors = VersionedTimeSeriesDataset.from_dataframe(
         data=pd.DataFrame({"available_at": timestamps, "feature1": range(len(timestamps))}, index=timestamps),
-        sample_interval=timedelta(hours=1),
+        sample_interval=sample_interval,
     )
 
     # Act
diff --git a/packages/openstef-beam/tests/unit/backtesting/test_batch_prediction.py b/packages/openstef-beam/tests/unit/backtesting/test_batch_prediction.py
index 09e6ac4f0..1863d8ec6 100644
--- a/packages/openstef-beam/tests/unit/backtesting/test_batch_prediction.py
+++ b/packages/openstef-beam/tests/unit/backtesting/test_batch_prediction.py
@@ -55,7 +55,6 @@ def predict(self, data: RestrictedHorizonVersionedTimeSeries) -> TimeSeriesDatas
         timestamps = pd.date_range(start=data.horizon, periods=2, freq="1h")
         return TimeSeriesDataset(
             data=pd.DataFrame({"quantile_P50": [0.5, 0.5]}, index=timestamps),
-            sample_interval=self.config.predict_sample_interval,
         )
 
     def predict_batch(self, batch: list[RestrictedHorizonVersionedTimeSeries]) -> list[TimeSeriesDataset]:
@@ -66,7 +65,6 @@ def predict_batch(self, batch: list[RestrictedHorizonVersionedTimeSeries]) -> li
             results.append(
                 TimeSeriesDataset(
                     data=pd.DataFrame({"quantile_P50": [0.5, 0.5]}, index=timestamps),
-                    sample_interval=self.config.predict_sample_interval,
                 )
             )
         return results
diff --git a/packages/openstef-beam/tests/unit/benchmarking/test_target_provider.py b/packages/openstef-beam/tests/unit/benchmarking/test_target_provider.py
index 1692f92b4..61b7f1b09 100644
--- a/packages/openstef-beam/tests/unit/benchmarking/test_target_provider.py
+++ b/packages/openstef-beam/tests/unit/benchmarking/test_target_provider.py
@@ -82,13 +82,13 @@ def test_get_predictors_for_target(tmp_path: Path, test_target: BenchmarkTarget)
     interval = timedelta(hours=1)
 
     weather = VersionedTimeSeriesDataset.from_dataframe(
-        pd.DataFrame({"temp": range(3), "available_at": index}, index=index), interval
+        pd.DataFrame({"temp": range(3), "available_at": index}, index=index), sample_interval=interval
     )
     profiles = VersionedTimeSeriesDataset.from_dataframe(
-        pd.DataFrame({"prof": range(3), "available_at": index}, index=index), interval
+        pd.DataFrame({"prof": range(3), "available_at": index}, index=index), sample_interval=interval
     )
     prices = VersionedTimeSeriesDataset.from_dataframe(
-        pd.DataFrame({"price": range(3), "available_at": index}, index=index), interval
+        pd.DataFrame({"price": range(3), "available_at": index}, index=index), sample_interval=interval
     )
 
     class TestProvider(SimpleTargetProvider[BenchmarkTarget, None]):
diff --git a/packages/openstef-core/src/openstef_core/datasets/timeseries_dataset.py b/packages/openstef-core/src/openstef_core/datasets/timeseries_dataset.py
index d4dfb63d2..c316ba7c5 100644
--- a/packages/openstef-core/src/openstef_core/datasets/timeseries_dataset.py
+++ b/packages/openstef-core/src/openstef_core/datasets/timeseries_dataset.py
@@ -99,7 +99,7 @@ class TimeSeriesDataset(TimeSeriesMixin, DatasetMixin):  # noqa: PLR0904 - impor
     def __init__(
         self,
         data: pd.DataFrame,
-        sample_interval: timedelta = timedelta(minutes=15),
+        sample_interval: timedelta | None = None,
         *,
         horizon_column: str = "horizon",
         available_at_column: str = "available_at",
@@ -122,10 +122,27 @@ def __init__(
         Raises:
             TypeError: If data index is not a pandas DatetimeIndex or if versioning
                 columns have incorrect types.
+            ValueError: If data frequency does not match sample_interval.
         """
         if not isinstance(data.index, pd.DatetimeIndex):
             raise TypeError("Data index must be a pandas DatetimeIndex.")
 
+        if sample_interval is None:
+            inferred_freq = pd.Timedelta(
+                self._infer_frequency(data.index) if data.index.freq is None else data.index.freq  # type: ignore
+            )
+            sample_interval = inferred_freq.to_pytimedelta()
+
+        # Check input data frequency matches sample_interval, only if there are enough data points to infer frequency
+        minimum_required_length = 2
+        if len(data) >= minimum_required_length:
+            input_sample_interval = self._infer_frequency(data.index) if data.index.freq is None else data.index.freq
+            if input_sample_interval != sample_interval:
+                msg = (
+                    f"Data frequency ({input_sample_interval}) does not match the sample_interval ({sample_interval})."
+                )
+                raise ValueError(msg)
+
         self.data = data
         self.horizon_column = horizon_column
         self.available_at_column = available_at_column
@@ -443,6 +460,44 @@ def copy_with(self, data: pd.DataFrame, *, is_sorted: bool = False) -> "TimeSeri
             is_sorted=is_sorted,
         )
 
+    @staticmethod
+    def _infer_frequency(index: pd.DatetimeIndex) -> pd.Timedelta:
+        """Infer the frequency of a pandas DatetimeIndex if the freq attribute is not set.
+
+        This method calculates the most common time difference between consecutive timestamps,
+        which is more permissive of missing chunks of data than the pandas infer_freq method.
+
+        Args:
+            index (pd.DatetimeIndex): The datetime index to infer the frequency from.
+
+        Returns:
+            pd.Timedelta: The inferred frequency as a pandas Timedelta.
+
+        Raises:
+            ValueError: If the index has fewer than 2 timestamps.
+        """
+        minimum_required_length = 2
+        if len(index) < minimum_required_length:
+            raise ValueError("Cannot infer frequency from an index with fewer than 2 timestamps.")
+
+        # Calculate the differences between consecutive timestamps
+        deltas = index.to_series().drop_duplicates().sort_values().diff().dropna()
+
+        # Find the most common difference
+        return deltas.mode().iloc[0]
+
+    def _frequency_matches(self, index: pd.DatetimeIndex) -> bool:
+        """Check if the frequency of the data matches the model frequency.
+
+        Args:
+            index (pd.DatetimeIndex): The data to check.
+
+        Returns:
+            bool: True if the frequencies match, False otherwise.
+        """
+        input_sample_interval = self._infer_frequency(index) if index.freq is None else index.freq
+        return input_sample_interval == self.sample_interval
+
 
 def validate_horizons_present(dataset: TimeSeriesDataset, horizons: list[LeadTime]) -> None:
     """Validate that the specified forecast horizons are present in the dataset.
diff --git a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
index 29a1ff7c0..b6449c78f 100644
--- a/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
+++ b/packages/openstef-models/src/openstef_models/models/forecasting/median_forecaster.py
@@ -178,44 +178,6 @@ def _fill_diagonal_with_median(lag_array: np.ndarray, start: int, end: int, medi
             np.fill_diagonal(view, updated_diagonal)
         return None
 
-    @staticmethod
-    def _infer_frequency(index: pd.DatetimeIndex) -> pd.Timedelta:
-        """Infer the frequency of a pandas DatetimeIndex if the freq attribute is not set.
-
-        This method calculates the most common time difference between consecutive timestamps,
-        which is more permissive of missing chunks of data than the pandas infer_freq method.
-
-        Args:
-            index (pd.DatetimeIndex): The datetime index to infer the frequency from.
-
-        Returns:
-            pd.Timedelta: The inferred frequency as a pandas Timedelta.
-
-        Raises:
-            ValueError: If the index has fewer than 2 timestamps.
-        """
-        minimum_required_length = 2
-        if len(index) < minimum_required_length:
-            raise ValueError("Cannot infer frequency from an index with fewer than 2 timestamps.")
-
-        # Calculate the differences between consecutive timestamps
-        deltas = index.to_series().diff().dropna()
-
-        # Find the most common difference
-        return deltas.mode().iloc[0]
-
-    def _frequency_matches(self, index: pd.DatetimeIndex) -> bool:
-        """Check if the frequency of the input data matches the model frequency.
-
-        Args:
-            index (pd.DatetimeIndex): The input data to check.
-
-        Returns:
-            bool: True if the frequencies match, False otherwise.
-        """
-        input_frequency = self._infer_frequency(index) if index.freq is None else index.freq
-        return input_frequency == self.frequency
-
     @override
     def predict(self, data: ForecastInputDataset) -> ForecastDataset:
         """Predict the median of the lag features for each time step in the context window.
@@ -245,14 +207,6 @@ def predict(self, data: ForecastInputDataset) -> ForecastDataset:
             msg = f"The input data is missing the following lag features: {missing_features}"
             raise ValueError(msg)
 
-        if not self._frequency_matches(data.input_data().index):  # type: ignore
-            msg = (
-                f"The frequency of the input data does not match the model frequency. "
-                f"Input data frequency: {data.input_data().index.freq}, "  # type: ignore
-                f"Model frequency: {self.frequency_}"
-            )
-            raise ValueError(msg)
-
         # Reindex the input data to ensure there are no gaps in the time series.
         # This is important for the autoregressive logic that follows.
         # Create a new date range with the expected frequency.
@@ -330,14 +284,6 @@ def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None
 
         """
         self.frequency_ = data.sample_interval
-        # Check that the frequency of the input data matches frequency of the lags
-        # Several training horizons give duplicates
-        if not self._frequency_matches(data.data.index.drop_duplicates()):  # type: ignore
-            msg = (
-                f"The input data frequency ({data.data.index.freq}) "  # type: ignore
-                f"does not match the model frequency ({self.frequency_})."
-            )
-            raise ValueError(msg)
 
         lag_prefix = f"{data.target_column}_lag_"
         self.feature_names_ = [
@@ -363,6 +309,16 @@ def fit(self, data: ForecastInputDataset, data_val: ForecastInputDataset | None
             )
             raise ValueError(msg)
 
+        # Check that lag frequency matches data frequency
+        expected_lag_interval = lag_intervals[0]
+        if expected_lag_interval != self.frequency_.total_seconds():
+            msg = (
+                f"Lag feature interval ({timedelta(seconds=expected_lag_interval)}) does not match "
+                f"data frequency ({self.frequency_}). "
+                "Please ensure lag features match the data frequency."
+            )
+            raise ValueError(msg)
+
         self.feature_names_ = sorted(self.feature_names_, key=lambda f: self.lags_to_time_deltas_[f])
 
         self.feature_importances_ = np.ones(len(self.feature_names_)) / (len(self.feature_names_) or 1.0)
diff --git a/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py b/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
index 35408de93..07ca2dfe5 100644
--- a/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
+++ b/packages/openstef-models/tests/unit/models/forecasting/test_median_forecaster.py
@@ -1,6 +1,8 @@
 # SPDX-FileCopyrightText: 2017-2025 Contributors to the OpenSTEF project <openstef@lfenergy.org>
 #
 # SPDX-License-Identifier: MPL-2.0
+import re
+from datetime import timedelta
 
 import numpy as np
 import pandas as pd
@@ -409,6 +411,7 @@ def test_median_fit_with_inconsistent_frequency_raises():
         load_lag_PT1H=[1.0, 2.0, 3.0],
         load_lag_PT2H=[4.0, 5.0, 6.0],
         load_lag_PT3H=[7.0, 8.0, 9.0],
+        sample_interval=timedelta(minutes=1),
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(
@@ -421,7 +424,12 @@ def test_median_fit_with_inconsistent_frequency_raises():
     model = MedianForecaster(config=config)
 
     # Act & Assert
-    with pytest.raises(ValueError, match=r"does not match the model frequency."):
+    with pytest.raises(
+        ValueError,
+        match=re.escape(
+            r"Lag feature interval (1:00:00) does not match data frequency (0:01:00). Please ensure lag features match the data frequency."
+        ),
+    ):
         model.fit(training_input_data)
 
 
@@ -434,6 +442,7 @@ def test_predicting_without_fitting_raises():
         load_lag_PT1M=[1.0, 2.0, 3.0],
         load_lag_PT2M=[4.0, 5.0, 6.0],
         load_lag_PT3M=[7.0, 8.0, 9.0],
+        sample_interval=timedelta(minutes=1),
     )
 
     training_input_data = ForecastInputDataset.from_timeseries(
diff --git a/packages/openstef-models/tests/unit/transforms/postprocessing/test_confidence_interval_applicator.py b/packages/openstef-models/tests/unit/transforms/postprocessing/test_confidence_interval_applicator.py
index c07399345..79a16dd17 100644
--- a/packages/openstef-models/tests/unit/transforms/postprocessing/test_confidence_interval_applicator.py
+++ b/packages/openstef-models/tests/unit/transforms/postprocessing/test_confidence_interval_applicator.py
@@ -32,7 +32,7 @@ def validation_predictions() -> ForecastDataset:
             },
             index=index,
         ),
-        sample_interval=timedelta(minutes=15),
+        sample_interval=timedelta(hours=1),
         forecast_start=datetime.fromisoformat("2018-01-01T00:00:00"),
     )
 
@@ -51,7 +51,7 @@ def predictions() -> ForecastDataset:
             },
             index=forecast_index,
         ),
-        sample_interval=timedelta(minutes=15),
+        sample_interval=timedelta(hours=1),
         forecast_start=forecast_start,
     )
 
diff --git a/packages/openstef-models/tests/unit/transforms/time_domain/test_datetime_features_adder.py b/packages/openstef-models/tests/unit/transforms/time_domain/test_datetime_features_adder.py
index 79b52f2e4..d4bad586a 100644
--- a/packages/openstef-models/tests/unit/transforms/time_domain/test_datetime_features_adder.py
+++ b/packages/openstef-models/tests/unit/transforms/time_domain/test_datetime_features_adder.py
@@ -48,7 +48,7 @@ def test_datetime_features_weekday_weekend():
         {"value": [1.0, 2.0]},
         index=pd.DatetimeIndex(["2025-01-06", "2025-01-11"]),  # Monday, Saturday
     )
-    input_data = TimeSeriesDataset(data, timedelta(days=1))
+    input_data = TimeSeriesDataset(data, timedelta(days=5))
 
     transform = DatetimeFeaturesAdder()
     result = transform.transform(input_data)
@@ -82,7 +82,7 @@ def test_datetime_features_month_quarter():
         {"value": [1.0, 2.0, 3.0]},
         index=pd.DatetimeIndex(["2025-01-15", "2025-04-15", "2025-10-15"]),  # Jan Q1, Apr Q2, Oct Q4
     )
-    input_data = TimeSeriesDataset(data, timedelta(days=1))
+    input_data = TimeSeriesDataset(data, timedelta(days=90))
 
     transform = DatetimeFeaturesAdder()
     result = transform.transform(input_data)
@@ -102,7 +102,7 @@ def test_datetime_features_onehot_encoding():
         {"value": [1.0, 2.0]},
         index=pd.DatetimeIndex(["2025-01-15", "2025-04-15"]),  # Jan, Apr
     )
-    input_data = TimeSeriesDataset(data, timedelta(days=1))
+    input_data = TimeSeriesDataset(data, timedelta(days=90))
 
     transform = DatetimeFeaturesAdder(onehot_encode=True)
     result = transform.transform(input_data)
diff --git a/packages/openstef-models/tests/unit/transforms/time_domain/test_rolling_aggregates_adder.py b/packages/openstef-models/tests/unit/transforms/time_domain/test_rolling_aggregates_adder.py
index b74a42317..fe7e0bd2c 100644
--- a/packages/openstef-models/tests/unit/transforms/time_domain/test_rolling_aggregates_adder.py
+++ b/packages/openstef-models/tests/unit/transforms/time_domain/test_rolling_aggregates_adder.py
@@ -92,7 +92,7 @@ def test_rolling_aggregate_features_missing_column_raises_error():
         {"not_load": [1.0, 2.0, 3.0]},
         index=pd.date_range("2023-01-01 00:00:00", periods=3, freq="1h"),
     )
-    dataset = TimeSeriesDataset(data, sample_interval=timedelta(minutes=15))
+    dataset = TimeSeriesDataset(data, sample_interval=timedelta(hours=1))
     transform = RollingAggregatesAdder(
         feature="load",
         horizons=[LeadTime.from_string("PT36H")],
diff --git a/packages/openstef-models/tests/unit/transforms/weather_domain/test_daylight_feature_adder.py b/packages/openstef-models/tests/unit/transforms/weather_domain/test_daylight_feature_adder.py
index 84b31ca38..a95de7aed 100644
--- a/packages/openstef-models/tests/unit/transforms/weather_domain/test_daylight_feature_adder.py
+++ b/packages/openstef-models/tests/unit/transforms/weather_domain/test_daylight_feature_adder.py
@@ -79,10 +79,10 @@ def test_daylight_features_nighttime_values():
 def test_daylight_features_different_coordinates():
     """Test daylight calculation with different geographical coordinates."""
     # Same time, different locations should give different results
-    time_index = pd.DatetimeIndex(["2025-06-21 12:00:00"], tz="UTC")
+    time_index = pd.DatetimeIndex(["2025-06-21 12:00:00", "2025-06-21 13:00:00"], tz="UTC")
 
     # Amsterdam vs Cape Town (southern hemisphere, winter)
-    data = pd.DataFrame({"value": [1.0]}, index=time_index)
+    data = pd.DataFrame({"value": [1.0, 2.0]}, index=time_index)
     input_data = TimeSeriesDataset(data, timedelta(hours=1))
 
     # Amsterdam (northern hemisphere, summer)

From ef9db79629a7e7c43aa3574887e95da881bf63ec Mon Sep 17 00:00:00 2001
From: Jan Maarten van Doorn <janmaarten.vandoorn@sigholm.se>
Date: Mon, 12 Jan 2026 10:37:13 +0100
Subject: [PATCH 13/13] fix doc test

---
 .../src/openstef_models/transforms/general/nan_dropper.py       | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/packages/openstef-models/src/openstef_models/transforms/general/nan_dropper.py b/packages/openstef-models/src/openstef_models/transforms/general/nan_dropper.py
index e9144bba9..cf85856f5 100644
--- a/packages/openstef-models/src/openstef_models/transforms/general/nan_dropper.py
+++ b/packages/openstef-models/src/openstef_models/transforms/general/nan_dropper.py
@@ -38,7 +38,7 @@ class NaNDropper(BaseConfig, TimeSeriesTransform):
         ...     'temperature': [20.0, 22.0, np.nan, 23.0],
         ...     'humidity': [60.0, 65.0, 70.0, 75.0]
         ... }, index=pd.date_range('2025-01-01', periods=4, freq='1h'))
-        >>> dataset = TimeSeriesDataset(data, timedelta(hours=1))
+        >>> dataset = TimeSeriesDataset(data)
         >>>
         >>> # Drop rows with NaN in load or temperature
         >>> dropper = NaNDropper(selection=FeatureSelection(include=['load', 'temperature']))