diff --git a/python/statsforecast/theta.py b/python/statsforecast/theta.py
index 35820ca25..f688d6aa0 100644
--- a/python/statsforecast/theta.py
+++ b/python/statsforecast/theta.py
@@ -272,6 +272,7 @@ def compute_pi_samples(
     error_distribution="normal",
     error_params=None,
     residuals=None,
+    modeltype=None,
 ):
     """
     Compute prediction interval samples for Theta model.
@@ -292,6 +293,15 @@ def compute_pi_samples(
     # states: level, meany, An, Bn, mu
     smoothed, _, A, B, _ = states[-1]
 
+    # The drift terms A and B only evolve for the dynamic models (DSTM/DOTM);
+    # for STM/OTM they are held constant, matching the update() reference in
+    # src/theta.cpp. modeltype=None keeps the previous behaviour for any
+    # external caller that does not pass it.
+    update_ab = modeltype is None or modeltype in (
+        _theta.ModelType.DSTM,
+        _theta.ModelType.DOTM,
+    )
+
     rng = np.random.default_rng(seed)
 
     for i in range(n, n + h):
@@ -309,9 +319,12 @@ def compute_pi_samples(
         )
         samples[i - n] += errors
         smoothed = alpha * samples[i - n] + (1 - alpha) * smoothed
-        mean_y = (i * mean_y + samples[i - n]) / (i + 1)
-        B = ((i - 1) * B + 6 * (samples[i - n] - mean_y) / (i + 1)) / (i + 2)
-        A = mean_y - B * (i + 2) / 2
+        new_mean_y = (i * mean_y + samples[i - n]) / (i + 1)
+        if update_ab:
+            # B uses the pre-update mean_y, matching src/theta.cpp.
+            B = ((i - 1) * B + 6 * (samples[i - n] - mean_y) / (i + 1)) / (i + 2)
+            A = new_mean_y - B * (i + 2) / 2
+        mean_y = new_mean_y
     return samples
 
 
@@ -375,6 +388,7 @@ def simulate_theta(
         error_distribution=error_distribution,
         error_params=error_params,
         residuals=residuals,
+        modeltype=switch_theta(model["modeltype"]),
     )
 
     res = samples.T
@@ -428,6 +442,7 @@ def forecast_theta(obj, h, level=None):
             error_distribution=dist,
             error_params=error_params_from_model(obj),
             residuals=residuals_tail,
+            modeltype=switch_theta(obj["modeltype"]),
         )
         for lv in level:
             min_q = (100 - lv) / 200
diff --git a/tests/test_theta.py b/tests/test_theta.py
index 7f153b756..bbc5bb148 100644
--- a/tests/test_theta.py
+++ b/tests/test_theta.py
@@ -5,6 +5,7 @@
 from statsforecast._lib import theta as _theta
 from statsforecast.theta import (
     auto_theta,
+    compute_pi_samples,
     forecast_theta,
     forward_theta,
     initparamtheta,
@@ -517,3 +518,37 @@ def test_theta_short_series_no_nan_intervals():
     assert not np.isnan(out["lo-95"]).any(), "NaN in lower prediction interval"
     assert not np.isnan(out["hi-95"]).any(), "NaN in upper prediction interval"
     assert not np.isnan(out["mean"]).any(), "NaN in point forecast"
+
+
+def test_theta_pi_stays_centered_for_static_models():
+    """The Theta prediction-interval simulation must stay centered on the analytic
+    point forecast.
+
+    For STM/OTM the drift terms ``A`` and ``B`` are held constant (matching the
+    ``update()`` reference in ``src/theta.cpp``); only DSTM/DOTM evolve them.
+    ``compute_pi_samples`` previously evolved them for every model type, which
+    miscentered the STM/OTM intervals with a bias that grew over the horizon.
+    """
+    rng = np.random.default_rng(0)
+    y = (
+        np.arange(1, 73) + 15 * np.sin(np.arange(72) / 6.0) + rng.normal(0, 1.5, 72)
+    ).astype(np.float64)
+    h = 24
+    for mtype in ["STM", "OTM", "DSTM", "DOTM"]:
+        res = auto_theta(y, m=12, model=mtype)
+        mean = forecast_theta(res, h, level=None)["mean"]
+        samples = compute_pi_samples(
+            n=res["n"],
+            h=h,
+            states=res["states"],
+            sigma=1.0,
+            alpha=res["par"]["alpha"],
+            theta=res["par"]["theta"],
+            mean_y=res["mean_y"],
+            n_samples=200_000,
+            modeltype=switch_theta(res["modeltype"]),
+        )
+        # With many samples the Monte Carlo error is ~1e-3, so a persistent gap
+        # above 0.05 indicates the simulation has drifted off the point forecast.
+        drift = float(np.abs(samples.mean(axis=1) - mean).max())
+        assert drift < 0.05, f"{mtype} PI drifted from the point forecast by {drift:.3f}"