Score forecasts against reality; expose forecast trajectories

There was no evidence of how accurate the thermal model is: the only
number was an in-sample 5-minute error.

At every poll, record each room's reading, the outdoor temperature, the
radiator-lag estimate and the model parameters. 15, 30, 60 and 120
minutes later, replay the room from that poll with the parameters it
had then (out of sample), the outdoor temperature held constant and the
radiator lag that actually followed from the relay's real history, and
compare with the reading now. A no-change forecast is scored alongside.
Only spans with regular polls and a well-known lag are scored. 7 days
of per-room statistics are saved with the model (optional "evaluation"
payload key; about 18 KB at most).

Four new sensors, Heating MPC Forecast Error 15/30/60/120 min, show the
mean absolute error with bias, persistence error, skill and per-room
detail. On the synthetic lagged plant the 60-minute error is 0.03 °C
against 0.39 °C for persistence.

Room forecast entities gain forecast_15min and plan_15min attributes
and are labelled counterfactual. Volatile attributes are unrecorded.
The model is now saved after scoring, so saved statistics are current.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

README.md

@@ -50,8 +50,13 @@ Entities (entity IDs may vary):
 - `sensor.heating_mpc_status` — learning/ready/paused status and diagnostics, including outside weather, active source, all valid alternative weather values, model coefficients, proposed 15-minute plan.
 - `sensor.heating_mpc_weighted_temperature` — own room-weighted average, informative only.
 - `sensor.mpc_outdoor_temperature` — selected external temp and alternatives.
-- Per-room predicted temperatures after 60 minutes.
+- Per-room predicted temperatures after 60 minutes. Attributes `forecast_15min` (the
+  temperature every 15 minutes over the horizon) and `plan_15min` (the proposed
+  ON/OFF plan). These assume the MPC's plan is executed; while Generic Thermostat
+  controls the relay they are **counterfactual**.
 - `binary_sensor.heating_mpc_proposed_heat` — **hypothetical** heating demand (unavailable until ready).
+- **Heating MPC Forecast Error 15/30/60/120 min** — how far the thermal model's
+  forecasts were from what actually happened (see *Forecast accuracy*).
 
 ## Known limitations / safety
 
@@ -126,6 +131,40 @@ After a day of running, set each limit comfortably above (about 1.5–2×) the l
 `max_age_seconds_24h` among the sensors it covers. The relay and thermostat have no
 age limit.
 
+## Forecast accuracy
+
+The four **Forecast Error** sensors measure the thermal model against reality,
+continuously and out of sample. At every poll the integration records each room's
+reading, the outdoor temperature, the radiator-lag estimate and the model's
+parameters. 15, 30, 60 and 120 minutes later it replays that room from the earlier
+reading using the parameters it had **then**, the outdoor temperature held at its
+starting value, and the radiator heat that actually followed from the relay's real
+ON/OFF history. The result is compared with the current reading.
+
+This tests the thermal model, not the MPC plan, because the plan is never
+executed. Each sensor's state is the mean absolute error (°C) over the last 7 days,
+across all rooms. Attributes, all in °C except counts:
+
+- `samples` — number of scored forecasts. Consecutive forecasts overlap, so they
+  are not independent; treat small counts as indicative only.
+- `bias` — mean error: negative means the model predicts too cold.
+- `persistence_mae` — error of simply assuming no change; the baseline to beat.
+- `skill` — `1 - mae / persistence_mae`: 1 is perfect, 0 is no better than
+  assuming no change, negative is worse. Empty while the temperature did not change.
+- `rooms` — the same numbers per room (not written to the Recorder).
+
+Forecasts are scored only where the inputs were known: polls 4.5–5.5 minutes
+apart, a usable reading at the start and the end, and a well-known radiator-lag
+estimate throughout. A missing relay, outdoor temperature or all rooms restarts the
+replay window, so scores resume 15–120 minutes later. Statistics are kept with the
+model and survive restarts; deleting the model resets them.
+
+What it does not cover: outdoor temperature changes during the horizon (held
+constant, as in the MPC), and periods when the model was still learning (included,
+so early numbers are pessimistic). The roadmap's working target is about
+0.3–0.4 °C at 60 minutes with clearly positive skill, before any control is
+considered.
+
 ## Release notes
 
 See [CHANGELOG.md](CHANGELOG.md) for changes and upgrade notes per version, and

custom_components/heating_mpc/coordinator.py

