diff --git a/README.md b/README.md
index 71cde3dd..9b4948fe 100644
--- a/README.md
+++ b/README.md
@@ -142,7 +142,7 @@ Important boundaries:
 - Stochastic exact stream parity is not expected because Python paths use NumPy random generation.
 - Factor and class ordering should be passed explicitly when ordering matters.
 - Direct raw-factor `nns_m_reg(..., factor_2_dummy=True)` is intentionally guarded. Use `prepare_factor_predictors(...)` before `nns_m_reg(...)`.
-- Compute functions' `plot` arguments are ignored and data is returned instead; visual plotting is a separate API in `nns.plotting`, color/element-faithful to R but not pixel-diffed.
+- Compute functions still return values, not figures; passing `plot=True` (where R has it) additionally renders a Matplotlib figure as a side effect via the `nns.plotting` layer, which is color/element-faithful to R but not pixel-diffed. The plot functions can also be called directly on a computed result.
 
 See [behavior conventions](docs/conventions.md) for detailed compatibility notes.
 
diff --git a/docs/plot_parity_policy.md b/docs/plot_parity_policy.md
index 4ac91129..d4e29555 100644
--- a/docs/plot_parity_policy.md
+++ b/docs/plot_parity_policy.md
@@ -50,9 +50,13 @@ colors and which element they sit on*, never rendered images.
   a plotting API, so Python parity calls pass the R `plot = FALSE` equivalent
   and assert only on returned values.
 
-When a ported function has an R `plot` argument, the Python API either omits the
-argument entirely or treats plotting as out of scope; only the value-bearing
-return is asserted in parity tests.
+When a ported function has an R `plot` argument, the Python function keeps its
+value-only **return** contract (parity asserts only on the returned value). As a
+side effect, passing `plot=True` renders a Matplotlib figure through the
+`nns.plotting` layer — `nns_reg`, `nns_m_reg`, `nns_arma`, `nns_arma_optim`,
+`nns_cdf`, and `nns_seas` are wired this way (plus `residual_plot=True` for the
+regression functions). The figures are color/element-faithful but never
+pixel-compared, and computation with the default `plot=False` opens no figure.
 
 ## Inventory of committed graphics artifacts
 
diff --git a/src/nns/arma.py b/src/nns/arma.py
index e3fb50f6..5b3501e6 100644
--- a/src/nns/arma.py
+++ b/src/nns/arma.py
@@ -32,7 +32,7 @@ def nns_arma_optim(
     plot: bool = False,
 ) -> dict[str, Any]:
     """Optimize seasonal factors for :func:`nns_arma` like R's ``NNS.ARMA.optim``."""
-    del ncores, print_trace, plot
+    del ncores, print_trace
 
     values = _as_variable(variable)
     original_values = values.copy()
