diff --git a/README.md b/README.md
index 4ccae6d..eeeadf5 100644
--- a/README.md
+++ b/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
diff --git a/custom_components/heating_mpc/coordinator.py b/custom_components/heating_mpc/coordinator.py
index bb076bc..06549bd 100644
--- a/custom_components/heating_mpc/coordinator.py
+++ b/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.
-                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"],
-                })
+                learned = True
 
         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
diff --git a/custom_components/heating_mpc/evaluation.py b/custom_components/heating_mpc/evaluation.py
new file mode 100644
index 0000000..cfe42b4
--- /dev/null
+++ b/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,
+    }
diff --git a/custom_components/heating_mpc/physics.py b/custom_components/heating_mpc/physics.py
index a93f0eb..9b073d3 100644
--- a/custom_components/heating_mpc/physics.py
+++ b/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),
     }
diff --git a/custom_components/heating_mpc/sensor.py b/custom_components/heating_mpc/sensor.py
index 7456f93..e2522a1 100644
--- a/custom_components/heating_mpc/sensor.py
+++ b/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"),
+        }
diff --git a/tests/ha/test_lifecycle.py b/tests/ha/test_lifecycle.py
index 236edab..f98a16b 100644
--- a/tests/ha/test_lifecycle.py
+++ b/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"
diff --git a/tests/test_coordinator.py b/tests/test_coordinator.py
index 80dfbeb..a74e2f6 100644
--- a/tests/test_coordinator.py
+++ b/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)
diff --git a/tests/test_evaluation.py b/tests/test_evaluation.py
new file mode 100644
index 0000000..a6eaf05
--- /dev/null
+++ b/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()
diff --git a/tests/test_physics.py b/tests/test_physics.py
index ecb873a..1ce14fa 100644
--- a/tests/test_physics.py
+++ b/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()