@@ -20,6 +20,7 @@ from .const import (
     CONF_BAND, DEFAULTS,
 )
 from .data_quality import ReportAgeHistory, freshness
+from .evaluation import ForecastEvaluator
 from .physics import (
     LAG_LEARN_MAX_UNCERTAINTY, LAG_TAU_SECONDS, Learner, advance_lag,
     lag_from_time_in_state, optimize, relay_segments,
@@ -45,6 +46,7 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
         self.learners: dict[str, Learner] = {room: Learner() for room in self.rooms}
         self._retired: dict[str, dict[str, Any]] = {}
         self._report_ages = ReportAgeHistory()
+        self._evaluator = ForecastEvaluator(LAG_LEARN_MAX_UNCERTAINTY)
         # How each poll was handled since start, to explain slow learning.
         self._sample_results: dict[str, int] = {}
         self._store = Store(hass, 1, f"{DOMAIN}.{entry.entry_id}")
@@ -68,6 +70,7 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
         if stored.get("heater_entity", self.config[CONF_HEATER]) != self.config[CONF_HEATER]:
             _LOGGER.warning("Heating MPC heater changed; starting a fresh thermal model")
             return
+        self._evaluator.load(stored.get("evaluation"))
         # Models carry room identity, so changing the room list is safe. Models of
         # rooms no longer configured are retained, most recently removed last.
         for key in ("retired_learners", "learners"):  # current learners take precedence
@@ -255,10 +258,13 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
         vals: dict[str, float] = {room: float(t) for room, t in temps.items() if t is not None}
         if not heater_valid or climate is None or outside is None or not vals:
             self._previous = None
+            self._evaluator.reset_window()
+            status["forecast_accuracy"] = self._evaluator.summary(self.rooms)
             return status
 
         # A room learns only across a poll pair where it was usable at both ends.
         pairs = [room for room in vals if self._previous and room in self._previous["temps"]]
+        learned = False
         if not pairs:
             status["sample_result"] = "first_sample"
         else:
@@ -282,16 +288,7 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
                 ]
                 status["sample_result"] = "learned" if all(accepted) else "partly_rejected"
                 self._update_counts()
-                # Persist at each sample so restart does not erase hours of training.
+                learned = True
-                await self._store.async_save({
-                    "schema_version": 2,
-                    "learners": {room: learner.as_dict() for room, learner in self.learners.items()},
-                    # Optional; older 0.2.x readers ignore it and keep working.
-                    "retired_learners": self._retired,
-                    "on_samples": self._on_samples, "off_samples": self._off_samples,
-                    "heater_entity": self.config[CONF_HEATER],
-                    "outside_source": active_sources["temperature"],
-                })
 
         result = status["sample_result"]
         self._sample_results[result] = self._sample_results.get(result, 0) + 1
@@ -300,6 +297,22 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
             "time": now, "temps": vals, "outside": outside, "heater_on": heater_on,
             "outside_source": active_sources["temperature"],
         }
+        # Parameters after this poll's learning: forecasts from here use only data so far.
+        self._evaluator.observe(now, vals, outside, self._lag["lag"], self._lag["uncertainty"],
+                                {room: tuple(self.learners[room].theta) for room in vals})
+        status["forecast_accuracy"] = self._evaluator.summary(self.rooms)
+        if learned:
+            # Persist at each sample so restart does not erase hours of training.
+            await self._store.async_save({
+                "schema_version": 2,
+                "learners": {room: learner.as_dict() for room, learner in self.learners.items()},
+                # Optional keys; older 0.2.x readers ignore them and keep working.
+                "retired_learners": self._retired,
+                "evaluation": self._evaluator.as_dict(),
+                "on_samples": self._on_samples, "off_samples": self._off_samples,
+                "heater_entity": self.config[CONF_HEATER],
+                "outside_source": active_sources["temperature"],
+            })
         weights = self.config.get(CONF_WEIGHTS, {})
         try:
             weights = {r: float(weights.get(r, 1.0)) for r in self.rooms}
@@ -381,6 +394,7 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
             "prediction_60m": output["predicted_60m"],
             "prediction_end": output["predicted_end"],
             "plan_15min": output["plan_15min"], "objective": output["cost"],
+            "prediction_trajectory": output["trajectory_15min"],
             "model_quality_note": "5-minute in-sample MAE; not a validated 60-minute forecast error",
         })
         return status

custom_components/heating_mpc/evaluation.py