@@ -326,7 +326,7 @@ def nns_arma_optim(
         lower_pi = np.maximum(0.0, lower_pi)
         upper_pi = np.maximum(0.0, upper_pi)
 
-    return {
+    result = {
         "periods": nns_periods,
         "weights": nns_weights,
         "obj.fn": nns_score,
@@ -339,6 +339,11 @@ def nns_arma_optim(
         "lower.pred.int": lower_pi,
         "upper.pred.int": upper_pi,
     }
+    if plot:
+        from nns.plotting.arma import plot_nns_arma_optim
+
+        plot_nns_arma_optim(result, original_values)
+    return result
 
 
 def nns_arma(
@@ -360,12 +365,22 @@ def nns_arma(
     random_seed: int | None = None,
 ) -> NDArray[np.float64] | dict[str, NDArray[np.float64]]:
     """Autoregressive NNS forecast matching R's installed NNS.ARMA behavior."""
-    del plot, seasonal_plot
+    del seasonal_plot
 
     horizon = int(h)
     if horizon < 1:
         raise ValueError("h must be a positive integer.")
     values = _as_variable(variable)
+
+    def _finish(
+        forecast: NDArray[np.float64] | dict[str, NDArray[np.float64]],
+    ) -> NDArray[np.float64] | dict[str, NDArray[np.float64]]:
+        if plot:
+            from nns.plotting.arma import plot_nns_arma
+
+            ts = int(training_set) if training_set is not None else int(values.size)
+            plot_nns_arma(forecast, values, training_set=ts)
+        return forecast
     if _is_numeric_seasonal(seasonal_factor) and dynamic:
         raise ValueError(
             'Hmmm...Seems you have "seasonal.factor" specified and "dynamic = TRUE".  '
@@ -389,12 +404,12 @@ def nns_arma(
 
     estimates = np.zeros(horizon, dtype=np.float64)
     if not _is_numeric_seasonal(seasonal_factor) and np.ptp(values) == 0.0:
-        return _with_prediction_intervals(
+        return _finish(_with_prediction_intervals(
             estimates,
             lin_residual=0.0,
             pred_int=pred_int,
             random_seed=random_seed,
-        )
+        ))
     lags, lag_weights = _resolve_lags_and_weights(
         values,
         seasonal_factor=seasonal_factor,
@@ -416,7 +431,7 @@ def nns_arma(
             method=method_l,
             shrink=shrink,
         )
-        return estimates
+        return _finish(estimates)
 
     current = values
     lin_regression_estimates = np.array([], dtype=np.float64)
@@ -490,12 +505,12 @@ def nns_arma(
         if not np.isfinite(lin_resid):
             lin_resid = 0.0
 
-    return _with_prediction_intervals(
+    return _finish(_with_prediction_intervals(
         estimates,
         lin_residual=lin_resid,
         pred_int=pred_int,
         random_seed=random_seed,
-    )
+    ))
 
 
 def _valid_arma_optim_seasonals(
diff --git a/src/nns/cdf.py b/src/nns/cdf.py
index a70748cf..fd0d803d 100644
--- a/src/nns/cdf.py
+++ b/src/nns/cdf.py
@@ -20,7 +20,6 @@ def nns_cdf(
     names: Sequence[str] | None = None,
 ) -> dict[str, object]:
     """Partial-moment CDF wrapper matching R's non-plotting NNS.CDF paths."""
-    del plot
     type_value = type.lower()
     if type_value not in {"cdf", "survival", "hazard", "cumulative hazard"}:
         raise ValueError("invalid type")
@@ -29,12 +28,28 @@ def nns_cdf(
     if values.ndim == 0:
         values = values.reshape(1)
     if values.ndim == 1 or (values.ndim == 2 and values.shape[1] == 1):
-        return _univariate_cdf(values.reshape(-1), float(degree), target, type_value)
+        result = _univariate_cdf(values.reshape(-1), float(degree), target, type_value)
+        if plot:
+            _render_cdf(result, target)
+        return result
     if values.ndim == 2:
+        # Multivariate CDF has no faithful single-Axes plot; plot is a no-op here.
         return _multivariate_cdf(values, float(degree), target, type_value, names)
     raise ValueError("variable must be a vector or 2D matrix.")
 
 
+def _render_cdf(result: dict[str, object], target: float | NDArray[np.float64] | None) -> None:
+    """Render the univariate NNS.CDF figure as a side effect of ``plot=True``."""
+    from nns.plotting.partial_moments import plot_nns_cdf
+
+    plot_target: float | None = None
+    if target is not None:
+        coords = np.asarray(target, dtype=np.float64).reshape(-1)
+        if coords.size == 1:
+            plot_target = float(coords[0])
+    plot_nns_cdf(result, target=plot_target)
+
+
 def _univariate_cdf(
     values: NDArray[np.float64],
     degree: float,
diff --git a/src/nns/multivariate_regression.py b/src/nns/multivariate_regression.py
index 7918b096..f699691e 100644
--- a/src/nns/multivariate_regression.py
+++ b/src/nns/multivariate_regression.py
@@ -39,8 +39,8 @@ def nns_m_reg(
     confidence_interval: float | None = None,
     class_levels: list[object] | None = None,
 ) -> MRegResult:
-    """Multivariate numeric regression matching R's non-plotting NNS.M.reg path."""
-    del plot, residual_plot, location, dist, return_values, plot_regions, ncores
+    """Multivariate numeric regression matching R's NNS.M.reg path."""
+    del location, dist, return_values, plot_regions, ncores
     type_value = _normalize_type(type)
     x_values, y_values = _validate_inputs(
         x,
@@ -107,7 +107,7 @@ def nns_m_reg(
         confidence_interval=confidence_interval,
     )
     r2 = _class_accuracy(y_values, fitted_y) if type_value == "class" else _r2(y_values, fitted_y)
-    return {
+    result: MRegResult = {
         "R2": r2,
         "rhs.partitions": _rhs_partitions_dict(reg_points_matrix),
         "RPM": _rpm_dict(rpm),
@@ -115,6 +115,41 @@ def nns_m_reg(
         "pred.int": pred_int,
         "Fitted.xy": fitted,
     }
+    if plot or residual_plot:
+        _render_m_reg(fitted)
+    return result
+
+
+def _render_m_reg(fitted: dict[str, NDArray[np.float64] | NDArray[np.str_]]) -> None:
+    """Render R's NNS.M.reg residual plot (Multivariate_Regression.R:367-377).
+
+    The multivariate plot output is the residual plot: actual ``y`` over the
+    observation index as ``steelblue`` open circles, fitted ``y.hat`` as a
+    ``red`` line, and a pink (alpha 0.375) confidence band when present.
+    """
+    from nns.plotting import palette
+    from nns.plotting._mpl import resolve_ax
+
+    y = np.asarray(fitted["y"], dtype=np.float64)
+    y_hat = np.asarray(fitted["y.hat"], dtype=np.float64)
+    if y.size == 0 or y.size != y_hat.size:
+        return
+    ax = resolve_ax(None)
+    index = np.arange(1, y.size + 1)
+    ax.scatter(index, y, facecolors="none", edgecolors="steelblue", marker="o")
+    ax.plot(index, y_hat, color="red", linewidth=2)
+    if "conf.int.pos" in fitted and "conf.int.neg" in fitted:
+        pos = np.asarray(fitted["conf.int.pos"], dtype=np.float64)
+        neg = np.asarray(fitted["conf.int.neg"], dtype=np.float64)
+        mask = np.isfinite(pos) & np.isfinite(neg)
+        if mask.any():
+            ax.fill_between(
+                index[mask], neg[mask], pos[mask],
+                color=palette.PINK, alpha=palette.CI_ALPHA_REG, linewidth=0.0,
+            )
+    ax.set_xlabel("Index")
+    ax.set_ylabel("y (blue)   y.hat (red)")
+    ax.set_title("NNS.M.reg Residual Plot")
 
 
 def _validate_inputs(
diff --git a/src/nns/regression.py b/src/nns/regression.py
index a955263c..711ee71c 100644
--- a/src/nns/regression.py
+++ b/src/nns/regression.py
@@ -57,10 +57,10 @@ def nns_reg(
     factor_levels: Sequence[object] | Sequence[Sequence[object] | None] | None = None,
 ) -> dict[str, Any]:
     """Univariate numeric port of R's NNS.reg."""
-    del return_values, plot, plot_regions, residual_plot, ncores
+    del return_values, ncores
 
     if dim_red_method is not None:
-        return _nns_reg_dimred(
+        result = _nns_reg_dimred(
             x,
             y,
             factor_2_dummy=factor_2_dummy,
@@ -80,6 +80,11 @@ def nns_reg(
             class_levels=class_levels,
             factor_levels=factor_levels,
         )
+        _maybe_render_reg(
+            result, plot=plot, plot_regions=plot_regions,
+            residual_plot=residual_plot, point_est=point_est,
+        )
+        return result
 
     type_value = _normalize_type(type)
     if type_value == "class":
@@ -101,7 +106,7 @@ def nns_reg(
         dispatch_n_best = n_best
         if type_value == "class" and dispatch_n_best is None:
             dispatch_n_best = 1
-        return nns_m_reg(
+        result = nns_m_reg(
             np.asarray(x_for_dispatch, dtype=np.float64),
             y_matrix_values,
             factor_2_dummy=False,
@@ -117,6 +122,11 @@ def nns_reg(
             confidence_interval=confidence_interval,
             class_levels=class_levels,
         )
+        _maybe_render_reg(
+            result, plot=plot, plot_regions=plot_regions,
+            residual_plot=residual_plot, point_est=point_est,
+        )
+        return result
 
     del tau, threshold, n_best, dist
     x_values, y_values = _validate_univariate_inputs(
@@ -137,7 +147,7 @@ def nns_reg(
     )
     noise = _validate_noise_reduction(noise_reduction)
     point_values = _as_point_est(point_for_dispatch)
-    return _nns_reg_univariate_core(
+    result = _nns_reg_univariate_core(
         x_values,
         y_values,
         order=order,
@@ -150,6 +160,47 @@ def nns_reg(
         equation=None,
         x_star=None,
     )
+    _maybe_render_reg(
+        result, plot=plot, plot_regions=plot_regions,
+        residual_plot=residual_plot, point_est=point_est,
+    )
+    return result
+
+
+def _maybe_render_reg(
+    result: dict[str, Any],
+    *,
+    plot: bool,
+    plot_regions: bool,
+    residual_plot: bool,
+    point_est: NDArray[np.float64] | float | None,
+) -> None:
+    """Render the NNS.reg figure(s) as a side effect when a plot flag is set.
+
+    Plotting is decoupled from computation: this only fires for the standard
+    univariate result (one that carries ``Fitted.xy``) so the value-only return
+    contract is unchanged. ``plot``/``plot_regions`` draw the regression figure;
+    ``residual_plot`` draws a residual scatter.
+    """
+    if not (plot or plot_regions or residual_plot):
+        return
+    fitted = result.get("Fitted.xy") if isinstance(result, dict) else None
+    if not isinstance(fitted, dict) or "x" not in fitted:
+        return
+    if plot or plot_regions:
+        from nns.plotting.regression import plot_nns_reg
+
+        plot_nns_reg(result, point_est=point_est)
+    if residual_plot:
+        from nns.plotting._mpl import resolve_ax
+
+        xs = np.asarray(fitted["x"], dtype=np.float64)
+        residuals = np.asarray(fitted.get("residuals", []), dtype=np.float64)
+        if residuals.size and residuals.size == xs.size:
+            ax = resolve_ax(None)
+            ax.scatter(xs, residuals, color="steelblue")
+            ax.axhline(0.0, color="red")
+            ax.set_title("NNS Residual Plot")
 
 
 def prepare_factor_predictors(
diff --git a/src/nns/seasonality.py b/src/nns/seasonality.py
index 6b89098f..907da374 100644
--- a/src/nns/seasonality.py
+++ b/src/nns/seasonality.py
@@ -19,8 +19,21 @@ def nns_seas(
     mod_only: bool = True,
     plot: bool = False,
 ) -> SeasonalityResult:
-    """Seasonality test matching R's NNS.seas non-plotting path."""
-    del plot
+    """Seasonality test matching R's NNS.seas; ``plot=True`` renders the figure."""
+    result = _nns_seas_compute(variable, modulo=modulo, mod_only=mod_only)
+    if plot:
+        from nns.plotting.seasonality import plot_nns_seas
+
+        plot_nns_seas(result)
+    return result
+
+
+def _nns_seas_compute(
+    variable: NDArray[np.float64],
+    *,
+    modulo: int | list[int] | NDArray[np.int64] | None = None,
+    mod_only: bool = True,
+) -> SeasonalityResult:
     values = _validate_variable(variable)
     modulo_values = None if modulo is None else _as_modulo(modulo)
     cache_key = _cache_key(values, modulo_values, mod_only)
diff --git a/tests/plotting/conftest.py b/tests/plotting/conftest.py
index 4a84b82b..7f4f2bcc 100644
--- a/tests/plotting/conftest.py
+++ b/tests/plotting/conftest.py
@@ -5,3 +5,7 @@
 import matplotlib
 
 matplotlib.use("Agg")
+
+# Plotting tests open many short-lived figures (closed per-test); don't warn
+# about the open-figure count if this session shares a process with others.
+matplotlib.rcParams["figure.max_open_warning"] = 0
diff --git a/tests/plotting/test_compute_plot_flag.py b/tests/plotting/test_compute_plot_flag.py
new file mode 100644
index 00000000..bdea5d88
--- /dev/null
+++ b/tests/plotting/test_compute_plot_flag.py
@@ -0,0 +1,120 @@
+"""`plot=True` on the compute functions now renders a figure as a side effect.
+
+This validates the wiring requested over the value-only default: each function
+still returns the same value it returns with ``plot=False`` (the contract the
+parity suite depends on), but ``plot=True`` additionally *creates* a Matplotlib
+figure via the ``nns.plotting`` layer.
+"""
+
+from __future__ import annotations
+
+import matplotlib
+
+matplotlib.use("Agg")
+
+from typing import Any
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pytest
+
+import nns
+
+
+@pytest.fixture(autouse=True)
+def _close_figs() -> object:
+    plt.close("all")
+    yield
+    plt.close("all")
+
+
+def _fig_count() -> int:
+    return len(plt.get_fignums())
+
+
+def test_nns_reg_plot_true_creates_figure() -> None:
+    rng = np.random.default_rng(0)
+    x = np.sort(rng.normal(size=40))
+    y = 2.0 * x + rng.normal(scale=0.3, size=40)
+    base = nns.nns_reg(x, y, confidence_interval=0.95)
+    plt.close("all")
+    drawn = nns.nns_reg(x, y, confidence_interval=0.95, plot=True)
+    assert _fig_count() > 0
+    assert base["R2"] == drawn["R2"]
+
+
+def test_nns_reg_residual_plot_creates_figure() -> None:
+    rng = np.random.default_rng(1)
+    x = np.sort(rng.normal(size=40))
+    y = 2.0 * x + rng.normal(scale=0.3, size=40)
+    nns.nns_reg(x, y, residual_plot=True)
+    assert _fig_count() > 0
+
+
+def test_nns_cdf_plot_true_creates_figure() -> None:
+    rng = np.random.default_rng(2)
+    v = np.cumsum(rng.normal(size=60)) + 20.0
+    base = nns.nns_cdf(v, target=20.0)
+    plt.close("all")
+    drawn = nns.nns_cdf(v, target=20.0, plot=True)
+    assert _fig_count() > 0
+    assert np.allclose(np.asarray(base["target.value"]), np.asarray(drawn["target.value"]))
+
+
+def test_nns_arma_plot_true_creates_figure() -> None:
+    rng = np.random.default_rng(3)
+    v = np.cumsum(rng.normal(size=60)) + 20.0
+    drawn = nns.nns_arma(v, h=6, pred_int=0.95, seasonal_factor=False, plot=True)
+    assert _fig_count() > 0
+    assert "Estimates" in drawn
+
+
+def test_nns_arma_optim_plot_true_creates_figure() -> None:
+    rng = np.random.default_rng(4)
+    v = np.cumsum(rng.normal(size=60)) + 20.0
+    drawn = nns.nns_arma_optim(
+        v[:40], h=6, seasonal_factor=[1], pred_int=0.95, print_trace=False, plot=True
+    )
+    assert _fig_count() > 0
+    assert "results" in drawn
+
+
+def test_nns_seas_plot_true_creates_figure() -> None:
+    rng = np.random.default_rng(5)
+    v = np.cumsum(rng.normal(size=60)) + 20.0
+    base = nns.nns_seas(v)
+    plt.close("all")
+    drawn = nns.nns_seas(v, plot=True)
+    assert _fig_count() > 0
+    assert base["best.period"] == drawn["best.period"]
+
+
+def test_nns_m_reg_plot_true_creates_figure() -> None:
+    import matplotlib.colors as mcolors
+
+    rng = np.random.default_rng(6)
+    x = np.sort(rng.normal(size=40))
+    y = 2.0 * x + rng.normal(scale=0.3, size=40)
+    features = np.column_stack([x, x**2])
+    nns.nns_m_reg(features, y, plot=True)
+    assert _fig_count() > 0
+    ax: Any = plt.gca()
+    # R's M.reg residual plot: actual y is steelblue, fitted y.hat is a red line.
+    edge_hexes = {
+        mcolors.to_hex(row)
+        for coll in ax.collections
+        for row in coll.get_edgecolor()
+        if len(row)
+    }
+    line_hexes = {mcolors.to_hex(line.get_color()) for line in ax.get_lines()}
+    assert "#4682b4" in edge_hexes
+    assert "#ff0000" in line_hexes
+
+
+def test_plot_false_creates_no_figure() -> None:
+    rng = np.random.default_rng(7)
+    x = np.sort(rng.normal(size=40))
+    y = 2.0 * x + rng.normal(scale=0.3, size=40)
+    plt.close("all")
+    nns.nns_reg(x, y)
+    assert _fig_count() == 0