diff --git a/agents/codex-961.md b/agents/codex-961.md
new file mode 100644
index 00000000..12a56418
--- /dev/null
+++ b/agents/codex-961.md
@@ -0,0 +1 @@
+<!-- bootstrap for codex on issue #961 -->
diff --git a/pa_core/sim/paths.py b/pa_core/sim/paths.py
index 38788ad1..f4a00aa9 100644
--- a/pa_core/sim/paths.py
+++ b/pa_core/sim/paths.py
@@ -359,10 +359,14 @@ def draw_returns(
     n_months: int,
     n_sim: int,
     params: Dict[str, Any],
+    seed: Optional[int] = None,
     rng: Optional[GeneratorLike] = None,
     shocks: Optional[Dict[str, Any]] = None,
 ) -> tuple[npt.NDArray[Any], npt.NDArray[Any], npt.NDArray[Any], npt.NDArray[Any]]:
-    """Vectorised draw of monthly returns for (beta, H, E, M)."""
+    """Vectorised draw of monthly returns for (beta, H, E, M).
+
+    ``seed`` is used to create a per-run generator when ``rng`` is not supplied.
+    """
     _assert_canonical_params(params)
     distribution = params.get("return_distribution", "normal")
     dist_overrides = (
@@ -439,7 +443,7 @@ def draw_returns(
                 sims = μ + shocks_out * σ
     else:
         if rng is None:
-            rng = spawn_rngs(None, 1)[0]
+            rng = spawn_rngs(seed, 1)[0]
         assert rng is not None
         if all(dist == "normal" for dist in distributions):
             Σ = corr * (σ[:, None] * σ[None, :])
@@ -491,6 +495,7 @@ def draw_joint_returns(
     n_months: int,
     n_sim: int,
     params: Dict[str, Any],
+    seed: Optional[int] = None,
     rng: Optional[GeneratorLike] = None,
     shocks: Optional[Dict[str, Any]] = None,
     regime_paths: Optional[npt.NDArray[Any]] = None,
@@ -501,11 +506,14 @@ def draw_joint_returns(
     When ``regime_params`` is provided, returns are drawn from regime-specific
     parameters based on ``regime_paths``.
     """
+    if rng is None and seed is not None:
+        rng = spawn_rngs(seed, 1)[0]
     if regime_params is None:
         return draw_returns(
             n_months=n_months,
             n_sim=n_sim,
             params=params,
+            seed=seed,
             rng=rng,
             shocks=shocks,
         )
@@ -534,6 +542,7 @@ def draw_joint_returns(
             n_months=n_months,
             n_sim=n_sim,
             params=regime,
+            seed=seed,
             rng=rng,
         )
         mask = regime_paths == regime_idx
diff --git a/pa_core/sim/regimes.py b/pa_core/sim/regimes.py
index 093f235b..a68f8d33 100644
--- a/pa_core/sim/regimes.py
+++ b/pa_core/sim/regimes.py
@@ -103,9 +103,13 @@ def simulate_regime_paths(
     n_months: int,
     transition: Sequence[Sequence[float]],
     start_state: int,
+    seed: int | None = None,
     rng: Any | None = None,
 ) -> npt.NDArray[Any]:
-    """Simulate regime paths using a Markov transition matrix."""
+    """Simulate regime paths using a Markov transition matrix.
+
+    ``seed`` is used to create a per-run generator when ``rng`` is not supplied.
+    """
     if n_sim <= 0 or n_months <= 0:
         raise ValueError("n_sim and n_months must be positive")
     transition_mat = np.asarray(transition, dtype=float)
@@ -115,7 +119,7 @@ def simulate_regime_paths(
     if not 0 <= start_state < n_regimes:
         raise ValueError("start_state must be within regime index range")
     if rng is None:
-        rng = spawn_rngs(None, 1)[0]
+        rng = spawn_rngs(seed, 1)[0]
 
     paths = np.empty((n_sim, n_months), dtype=int)
     paths[:, 0] = start_state
diff --git a/tests/test_regime_switching.py b/tests/test_regime_switching.py
index 3ae93c09..5342cbee 100644
--- a/tests/test_regime_switching.py
+++ b/tests/test_regime_switching.py
@@ -110,3 +110,72 @@ def test_regime_switching_increases_corr_and_vol() -> None:
 
     assert stress_corr > calm_corr
     assert stress_vol > calm_vol
+
+
+def test_regime_switching_seed_is_deterministic() -> None:
+    cfg = ModelConfig(
+        N_SIMULATIONS=128,
+        N_MONTHS=12,
+        financing_mode="broadcast",
+        return_unit="monthly",
+        sigma_H=0.02,
+        sigma_E=0.02,
+        sigma_M=0.02,
+        rho_idx_H=0.1,
+        rho_idx_E=0.1,
+        rho_idx_M=0.1,
+        rho_H_E=0.2,
+        rho_H_M=0.2,
+        rho_E_M=0.2,
+        regimes=[
+            RegimeConfig(name="calm"),
+            RegimeConfig(
+                name="stress",
+                idx_sigma_multiplier=1.5,
+                sigma_H=0.04,
+                sigma_E=0.04,
+                sigma_M=0.04,
+                rho_idx_H=0.7,
+                rho_idx_E=0.7,
+                rho_idx_M=0.7,
+                rho_H_E=0.75,
+                rho_H_M=0.75,
+                rho_E_M=0.75,
+            ),
+        ],
+        regime_transition=[[0.8, 0.2], [0.2, 0.8]],
+        regime_start="calm",
+    )
+    params, _labels = build_regime_draw_params(
+        cfg,
+        mu_idx=0.0,
+        idx_sigma=0.015,
+        n_samples=120,
+    )
+
+    def _run(seed: int) -> tuple[np.ndarray, np.ndarray]:
+        paths = simulate_regime_paths(
+            n_sim=cfg.N_SIMULATIONS,
+            n_months=cfg.N_MONTHS,
+            transition=cfg.regime_transition or [],
+            start_state=resolve_regime_start(cfg),
+            seed=seed,
+        )
+        _r_beta, r_H, r_E, _r_M = draw_joint_returns(
+            n_months=cfg.N_MONTHS,
+            n_sim=cfg.N_SIMULATIONS,
+            params=params[0],
+            seed=seed + 1,
+            regime_paths=paths,
+            regime_params=params,
+        )
+        return paths, np.stack([r_H, r_E])
+
+    paths_a, returns_a = _run(101)
+    paths_b, returns_b = _run(202)
+    paths_c, returns_c = _run(101)
+
+    assert np.array_equal(paths_a, paths_c)
+    assert np.array_equal(returns_a, returns_c)
+    assert not np.array_equal(paths_a, paths_b)
+    assert not np.array_equal(returns_a, returns_b)