@@ -0,0 +1,161 @@
+"""Walk-forward evaluation of the thermal model against what actually happened.
+
+Pure Python; no Home Assistant imports.
+
+At every poll the readings, the outdoor temperature, the radiator-lag estimate
+and each room's model parameters are recorded. When 15/30/60/120 minutes have
+passed, each room is replayed from that earlier poll with the parameters it had
+then (so later data is never used), the outdoor temperature held at its starting
+value, and the radiator lag that actually followed from the relay's real ON/OFF
+history. The replayed temperature is compared with the reading now.
+
+This scores the thermal model, not the MPC's proposed plan: Generic Thermostat
+controls the relay, so the plan is never executed. A persistence forecast (no
+change) is scored alongside as the baseline to beat.
+"""
+from __future__ import annotations
+
+from collections import deque
+from datetime import datetime, timedelta
+from math import isfinite
+from typing import Any
+
+from .physics import Learner
+
+HORIZONS = (15, 30, 60, 120)  # minutes
+WINDOW_DAYS = 7
+STEP_SECONDS = 300
+
+
+def replay(temp: float, outside: float, theta: tuple[float, float, float],
+           steps: list[tuple[float, float]]) -> float:
+    """Advance one room through (end-of-step lag, step seconds) with fixed parameters."""
+    model = Learner(theta=list(theta))
+    for lag, seconds in steps:
+        # Same 5-minute model and jitter normalization as learning.
+        temp += model.predict_delta(temp, outside, lag) * seconds / STEP_SECONDS
+    return temp
+
+
+class ForecastEvaluator:
+    """Score replayed forecasts per room and horizon over the last WINDOW_DAYS days."""
+
+    def __init__(self, max_lag_uncertainty: float) -> None:
+        self._max_lag_uncertainty = max_lag_uncertainty
+        self._history: deque[dict[str, Any]] = deque(maxlen=HORIZONS[-1] * 60 // STEP_SECONDS + 1)
+        # UTC day -> horizon -> room -> [count, sum error, sum |error|, sum |persistence error|]
+        self._days: dict[str, dict[str, dict[str, list[float]]]] = {}
+
+    def reset_window(self) -> None:
+        """Inputs were missing: no forecast may be scored across this poll."""
+        self._history.clear()
+
+    def observe(self, time: datetime, temps: dict[str, float], outside: float, lag: float,
+                lag_uncertainty: float, theta: dict[str, tuple[float, float, float]]) -> None:
+        """Record a poll, then score every forecast that has reached its horizon."""
+        if self._history and not 270 <= (time - self._history[-1]["time"]).total_seconds() <= 330:
+            self._history.clear()  # irregular poll spacing: replays would be unreliable
+        self._history.append({"time": time, "temps": dict(temps), "outside": outside, "lag": lag,
+                              "lag_uncertainty": lag_uncertainty, "theta": dict(theta)})
+        self._score()
+        self._prune(time)
+
+    def _score(self) -> None:
+        history = list(self._history)
+        latest = history[-1]
+        day = latest["time"].date().isoformat()
+        for horizon in HORIZONS:
+            steps = horizon * 60 // STEP_SECONDS
+            if len(history) < steps + 1:
+                continue
+            path = history[-(steps + 1):]
+            # The actual heat input must be known along the way.
+            if any(poll["lag_uncertainty"] > self._max_lag_uncertainty for poll in path[1:]):
+                continue
+            origin = path[0]
+            inputs = [(poll["lag"], (poll["time"] - before["time"]).total_seconds())
+                      for before, poll in zip(path, path[1:])]
+            for room, actual in latest["temps"].items():
+                start, theta = origin["temps"].get(room), origin["theta"].get(room)
+                if start is None or theta is None:
+                    continue
+                error = replay(start, origin["outside"], theta, inputs) - actual
+                bucket = self._days.setdefault(day, {}).setdefault(str(horizon), {}).setdefault(
+                    room, [0, 0.0, 0.0, 0.0])
+                bucket[0] += 1
+                bucket[1] += error
+                bucket[2] += abs(error)
+                bucket[3] += abs(start - actual)
+
+    def _prune(self, time: datetime) -> None:
+        oldest = (time.date() - timedelta(days=WINDOW_DAYS - 1)).isoformat()
+        for day in [day for day in self._days if day < oldest]:
+            del self._days[day]
+
+    def summary(self, rooms: list[str]) -> dict[int, dict[str, Any]]:
+        """Error statistics per horizon, overall and per configured room."""
+        result = {}
+        for horizon in HORIZONS:
+            per_room: dict[str, list[float]] = {}
+            for day in self._days.values():
+                for room, values in day.get(str(horizon), {}).items():
+                    if room in rooms:
+                        totals = per_room.setdefault(room, [0, 0.0, 0.0, 0.0])
+                        for i, value in enumerate(values):
+                            totals[i] += value
+            overall = [sum(values[i] for values in per_room.values()) for i in range(4)]
+            result[horizon] = {**_stats(overall),
+                               "rooms": {room: _stats(values) for room, values in per_room.items()}}
+        return result
+
+    def as_dict(self) -> dict[str, Any]:
+        # Rounded copy: shorter JSON on every save; in-memory sums stay exact.
+        return {"days": {day: {horizon: {room: [v[0], *(round(x, 6) for x in v[1:])]
+                                         for room, v in by_room.items()}
+                               for horizon, by_room in by_horizon.items()}
+                         for day, by_horizon in self._days.items()}}
+
+    def load(self, data: Any) -> None:
+        """Restore stored statistics, keeping only well-formed entries."""
+        self._days = {}
+        days = data.get("days") if isinstance(data, dict) else None
+        if not isinstance(days, dict):
+            return
+        horizons = {str(horizon) for horizon in HORIZONS}
+        for day, by_horizon in days.items():
+            try:
+                datetime.strptime(day, "%Y-%m-%d")
+            except (TypeError, ValueError):
+                continue
+            if not isinstance(by_horizon, dict):
+                continue
+            for horizon, by_room in by_horizon.items():
+                if horizon not in horizons or not isinstance(by_room, dict):
+                    continue
+                for room, values in by_room.items():
+                    if _valid_bucket(room, values):
+                        self._days.setdefault(day, {}).setdefault(horizon, {})[room] = [
+                            values[0], float(values[1]), float(values[2]), float(values[3])]
+
+
+def _valid_bucket(room: Any, values: Any) -> bool:
+    if not isinstance(room, str) or not isinstance(values, list) or len(values) != 4:
+        return False
+    count, *sums = values
+    return (type(count) is int and count >= 0
+            and all(type(v) in (int, float) and isfinite(v) for v in sums)
+            and sums[1] >= 0 and sums[2] >= 0)
+
+
+def _stats(values: list[float]) -> dict[str, Any]:
+    count, error, absolute, persistence = values
+    if not count:
+        return {"samples": 0, "mae": None, "bias": None, "persistence_mae": None, "skill": None}
+    return {
+        "samples": int(count),
+        "mae": round(absolute / count, 3),
+        "bias": round(error / count, 3),
+        "persistence_mae": round(persistence / count, 3),
+        # Share of the persistence error removed by the model (1 = perfect, < 0 = worse).
+        "skill": round(1 - absolute / persistence, 3) if persistence > 0 else None,
+    }

custom_components/heating_mpc/physics.py

@@ -278,10 +278,19 @@ def optimize(
     if result is None:
         raise RuntimeError("No feasible schedule")
     plan, temps_60m, temps_end = result
+    # Re-simulate the chosen plan for its trajectory, keeping the search loop lean.
+    trajectory: dict[str, list[float]] = {room: [] for room in temps}
+    state, lag = dict(temps), effective_heat
+    for bit in plan:
+        for _ in range(nsteps_per_block):
+            state, lag = simulate_step(state, outside, lag, bit, learners)
+        for room, value in state.items():
+            trajectory[room].append(round(value, 3))
     return {
         "proposal": bool(plan[0]),
         "plan_15min": [int(v) for v in plan],
         "predicted_60m": {r: round(t, 3) for r, t in temps_60m.items()},
         "predicted_end": {r: round(t, 3) for r, t in temps_end.items()},
+        "trajectory_15min": trajectory,
         "cost": round(best, 4),
     }

custom_components/heating_mpc/sensor.py

@@ -1,17 +1,19 @@
 """Diagnostics and predicted temperatures exposed as HA sensors."""
 from __future__ import annotations
 
-from homeassistant.components.sensor import SensorEntity, SensorDeviceClass
+from homeassistant.components.sensor import SensorDeviceClass, SensorEntity, SensorStateClass
 from homeassistant.const import UnitOfTemperature
 from homeassistant.helpers.update_coordinator import CoordinatorEntity
 
 from .const import DOMAIN
+from .evaluation import HORIZONS, WINDOW_DAYS
 
 
 async def async_setup_entry(hass, entry, async_add_entities) -> None:
     coordinator = hass.data[DOMAIN][entry.entry_id]
     entities = [MPCStatus(coordinator, entry.entry_id), MPCWeightedTemp(coordinator, entry.entry_id), MPCOutdoorTemperature(coordinator, entry.entry_id)]
     entities += [MPCRoomPrediction(coordinator, entry.entry_id, room) for room in coordinator.rooms]
+    entities += [MPCForecastError(coordinator, entry.entry_id, horizon) for horizon in HORIZONS]
     async_add_entities(entities)
 
 
@@ -34,7 +36,9 @@ class MPCStatus(CoordinatorEntity, SensorEntity):
     @property
     def extra_state_attributes(self):
         data = self.coordinator.data
-        return {key: value for key, value in data.items() if key not in ("status", "prediction_60m", "prediction_end")}
+        # Forecast details live on their own entities.
+        return {key: value for key, value in data.items() if key not in (
+            "status", "prediction_60m", "prediction_end", "prediction_trajectory", "forecast_accuracy")}
 
 
 class MPCWeightedTemp(CoordinatorEntity, SensorEntity):
@@ -64,10 +68,22 @@ class MPCRoomPrediction(CoordinatorEntity, SensorEntity):
         self._attr_unique_id = f"{entry_id}_{room}_predicted_60min"
         self._attr_name = f"MPC {room} 60min Prediction"
 
+    # Change every poll; the 60-minute value itself is recorded.
+    _unrecorded_attributes = frozenset({"forecast_15min", "plan_15min"})
+
     @property
     def native_value(self):
         return self.coordinator.data.get("prediction_60m", {}).get(self.room)
 
+    @property
+    def extra_state_attributes(self):
+        data = self.coordinator.data
+        return {
+            "forecast_15min": data.get("prediction_trajectory", {}).get(self.room),
+            "plan_15min": data.get("plan_15min"),
+            "basis": "proposed MPC plan; counterfactual while Generic Thermostat controls the relay",
+        }
+
 class MPCOutdoorTemperature(CoordinatorEntity, SensorEntity):
     """Display chosen outdoor temperature and all currently valid alternatives."""
     _attr_has_entity_name = True
@@ -91,3 +107,43 @@ class MPCOutdoorTemperature(CoordinatorEntity, SensorEntity):
             "active_source": data.get("weather_active_sources", {}).get("temperature"),
             "alternatives": data.get("weather_comparison", {}).get("temperature", {}),
         }
+
+
+class MPCForecastError(CoordinatorEntity, SensorEntity):
+    """Mean absolute error of the thermal model at one horizon, replayed against reality."""
+    _attr_has_entity_name = True
+    _attr_device_class = SensorDeviceClass.TEMPERATURE_DELTA
+    _attr_native_unit_of_measurement = UnitOfTemperature.CELSIUS
+    _attr_state_class = SensorStateClass.MEASUREMENT
+    _attr_suggested_display_precision = 2
+    _attr_icon = "mdi:target"
+    _unrecorded_attributes = frozenset({"rooms"})
+
+    def __init__(self, coordinator, entry_id, horizon):
+        super().__init__(coordinator)
+        self.horizon = horizon
+        self._attr_unique_id = f"{entry_id}_forecast_error_{horizon}min"
+        self._attr_name = f"Heating MPC Forecast Error {horizon} min"
+
+    @property
+    def _accuracy(self):
+        return self.coordinator.data.get("forecast_accuracy", {}).get(self.horizon, {})
+
+    @property
+    def native_value(self):
+        return self._accuracy.get("mae")
+
+    @property
+    def extra_state_attributes(self):
+        accuracy = self._accuracy
+        return {
+            "samples": accuracy.get("samples", 0),
+            "bias": accuracy.get("bias"),
+            "persistence_mae": accuracy.get("persistence_mae"),
+            "skill": accuracy.get("skill"),
+            "rooms": accuracy.get("rooms", {}),
+            "window_days": WINDOW_DAYS,
+            "unit_of_attributes": "°C",
+            "basis": ("thermal model replayed with the actual relay history; outdoor "
+                      "temperature held at the start value"),
+        }

tests/ha/test_lifecycle.py

@@ -24,7 +24,7 @@ async def test_setup_reload_unload(hass, entry, inputs):
     assert coordinator.config_entry is entry
     registry = er.async_get(hass)
     entities = er.async_entries_for_config_entry(registry, entry.entry_id)
-    assert len(entities) == 5
+    assert len(entities) == 9  # 5 core + 4 forecast-error sensors
     ids = {e.unique_id: e.entity_id for e in entities}
     assert coordinator.data["status"] == "learning"
     status = hass.states.get(ids[f"{entry.entry_id}_status"])
@@ -55,7 +55,7 @@ async def test_options_reload_removes_only_old_room(hass, entry, inputs):
     assert hass.data[DOMAIN][entry.entry_id].rooms == ["sensor.new"]
     assert registry.async_get(old) is None
     assert registry.async_get_entity_id("sensor", DOMAIN, f"{entry.entry_id}_sensor.new_predicted_60min")
-    assert len(er.async_entries_for_config_entry(registry, entry.entry_id)) == 5
+    assert len(er.async_entries_for_config_entry(registry, entry.entry_id)) == 9
 
 
 async def test_storage_restore_and_remove(hass, entry, inputs):
@@ -69,6 +69,12 @@ async def test_storage_restore_and_remove(hass, entry, inputs):
     restored = hass.data[DOMAIN][entry.entry_id]
     assert restored.learners["sensor.room"].count == 250
     assert restored.data["status"] == "ready"
+    forecast_id = er.async_get(hass).async_get_entity_id(
+        "sensor", DOMAIN, f"{entry.entry_id}_sensor.room_predicted_60min")
+    forecast = hass.states.get(forecast_id)
+    assert len(forecast.attributes["forecast_15min"]) == len(forecast.attributes["plan_15min"]) == 4
+    assert forecast.attributes["forecast_15min"][3] == float(forecast.state)
+    assert {"forecast_15min", "plan_15min"} <= forecast.state_info["unrecorded_attributes"]
     assert await hass.config_entries.async_remove(entry.entry_id)
     await hass.async_block_till_done()
     assert await Store(hass, 1, f"{DOMAIN}.{entry.entry_id}").async_load() is None
@@ -210,3 +216,25 @@ async def test_removed_room_model_survives_options_round_trip(hass, entry, input
         entry, options={"room_entities": ["sensor.room", "sensor.new"]})
     await hass.async_block_till_done()
     assert hass.data[DOMAIN][entry.entry_id].learners["sensor.room"].count == 250
+
+
+async def test_forecast_error_sensors_score_scheduled_polls(hass, entry, inputs, freezer):
+    coordinator = await setup(hass, entry)
+    registry = er.async_get(hass)
+    ids = {horizon: registry.async_get_entity_id("sensor", DOMAIN, f"{entry.entry_id}_forecast_error_{horizon}min")
+           for horizon in (15, 30, 60, 120)}
+    assert all(hass.states.get(entity).state == "unknown" for entity in ids.values())
+    future = dt_util.utcnow()
+    for _ in range(3):  # 15 minutes of regular polls
+        future += timedelta(minutes=5, seconds=1)
+        freezer.move_to(future)
+        async_fire_time_changed(hass, future)
+        await hass.async_block_till_done()
+    state = hass.states.get(ids[15])
+    assert coordinator.data["forecast_accuracy"][15]["samples"] == 1
+    assert float(state.state) == coordinator.data["forecast_accuracy"][15]["mae"] > 0
+    assert state.attributes["samples"] == 1
+    assert state.attributes["rooms"]["sensor.room"]["samples"] == 1
+    assert state.attributes["unit_of_measurement"] == "°C"
+    assert "rooms" in state.state_info["unrecorded_attributes"]
+    assert hass.states.get(ids[30]).state == "unknown"

tests/test_coordinator.py

@@ -423,7 +423,7 @@ class BoundaryTests(unittest.IsolatedAsyncioTestCase):
             for entity, value in (('sensor.room', temp), ('sensor.outside', outside)):
                 self.states[entity].state = f'{value:.5f}'
                 self.states[entity].last_reported = self.now
-            await self.sample()
+            data = await self.sample()
             command = temp < 20.7 or (on and temp < 21.3)
             if command != on:  # relay switches 120 s after this poll
                 self.switch('on' if command else 'off', 180)
@@ -436,6 +436,11 @@ class BoundaryTests(unittest.IsolatedAsyncioTestCase):
         self.assertGreater(model.on_count, 100)
         self.assertAlmostEqual(model.theta[1], b, delta=0.1 * b)
         self.assertAlmostEqual(model.theta[0], a, delta=0.25 * a)
+        # Walk-forward replay with the real relay history beats persistence clearly.
+        accuracy = data['forecast_accuracy'][60]
+        self.assertGreater(accuracy['samples'], 1400)
+        self.assertLess(accuracy['mae'], 0.06)
+        self.assertGreater(accuracy['skill'], 0.85)
 
     async def test_several_relay_changes_between_polls_are_integrated(self):
         advance = coordinator.advance_lag
@@ -654,3 +659,30 @@ class BoundaryTests(unittest.IsolatedAsyncioTestCase):
             self.assertIs(type(data['count']), int)
             float(data['error_sum'])
             [float(v) for v in data['theta'] + [v for row in data['p'] for v in row]]
+
+    async def test_forecasts_are_scored_saved_and_restored(self):
+        for _ in range(5):
+            self.fresh('sensor.room', 'sensor.outside')
+            data = await self.sample()
+        accuracy = data['forecast_accuracy']
+        self.assertEqual((accuracy[15]['samples'], accuracy[30]['samples'], accuracy[120]['samples']),
+                         (2, 0, 0))
+        # The default parameters predict cooling while the room held 20 °C.
+        self.assertLess(accuracy[15]['rooms']['sensor.room']['bias'], 0)
+        self.assertEqual(accuracy[15]['persistence_mae'], 0)
+        payload = self.c._store.async_save.call_args.args[0]
+        restored = coordinator.HeatingMPCCoordinator(self.hass, self.entry)
+        restored._store.async_load.return_value = payload
+        await restored.async_restore_model()
+        self.assertEqual(restored._evaluator.summary(['sensor.room']), self.c._evaluator.summary(['sensor.room']))
+        # A poll without a usable relay state breaks the replay window.
+        self.states['switch.boiler'].state = 'unavailable'
+        self.fresh('sensor.room', 'sensor.outside')
+        await self.sample()
+        self.states['switch.boiler'].state = 'off'
+        for _ in range(3):
+            self.fresh('sensor.room', 'sensor.outside')
+            data = await self.sample()
+        self.assertEqual(data['forecast_accuracy'][15]['samples'], 2)
+        self.fresh('sensor.room', 'sensor.outside')
+        self.assertEqual((await self.sample())['forecast_accuracy'][15]['samples'], 3)

tests/test_evaluation.py

@@ -0,0 +1,97 @@
+"""Walk-forward forecast evaluation, without Home Assistant."""
+import sys
+import types
+import unittest
+from datetime import datetime, timedelta, timezone
+from pathlib import Path
+
+COMPONENTS = Path(__file__).resolve().parents[1] / 'custom_components' / 'heating_mpc'
+# Fake only the package container, to avoid importing Home Assistant from __init__.
+sys.modules.setdefault('heating_mpc', types.ModuleType('heating_mpc'))
+sys.modules['heating_mpc'].__path__ = [str(COMPONENTS)]
+from heating_mpc import evaluation as ev  # noqa: E402
+
+THETA = (0.004, 0.12, 0.0005)
+START = datetime(2026, 10, 26, tzinfo=timezone.utc)
+
+
+def plant(polls, theta=THETA, model_theta=None, outside=5.0, lag_at=lambda k: 1.0 if k % 30 < 12 else 0.0,
+          evaluator=None, uncertainty_at=lambda k: 0.0, start=START, step=300):
+    """Feed polls of a room that follows the 5-minute model exactly with `theta`."""
+    a, b, c = theta
+    evaluator = evaluator or ev.ForecastEvaluator(0.1)
+    temp = 20.0
+    for k in range(polls):
+        lag = lag_at(k)
+        if k:
+            temp += a * (outside - temp) + b * lag + c
+        evaluator.observe(start + timedelta(seconds=k * step), {'sensor.room': temp}, outside, lag,
+                          uncertainty_at(k), {'sensor.room': model_theta or theta})
+    return evaluator
+
+
+class EvaluationTests(unittest.TestCase):
+    def test_exact_model_scores_zero_error_and_beats_persistence(self):
+        summary = plant(60).summary(['sensor.room'])
+        for horizon, expected_samples in ((15, 57), (30, 54), (60, 48), (120, 36)):
+            with self.subTest(horizon=horizon):
+                stats = summary[horizon]
+                self.assertEqual(stats['samples'], expected_samples)
+                self.assertEqual(stats['rooms']['sensor.room']['samples'], expected_samples)
+                self.assertAlmostEqual(stats['mae'], 0, places=6)
+                self.assertGreater(stats['persistence_mae'], 0.05)
+                self.assertEqual(stats['skill'], 1)
+
+    def test_parameters_from_the_forecast_start_are_used(self):
+        # A model without heat gain under-predicts warming: negative bias.
+        stats = plant(10, model_theta=(0.004, 0.0, 0.0005), lag_at=lambda k: 1.0).summary(['sensor.room'])
+        self.assertLess(stats[15]['bias'], -0.3)
+        self.assertEqual(stats[15]['bias'], -stats[15]['mae'])
+
+    def test_replay_matches_hand_calculation(self):
+        temp = 20.0
+        for _ in range(3):
+            temp += 0.004 * (5 - temp) + 0.12 * 0.5 + 0.0005
+        self.assertAlmostEqual(ev.replay(20.0, 5.0, THETA, [(0.5, 300)] * 3), temp)
+        # Jitter is normalized like learning: 330 s counts as 1.1 steps.
+        self.assertAlmostEqual(ev.replay(20.0, 5.0, (0, 0.1, 0), [(1.0, 330)]), 20.11)
+
+    def test_breaks_uncertain_lag_and_missing_rooms_are_not_scored(self):
+        evaluator = plant(4)
+        self.assertEqual(evaluator.summary(['sensor.room'])[15]['samples'], 1)
+        evaluator.reset_window()
+        plant(3, evaluator=evaluator, start=START + timedelta(minutes=20))
+        self.assertEqual(evaluator.summary(['sensor.room'])[15]['samples'], 1)
+        # Irregular spacing restarts the window as well.
+        jittery = plant(8, step=400)
+        self.assertEqual(jittery.summary(['sensor.room'])[15]['samples'], 0)
+        uncertain = plant(8, uncertainty_at=lambda k: 0.5 if k == 4 else 0.0)
+        self.assertEqual(uncertain.summary(['sensor.room'])[15]['samples'], 2)  # k = 3 and 7 only
+        self.assertEqual(plant(8).summary(['sensor.other'])[15]['samples'], 0)
+
+    def test_window_keeps_seven_days(self):
+        evaluator = plant(10)
+        before = evaluator.summary(['sensor.room'])[15]['samples']
+        plant(10, evaluator=evaluator, start=START + timedelta(days=7))
+        self.assertEqual(evaluator.summary(['sensor.room'])[15]['samples'], 7)
+        self.assertEqual(before, 7)
+
+    def test_saved_statistics_round_trip_and_reject_garbage(self):
+        evaluator = plant(10)
+        restored = ev.ForecastEvaluator(0.1)
+        restored.load(evaluator.as_dict())
+        self.assertEqual(restored.summary(['sensor.room']), evaluator.summary(['sensor.room']))
+        day = START.date().isoformat()
+        restored.load({'days': {
+            'not-a-day': {'15': {'sensor.room': [1, 0.0, 0.1, 0.1]}},
+            day: {'45': {'sensor.room': [1, 0.0, 0.1, 0.1]},
+                  '15': {'sensor.a': [1, 0.0, -0.1, 0.1], 'sensor.b': [True, 0, 0, 0],
+                         'sensor.c': [1, float('nan'), 0, 0], 'sensor.d': [2, 0.1, 0.2, 0.4]}}}})
+        self.assertEqual(restored.as_dict(), {'days': {day: {'15': {'sensor.d': [2, 0.1, 0.2, 0.4]}}}})
+        for garbage in (None, [], {'days': []}):
+            restored.load(garbage)
+            self.assertEqual(restored.as_dict(), {'days': {}})
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_physics.py

@@ -76,7 +76,9 @@ class ModelTests(unittest.TestCase):
                 temps, lag = m.simulate_step(temps, 5, lag, bit, learners)
             if index == 3:
                 self.assertEqual(out['predicted_60m']['r'], round(temps['r'], 3))
+            self.assertEqual(out['trajectory_15min']['r'][index], round(temps['r'], 3))
         self.assertEqual(out['predicted_end']['r'], round(temps['r'], 3))
+        self.assertEqual(len(out['trajectory_15min']['r']), 6)
 
     def test_jitter_normalized_and_rejection_counted(self):
         nominal, delayed = m.Learner(), m.Learner()