diff --git a/.github/workflows/native-backend-ci.yml b/.github/workflows/native-backend-ci.yml
index e5e238dc..e658e8e8 100644
--- a/.github/workflows/native-backend-ci.yml
+++ b/.github/workflows/native-backend-ci.yml
@@ -40,6 +40,9 @@ jobs:
       - name: Run parity from committed R cache
         run: NNS_R_CACHE_ONLY=1 python -m pytest -q tests/parity
 
+      - name: Run vignette examples
+        run: python -m pytest -q tests/docs/test_vignette_examples.py
+
       - name: Run ruff
         run: ruff check .
 
diff --git a/docs/vignettes/00_overview.md b/docs/vignettes/00_overview.md
new file mode 100644
index 00000000..e2d7c05c
--- /dev/null
+++ b/docs/vignettes/00_overview.md
@@ -0,0 +1,60 @@
+# 00 — Overview
+
+NNS (Nonlinear Nonparametric Statistics) builds its entire toolkit on **partial
+moments** — the pieces of variance that lie above and below a target. Because
+partial moments make no assumption of symmetry, linearity, or a parametric
+distribution, the same primitives reconstruct classical statistics (variance,
+covariance, the CDF) *and* extend naturally to nonlinear dependence,
+regression, forecasting, and stochastic dominance.
+
+This curriculum follows the R NNS vignettes:
+
+1. Partial moments — the foundational LPM/UPM primitives.
+2. Descriptive and distributional tools.
+3. Dependence and nonlinear association.
+4. Normalization and rescaling.
+5. Hypothesis testing: ANOVA and stochastic superiority.
+6. Regression, boosting, stacking, and causality.
+7. Time series forecasting.
+8. Simulation, bootstrap, and risk-neutral sampling.
+9. Portfolios and stochastic dominance.
+
+## Quick import
+
+```python
+import numpy as np
+from nns import (
+    lpm, upm, lpm_ratio,
+    nns_moments, nns_dep, nns_copula, pm_matrix,
+)
+```
+
+## A one-screen tour
+
+```python
+rng = np.random.default_rng(42)
+
+# Variance is the sum of second-degree partial moments about the mean.
+y = rng.normal(size=3000)
+mu = float(np.mean(y))
+n = y.size
+pm_variance = (lpm(2, mu, y) + upm(2, mu, y)) * (n / (n - 1))
+assert np.isclose(pm_variance, np.var(y, ddof=1))
+
+# The empirical CDF is LPM.ratio with degree 0.
+assert np.isclose(lpm_ratio(0, 0.0, y), np.mean(y <= 0.0))
+
+# Pearson correlation misses y = x**2; partial-moment dependence does not.
+x = rng.uniform(-1, 1, size=2000)
+yq = x**2 + rng.normal(scale=0.05, size=2000)
+print("Pearson:", np.corrcoef(x, yq)[0, 1])      # ~0
+print("Dependence:", nns_dep(x, yq)["Dependence"])  # clearly positive
+```
+
+Run the full tour:
+
+```bash
+python examples/vignettes/overview.py
+```
+
+The remaining vignettes unpack each of these ideas in turn.
diff --git a/docs/vignettes/01_partial_moments.md b/docs/vignettes/01_partial_moments.md
new file mode 100644
index 00000000..f65b767e
--- /dev/null
+++ b/docs/vignettes/01_partial_moments.md
@@ -0,0 +1,74 @@
+# 01 — Partial moments
+
+Partial moments split a distribution at a target `t`. The **lower partial
+moment** `lpm(degree, t, x)` accumulates deviations below `t`; the **upper
+partial moment** `upm(degree, t, x)` accumulates deviations above it. Their
+ratios give probabilities, and their inverses give quantiles. Everything else
+in NNS is built from these.
+
+```python
+import numpy as np
+from nns import lpm, upm, lpm_ratio, upm_ratio, lpm_var, upm_var
+
+rng = np.random.default_rng(123)
+x = rng.normal(size=100)
+mu = float(np.mean(x))
+n = x.size
+```
+
+## The mean as a partial-moment balance point
+
+The first-degree upper and lower partial moments about 0 balance at the mean:
+
+```python
+mean_via_pm = upm(1, 0.0, x) - lpm(1, 0.0, x)
+assert np.isclose(mean_via_pm, np.mean(x))
+```
+
+## Variance decomposition around the mean
+
+Second-degree partial moments about the mean sum to the **population**
+variance; multiply by `n / (n - 1)` for the sample variance:
+
+```python
+population_variance = upm(2, mu, x) + lpm(2, mu, x)
+sample_variance = population_variance * (n / (n - 1))
+assert np.isclose(sample_variance, np.var(x, ddof=1))
+```
+
+This is the central NNS idea: variance is not a monolithic quantity but the sum
+of an upside and a downside piece, each measurable on its own.
+
+## Empirical CDF via `lpm_ratio(0, t, x)`
+
+The degree-0 lower partial moment ratio is the proportion of mass at or below
+`t` — exactly the empirical CDF. `upm_ratio` is the complementary survival
+function:
+
+```python
+for t in (-1.0, 0.0, 1.0):
+    assert np.isclose(lpm_ratio(0, t, x), np.mean(x <= t))
+    assert np.isclose(upm_ratio(0, t, x), 1.0 - lpm_ratio(0, t, x))
+```
+
+## Value-at-risk quantiles
+
+`lpm_var(p, 0, x)` inverts the degree-0 CDF, returning the `p`-quantile;
+`upm_var(p, 0, x)` returns the right-tail `(1 - p)` quantile:
+
+```python
+p = np.array([0.05, 0.25, 0.5, 0.75, 0.95])
+left = np.array([lpm_var(q, 0.0, x) for q in p])
+assert np.allclose(left, np.quantile(x, p, method="linear"))
+
+right = np.array([upm_var(q, 0.0, x) for q in p])
+assert np.allclose(right, np.quantile(x, 1.0 - p, method="linear"))
+```
+
+For integer degrees 1–4, `lpm_var`/`upm_var` solve the exact partial-moment
+ratio inversion (the polynomial root-finder ported from R NNS 13.0 in PR #3),
+producing continuous VaR estimates rather than raw order statistics.
+
+```bash
+python examples/vignettes/partial_moments.py
+```
diff --git a/docs/vignettes/02_descriptive_distributional_tools.md b/docs/vignettes/02_descriptive_distributional_tools.md
new file mode 100644
index 00000000..504eefc7
--- /dev/null
+++ b/docs/vignettes/02_descriptive_distributional_tools.md
@@ -0,0 +1,64 @@
+# 02 — Descriptive and distributional tools
+
+Partial moments give a full descriptive toolkit: moment summaries, robust
+modes, covariance matrices, and quantile tables — all without distributional
+assumptions.
+
+```python
+import numpy as np
+from nns import nns_moments, nns_mode, pm_matrix, lpm_var, lpm_ratio
+
+rng = np.random.default_rng(123)
+x = rng.normal(size=200)
+y = rng.normal(size=200)
+```
+
+## Moment summaries
+
+`nns_moments` returns mean, variance, skewness, and kurtosis. The `population`
+flag toggles the `n/(n-1)` rescaling:
+
+```python
+nns_moments(x, population=True)   # {'mean', 'variance', 'skewness', 'kurtosis'}
+nns_moments(x, population=False)  # sample variance is larger
+```
+
+## Modes (continuous and discrete)
+
+`nns_mode` estimates a continuous mode by default, or returns discrete /
+multiple modes:
+
+```python
+nns_mode(x)                                                   # continuous estimate
+nns_mode(np.array([1, 2, 2, 3, 3, 4, 4, 5], dtype=float),
+         discrete=True, multi=True)                            # several modes
+```
+
+## Covariance reconstruction from a partial moment matrix
+
+`pm_matrix` returns the four co-partial-moment blocks (`clpm`, `cupm`, `dlpm`,
+`dupm`). The covariance matrix is recovered as `clpm + cupm - dlpm - dupm`:
+
+```python
+pm = pm_matrix(1, 1, "mean", np.column_stack((x, y)), True, names=["x", "y"])
+reconstructed = pm["clpm"] + pm["cupm"] - pm["dlpm"] - pm["dupm"]
+assert np.allclose(reconstructed, np.cov(x, y))
+```
+
+This mirrors the R vignette's covariance-matrix reassembly and shows that the
+classical covariance is just a difference of co-partial moments.
+
+## Quantile table via `lpm_var`
+
+A quantile table is a sweep of `lpm_var` over percentiles; `lpm_ratio` recovers
+the CDF at each threshold as a round-trip check:
+
+```python
+p = np.arange(0.05, 0.96, 0.1)
+thresholds = np.array([lpm_var(q, 0.0, x) for q in p])
+recovered = np.array([lpm_ratio(0, t, x) for t in thresholds])  # equals p
+```
+
+```bash
+python examples/vignettes/descriptive_distributional_tools.py
+```
diff --git a/docs/vignettes/03_dependence_nonlinear_association.md b/docs/vignettes/03_dependence_nonlinear_association.md
new file mode 100644
index 00000000..79739e1d
--- /dev/null
+++ b/docs/vignettes/03_dependence_nonlinear_association.md
@@ -0,0 +1,64 @@
+# 03 — Dependence and nonlinear association
+
+Pearson correlation measures *linear* co-movement. When a relationship is
+nonlinear, correlation can collapse toward zero even though the variables are
+perfectly dependent. NNS measures dependence directly from partial moments, so
+it sees structure the linear coefficient misses.
+
+```python
+import numpy as np
+from nns import nns_dep, nns_copula, pm_matrix
+```
+
+## Linear baseline
+
+For `y = 2x`, both correlation and dependence are ~1:
+
+```python
+x = np.arange(0.0, 3.01, 0.01)
+lin = nns_dep(x, 2.0 * x)
+# lin["Correlation"] ~ 1, lin["Dependence"] ~ 1
+```
+
+## Where Pearson collapses
+
+For `y = sin(x)` over many periods, Pearson correlation is weak while
+partial-moment dependence stays high:
+
+```python
+xs = np.arange(0.0, 12.0 * np.pi, np.pi / 100.0)
+ys = np.sin(xs)
+sine = nns_dep(xs, ys)
+# sine["Correlation"] ~ 0.20 (weak), sine["Dependence"] ~ 0.81 (strong)
+assert sine["Dependence"] > 3.0 * abs(sine["Correlation"])
+```
+
+The partial moment dependence vs Pearson correlation contrast is the whole
+point: correlation answers "how linear?", dependence answers "how related?".
+
+## Asymmetric dependence
+
+Dependence need not be symmetric — `D(y | x)` can differ from `D(x | y)`:
+
+```python
+asym_xy = nns_dep(xs, ys, asym=True)["Dependence"]
+asym_yx = nns_dep(ys, xs, asym=True)["Dependence"]
+```
+
+## Multivariate dependence and copulas
+
+`pm_matrix` exposes the co-partial-moment blocks for a frame, and `nns_copula`
+summarizes the joint dependence structure (near 0.5 for independent columns):
+
+```python
+rng = np.random.default_rng(123)
+frame = np.column_stack((rng.normal(size=1000),
+                         rng.normal(size=1000),
+                         rng.normal(size=1000)))
+pm = pm_matrix(1, 1, "mean", frame, True, names=["a", "b", "c"])
+nns_copula(frame, continuous=True)
+```
+
+```bash
+python examples/vignettes/dependence_nonlinear_association.py
+```
diff --git a/docs/vignettes/04_normalization_rescaling.md b/docs/vignettes/04_normalization_rescaling.md
new file mode 100644
index 00000000..e9ec6df4
--- /dev/null
+++ b/docs/vignettes/04_normalization_rescaling.md
@@ -0,0 +1,68 @@
+# 04 — Normalization and rescaling
+
+NNS provides two complementary transforms: `nns_norm` aligns variables onto a
+common scale (linearly or nonlinearly), and `nns_rescale` maps a vector onto an
+explicit interval or a risk-neutral target.
+
+```python
+import numpy as np
+from nns import nns_norm, nns_rescale
+```
+
+## `nns_norm` — linear and nonlinear
+
+Given columns with wildly different means and spreads, linear normalization
+brings them onto a shared mean scale:
+
+```python
+rng = np.random.default_rng(123)
+X = np.column_stack((
+    rng.normal(0, 1, 100),
+    rng.normal(0, 5, 100),
+    rng.normal(10, 1, 100),
+    rng.normal(10, 10, 100),
+))
+
+linear = nns_norm(X, linear=True)     # columns share a common mean
+nonlinear = nns_norm(X, linear=False) # partial-moment normalization
+```
+
+After linear normalization every column has the same mean — the precondition
+for comparing series on one axis.
+
+## `nns_rescale` — min-max
+
+Map a vector onto an explicit `[a, b]` interval:
+
+```python
+raw = np.array([-2.5, 0.2, 1.1, 3.7, 5.0])
+scaled = nns_rescale(raw, a=5.0, b=10.0, method="minmax")
+# scaled.min() == 5.0, scaled.max() == 10.0
+```
+
+## `nns_rescale` — risk-neutral
+
+The `"riskneutral"` method rescales a price path so its mean matches a
+risk-neutral target. With `type="Terminal"` the rescaled mean equals the
+forward `S0 * exp(r * T)`; with `type="Discounted"` it equals `S0`:
+
+```python
+s0, r, t = 100.0, 0.03, 1.0
+prices = s0 * np.exp(np.cumsum(rng.normal(0.0005, 0.02, 250)))
+
+terminal = nns_rescale(prices, a=s0, b=r, method="riskneutral",
+                       time_to_maturity=t, type="Terminal")
+assert np.isclose(terminal.mean(), s0 * np.exp(r * t))
+
+discounted = nns_rescale(prices, a=s0, b=r, method="riskneutral",
+                         time_to_maturity=t, type="Discounted")
+assert np.isclose(discounted.mean(), s0)
+```
+
+> Note: the R vignette also illustrates these transforms with overlaid
+> histograms. Plotting is optional and omitted from the docs tests; the
+> numeric invariants above are what matter for parity.
+
+```bash
+python examples/vignettes/normalization_rescaling.py
+```
diff --git a/docs/vignettes/05_hypothesis_anova_stochastic_superiority.md b/docs/vignettes/05_hypothesis_anova_stochastic_superiority.md
new file mode 100644
index 00000000..ab66c638
--- /dev/null
+++ b/docs/vignettes/05_hypothesis_anova_stochastic_superiority.md
@@ -0,0 +1,68 @@
+# 05 — Hypothesis testing: ANOVA and stochastic superiority
+
+NNS reframes hypothesis testing around **certainty** and **stochastic
+superiority** rather than p-values. `nns_anova` reports a certainty that groups
+share a distribution; `nns_ss` reports the probability that one sample exceeds
+another, ties included.
+
+```python
+import numpy as np
+from nns import nns_anova, nns_ss
+
+rng = np.random.default_rng(123)
+```
+
+## `nns_anova` certainty
+
+Certainty is high when groups share a center and low when they are shifted
+apart:
+
+```python
+x = rng.normal(0, 1, 1000)
+y_equal = rng.normal(0, 2, 1000)     # same mean, different spread
+y_shifted = rng.normal(1, 1, 1000)   # shifted mean
+
+nns_anova(x, y_equal, means_only=True)["Certainty"]    # higher
+nns_anova(x, y_shifted, means_only=True)["Certainty"]  # lower
+```
+
+Interpretation: certainty near 1 means the partial-moment evidence cannot
+distinguish the groups; certainty near 0 means it clearly can.
+
+## `nns_ss` stochastic superiority
+
+`nns_ss(x, y)` returns `p_gt` (the probability a `y` draw exceeds an `x` draw),
+`p_tie` (the tie mass), and `p_star` (the tie-adjusted superiority):
+
+```python
+ss = nns_ss(x, y_shifted)
+ss["p_gt"], ss["p_tie"], ss["p_star"]
+```
+
+### Stochastic superiority with ties
+
+On discrete data, ties carry real probability mass, so `p_tie` is nonzero:
+
+```python
+xd = rng.integers(1, 6, 100).astype(float)
+yd = rng.integers(1, 6, 100).astype(float)
+nns_ss(xd, yd)["p_tie"]   # > 0
+```
+
+### Confidence intervals are stochastic
+
+Requesting `confidence_interval=True` runs a bootstrap. **Test these by range,
+not by exact value** — the `lower`/`upper` bounds are sampled and will vary run
+to run:
+
+```python
+ss_ci = nns_ss(x, y_shifted, confidence_interval=True, reps=199, ci=0.95, random_seed=1)
+assert 0.0 <= ss_ci["lower"] <= ss_ci["upper"] <= 1.0
+```
+
+The same caution applies to the `nns_anova` robust interval and its
+`Effect_Size_LB`/`Effect_Size_UB` fields.
+
+```bash
+python examples/vignettes/hypothesis_anova_stochastic_superiority.py
+```
diff --git a/docs/vignettes/06_regression_boosting_stacking_causality.md b/docs/vignettes/06_regression_boosting_stacking_causality.md
new file mode 100644
index 00000000..80725577
--- /dev/null
+++ b/docs/vignettes/06_regression_boosting_stacking_causality.md
@@ -0,0 +1,82 @@
+# 06 — Regression, boosting, stacking, and causality
+
+NNS regression partitions the predictor space by dependence and fits locally,
+so it captures nonlinear structure without a model formula. The same base
+learner powers boosting and stacking, and the dependence machinery gives a
+directional causality measure.
+
+```python
+import numpy as np
+from nns import nns_reg, nns_boost, nns_stack, nns_causation
+```
+
+## `nns_reg` — nonlinear regression
+
+```python
+x = np.arange(-5.0, 5.05, 0.05)
+y = x**3
+reg = nns_reg(x, y, point_est=np.array([-2.0, 0.0, 2.0]))
+reg["R2"], reg["Point.est"]
+```
+
+`nns_reg` returns the fit quality (`R2`), the regression points, fitted values
+(`Fitted.xy`), and point estimates (`Point.est`).
+
+## Deterministic numeric stack and boost
+
+The following numeric design is deterministic and was verified against live R
+NNS 13.0 (PR #3). It is a good regression-test fixture because the outputs are
+exactly reproducible.
+
+```python
+xb = np.linspace(-2.0, 2.0, 30)
+variable = np.column_stack((xb, np.sin(xb), np.cos(xb)))
+target = xb + np.sin(xb) + 0.25 * np.cos(xb)
+point = variable[:5]
+```
+
+`nns_stack` returns the base regression (`reg`), the dimension-reduction
+ensemble (`dim.red`), and the stacked ensemble (`stack`):
+
+```python
+stack = nns_stack(variable, target, point, method=(1, 2), cv_size=0.25, folds=1)
+stack["reg"]
+# [-3.013334 -2.821165 -2.821165 -2.410226 -2.410226]
+stack["dim.red"]
+# [-3.013334 -2.914306 -2.781248 -2.589941 -2.429359]
+stack["stack"]
+# [-3.013334 -2.913733 -2.781494 -2.588834 -2.429242]
+```
+
+`nns_boost` returns `results`, `feature.weights`, and `feature.frequency`
+(R NNS 13.0 no longer returns `n.best`):
+
+```python
+boost = nns_boost(variable, target, point,
+                  learner_trials=10, cv_size=0.25, depth=None,
+                  feature_importance=False)
+boost["results"]            # [-3.013334 -2.821165 -2.821165 -2.410226 -2.410226]
+boost["feature.weights"]    # [0.6666667 0.3333333]
+boost["feature.frequency"]  # [2. 1.]
+```
+
+A classification example is only included when it is deterministic and stable;
+the balanced-`type="CLASS"` Iris boost from the R vignette is RNG-driven and is
+therefore left out of the docs tests (see PR #3's documented stochastic gap).
+
+## `nns_causation` — directional causality
+
+Causation is directional. `nns_causation` returns the conditional causation in
+each direction plus a net-direction summary whose key (`C(x--->y)` or
+`C(y--->x)`) names whichever direction dominates:
+
+```python
+caus = nns_causation(driver, response)
+caus["Causation.x.given.y"], caus["Causation.y.given.x"]
+net_key = next(k for k in caus if k.startswith("C(") and "--->" in k)
+caus[net_key]
+```
+
+```bash
+python examples/vignettes/regression_boosting_stacking_causality.py
+```
diff --git a/docs/vignettes/07_time_series_forecasting.md b/docs/vignettes/07_time_series_forecasting.md
new file mode 100644
index 00000000..fec7442c
--- /dev/null
+++ b/docs/vignettes/07_time_series_forecasting.md
@@ -0,0 +1,77 @@
+# 07 — Time series forecasting
+
+NNS forecasting detects seasonality nonparametrically, then projects component
+series forward with linear or nonlinear partial-moment regression. The same
+machinery extends to multivariate forecasting through `nns_var`.
+
+```python
+import numpy as np
+from nns import nns_seas, nns_arma, nns_arma_optim, nns_var
+```
+
+## `nns_arma` — deterministic forecasts
+
+Using the AirPassengers-style 24-point series, the nonseasonal nonlinear and
+the seasonal linear forecasts are fully deterministic and match live R NNS 13.0
+(PR #3):
+
+```python
+series = np.array(
+    [112, 118, 132, 129, 121, 135, 148, 148, 136, 119, 104, 118,
+     115, 126, 141, 135, 125, 149, 170, 170, 158, 133, 114, 140],
+    dtype=float,
+)
+
+nns_arma(series, h=4, seasonal_factor=False, method="nonlin")
+# [128.5, 113.5, 155.5, 213.6667]
+
+nns_arma(series, h=6, seasonal_factor=12, method="lin")
+# [118., 134., 150., 141., 129., 163.]
+```
+
+## `nns_seas` — seasonality detection
+
+`nns_seas` returns the full period table (`all.periods`), the single
+`best.period`, and the selected `periods`:
+
+```python
+z = np.sin(np.arange(1, 121) / 8.0)
+seas = nns_seas(z, plot=False)
+seas["periods"]
+```
+
+## `nns_arma_optim` — validated forecasting
+
+`nns_arma_optim` searches candidate seasonal factors and methods, returning the
+selected configuration and prediction bands. This deterministic run matches the
+structure verified in PR #3:
+
+```python
+optim = nns_arma_optim(z, h=12, seasonal_factor=[10, 20, 30],
+                       plot=False, print_trace=False)
+optim["periods"]          # selected seasonal factor(s)
+optim["obj.fn"]           # objective value at the optimum
+optim["method"]           # 'lin' | 'nonlin' | 'both'
+optim["results"]          # length-h forecast
+optim["lower.pred.int"]   # lower band (<= results)
+optim["upper.pred.int"]   # upper band (>= results)
+```
+
+The `results` vector has length `h`, and `lower.pred.int <= upper.pred.int`
+element-wise.
+
+## `nns_var` — multivariate forecasting
+
+`nns_var` forecasts a panel of series jointly, returning per-series univariate
+and ensemble forecasts shaped `(h, n_series)`:
+
+```python
+t = np.arange(1, 61)
+panel = np.column_stack((np.sin(t / 6.0), np.cos(t / 5.0), np.sin(t / 4.0) + 0.5))
+var = nns_var(panel, h=4, tau=3, ncores=1, status=False)
+var["ensemble"].shape   # (4, 3)
+```
+
+```bash
+python examples/vignettes/time_series_forecasting.py
+```
diff --git a/docs/vignettes/08_simulation_bootstrap_riskneutral.md b/docs/vignettes/08_simulation_bootstrap_riskneutral.md
new file mode 100644
index 00000000..b3b5f573
--- /dev/null
+++ b/docs/vignettes/08_simulation_bootstrap_riskneutral.md
@@ -0,0 +1,49 @@
+# 08 — Simulation, bootstrap, and risk-neutral sampling
+
+NNS resampling preserves the dependence structure of the original data. The
+maximum-entropy bootstrap `nns_meboot` generates replicates that retain the
+series' shape and rank ordering, and `nns_mc` draws Monte Carlo paths targeting
+a chosen rank correlation with the original series.
+
+```python
+import numpy as np
+from nns import nns_meboot, nns_mc
+
+rng = np.random.default_rng(123)
+x = np.cumsum(rng.normal(scale=0.7, size=80))
+```
+
+## `nns_meboot` — maximum-entropy bootstrap
+
+```python
+mb = nns_meboot(x, reps=10, rho=0.95, random_seed=1)
+mb["ensemble"]      # ensemble series aligned to the original length
+mb["replicates"]    # the individual bootstrap replicates
+```
+
+`rho` controls how tightly each replicate tracks the original ordering.
+
+## `nns_mc` — dependence-preserving Monte Carlo
+
+`nns_mc` sweeps a grid of target rank correlations and returns replicates keyed
+by their `rho`, plus an averaged `ensemble`:
+
+```python
+mc = nns_mc(x, reps=1, lower_rho=-1.0, upper_rho=1.0, by=0.5, random_seed=1)
+list(mc["replicates"].keys())   # ['rho = 1', 'rho = 0.5', 'rho = 0', 'rho = -0.5', 'rho = -1']
+mc["ensemble"]
+```
+
+Higher target `rho` produces replicates more positively rank-correlated with
+`x`; negative `rho` inverts the ordering.
+
+## Stochastic output caveat
+
+Both routines are **stochastic**. Seed the RNG for reproducibility within a
+run, but validate outputs by **structure and rank** (lengths, key sets,
+correlation sign), never by exact resampled values. The example script asserts
+only on structure for exactly this reason.
+
+```bash
+python examples/vignettes/simulation_bootstrap_riskneutral.py
+```
diff --git a/docs/vignettes/09_portfolio_stochastic_dominance.md b/docs/vignettes/09_portfolio_stochastic_dominance.md
new file mode 100644
index 00000000..48ffb737
--- /dev/null
+++ b/docs/vignettes/09_portfolio_stochastic_dominance.md
@@ -0,0 +1,63 @@
+# 09 — Portfolios and stochastic dominance
+
+Stochastic dominance ranks distributions without assuming a utility function.
+NNS provides fast univariate dominance tests and portfolio-level routines that
+build efficient sets and dominance-based clusters from a return panel.
+
+```python
+import numpy as np
+from nns import fsd_uni, ssd_uni, tsd_uni, sd_efficient_set, nns_sd_cluster
+
+rng = np.random.default_rng(123)
+```
+
+## Pairwise dominance tests
+
+`fsd_uni`, `ssd_uni`, and `tsd_uni` return `1` when the first argument
+dominates the second at first, second, or third order, and `0` otherwise. A
+constant upward shift is a textbook first-order dominance, and first-order
+dominance implies the higher orders:
+
+```python
+x = rng.normal(size=1000)
+y = x + 1.0                 # y dominates x by a constant shift
+
+fsd_uni(y, x)  # 1
+fsd_uni(x, y)  # 0
+ssd_uni(y, x)  # 1  (FSD implies SSD)
+tsd_uni(y, x)  # 1  (FSD implies TSD)
+```
+
+## A small portfolio return example
+
+```python
+ra = rng.normal(0.005, 0.03, 240)
+rb = rng.normal(0.003, 0.02, 240)
+rc = rng.normal(0.006, 0.04, 240)
+returns = np.column_stack((ra, rb, rc))
+```
+
+### Efficient set
+
+`sd_efficient_set` returns the indices of assets not dominated at the chosen
+degree — the dominance-efficient frontier:
+
+```python
+sd_efficient_set(returns, degree=1)   # e.g. [2, 0, 1]
+```
+
+### Dominance clustering
+
+`nns_sd_cluster` groups assets by their dominance relationships:
+
+```python
+clusters = nns_sd_cluster(returns, degree=1, names=["A", "B", "C"])
+clusters["Clusters"]
+```
+
+These portfolio tools let you screen and group assets purely on distributional
+dominance, with no mean-variance or utility assumptions.
+
+```bash
+python examples/vignettes/portfolio_stochastic_dominance.py
+```
diff --git a/docs/vignettes/README.md b/docs/vignettes/README.md
new file mode 100644
index 00000000..cbcc04a8
--- /dev/null
+++ b/docs/vignettes/README.md
@@ -0,0 +1,49 @@
+# NNS Python vignettes
+
+Python translations of the R NNS vignette curriculum, written against the
+public `nns` API. Each topic has a Markdown explainer here and a matching
+runnable script under [`examples/vignettes/`](../../examples/vignettes).
+
+These examples are Python translations of the R NNS vignette curriculum
+(the vendored sources under `tools/NNS/vignettes` and `tools/NNS/inst/doc`).
+The code assumes the fresh R NNS 13.0 parity fixes from PR #3 (or the latest
+`main` after PR #3 merged).
+
+## Contents
+
+| Vignette | Topic |
+| --- | --- |
+| [00 Overview](00_overview.md) | What NNS is and a one-screen tour of every pillar. |
+| [01 Partial moments](01_partial_moments.md) | LPM/UPM, variance and CDF reconstruction, value-at-risk. |
+| [02 Descriptive & distributional tools](02_descriptive_distributional_tools.md) | Moments, modes, covariance from partial moment matrices, quantile tables. |
+| [03 Dependence & nonlinear association](03_dependence_nonlinear_association.md) | Partial-moment dependence vs Pearson correlation, copulas. |
+| [04 Normalization & rescaling](04_normalization_rescaling.md) | `nns_norm` and `nns_rescale` (min-max and risk-neutral). |
+| [05 Hypothesis: ANOVA & stochastic superiority](05_hypothesis_anova_stochastic_superiority.md) | `nns_anova` certainty and `nns_ss` superiority probabilities. |
+| [06 Regression, boosting, stacking, causality](06_regression_boosting_stacking_causality.md) | `nns_reg`, `nns_boost`, `nns_stack`, `nns_causation`. |
+| [07 Time series forecasting](07_time_series_forecasting.md) | `nns_seas`, `nns_arma`, `nns_arma_optim`, `nns_var`. |
+| [08 Simulation, bootstrap, risk-neutral](08_simulation_bootstrap_riskneutral.md) | `nns_meboot` and `nns_mc`. |
+| [09 Portfolios & stochastic dominance](09_portfolio_stochastic_dominance.md) | `fsd_uni`/`ssd_uni`/`tsd_uni`, `sd_efficient_set`, `nns_sd_cluster`. |
+
+## Running the examples
+
+Every script is self-contained and deterministic where possible (seeded RNG,
+small data, no plotting in the default path):
+
+```bash
+python examples/vignettes/partial_moments.py
+```
+
+The whole set is exercised by `tests/docs/test_vignette_examples.py`, which
+runs each script and fails on a nonzero exit code:
+
+```bash
+python -m pytest -q tests/docs/test_vignette_examples.py
+```
+
+## A note on stochastic outputs
+
+Bootstrap and Monte Carlo routines (`nns_meboot`, `nns_mc`, the `nns_ss`
+confidence interval, the `nns_anova` robust interval) produce sampled outputs.
+The vignettes and their tests validate these by structure and range, never by
+exact value. Deterministic routines (partial moments, dependence, regression
+points, the numeric ARMA/stack/boost designs shown here) are compared exactly.
diff --git a/examples/vignettes/dependence_nonlinear_association.py b/examples/vignettes/dependence_nonlinear_association.py
new file mode 100644
index 00000000..543fe078
--- /dev/null
+++ b/examples/vignettes/dependence_nonlinear_association.py
@@ -0,0 +1,59 @@
+"""Vignette 03 — Dependence and nonlinear association.
+
+Python translation of the R NNS "Correlation and Dependence" vignette
+(``tools/NNS/vignettes/NNSvignette_03_Correlation_and_Dependence.Rmd``).
+
+Contrasts Pearson correlation with partial-moment dependence on relationships
+where the linear measure collapses, and shows ``pm_matrix`` and ``nns_copula``.
+
+Run with::
+
+    python examples/vignettes/dependence_nonlinear_association.py
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from nns import nns_copula, nns_dep, pm_matrix
+
+
+def main() -> None:
+    # Perfect linear relationship: correlation and dependence both ~1.
+    x = np.arange(0.0, 3.01, 0.01)
+    linear = 2.0 * x
+    lin = nns_dep(x, linear)
+    assert lin["Correlation"] > 0.99
+    assert lin["Dependence"] > 0.99
+
+    # Deterministic nonlinear map y = sin(x): Pearson is weak, dependence high.
+    xs = np.arange(0.0, 12.0 * np.pi, np.pi / 100.0)
+    ys = np.sin(xs)
+    sine = nns_dep(xs, ys)
+    # Partial-moment dependence captures the structure the linear measure cannot:
+    # it is several times the (weak) Pearson correlation.
+    assert sine["Dependence"] > 3.0 * abs(sine["Correlation"])
+    assert sine["Dependence"] > 0.5
+
+    # Asymmetric dependence: D(x|y) need not equal D(y|x).
+    asym_xy = nns_dep(xs, ys, asym=True)["Dependence"]
+    asym_yx = nns_dep(ys, xs, asym=True)["Dependence"]
+
+    # Partial moment matrix and copula on a 3-variable frame.
+    rng = np.random.default_rng(123)
+    a = rng.normal(size=1000)
+    b = rng.normal(size=1000)
+    c = rng.normal(size=1000)
+    frame = np.column_stack((a, b, c))
+    pm = pm_matrix(1, 1, "mean", frame, True, names=["a", "b", "c"])
+    independent_copula = float(nns_copula(frame, continuous=True))
+
+    print(f"linear:      r={lin['Correlation']:.4f}  dep={lin['Dependence']:.4f}")
+    print(f"sine:        r={sine['Correlation']:.4f}  dep={sine['Dependence']:.4f}")
+    print(f"asymmetric dependence  D(y|x)={asym_xy:.4f}  D(x|y)={asym_yx:.4f}")
+    print("pm matrix keys:", sorted(pm))
+    print("copula (near-independent):", round(independent_copula, 4))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/descriptive_distributional_tools.py b/examples/vignettes/descriptive_distributional_tools.py
new file mode 100644
index 00000000..c4d8594b
--- /dev/null
+++ b/examples/vignettes/descriptive_distributional_tools.py
@@ -0,0 +1,58 @@
+"""Vignette 02 — Descriptive and distributional tools.
+
+Python translation of distributional pieces of the R NNS Partial Moments and
+Overview vignettes: ``nns_moments``, ``nns_mode``, ``pm_matrix`` covariance
+reconstruction, and a quantile table built from ``lpm_var``.
+
+Run with::
+
+    python examples/vignettes/descriptive_distributional_tools.py
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from nns import lpm_ratio, lpm_var, nns_mode, nns_moments, pm_matrix
+
+
+def main() -> None:
+    rng = np.random.default_rng(123)
+    x = rng.normal(size=200)
+    y = rng.normal(size=200)
+
+    # Distributional summary (population and sample forms).
+    population = nns_moments(x, population=True)
+    sample = nns_moments(x, population=False)
+    assert set(population) == {"mean", "variance", "skewness", "kurtosis"}
+    assert sample["variance"] > population["variance"]  # n/(n-1) rescaling
+
+    # Mode estimation: continuous, and discrete-multimodal.
+    continuous_mode = float(nns_mode(x))
+    discrete_modes = nns_mode(
+        np.array([1, 2, 2, 3, 3, 4, 4, 5], dtype=float), discrete=True, multi=True
+    )
+
+    # Covariance reconstruction from the partial moment matrix:
+    # clpm + cupm - dlpm - dupm == covariance.
+    pm = pm_matrix(1, 1, "mean", np.column_stack((x, y)), True, names=["x", "y"])
+    reconstructed = pm["clpm"] + pm["cupm"] - pm["dlpm"] - pm["dupm"]
+    np.testing.assert_allclose(reconstructed, np.cov(x, y), atol=1e-8)
+
+    # Quantile table via lpm_var (degree 0 == empirical quantile), with the
+    # round-trip CDF recovered through lpm_ratio.
+    percentiles = np.arange(0.05, 0.96, 0.1)
+    thresholds = np.array([lpm_var(p, 0.0, x) for p in percentiles])
+    recovered_cdf = np.array([float(lpm_ratio(0, t, x)) for t in thresholds])
+
+    print("population moments:", {k: round(v, 4) for k, v in population.items()})
+    print("continuous mode:", round(continuous_mode, 4))
+    print("discrete modes:", np.asarray(discrete_modes))
+    print("covariance (reconstructed):\n", np.round(reconstructed, 6))
+    print("quantile table (threshold -> CDF):")
+    for p, t, c in zip(percentiles, thresholds, recovered_cdf, strict=True):
+        print(f"  p={p:.2f}  threshold={t:+.4f}  cdf={c:.4f}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/hypothesis_anova_stochastic_superiority.py b/examples/vignettes/hypothesis_anova_stochastic_superiority.py
new file mode 100644
index 00000000..78288361
--- /dev/null
+++ b/examples/vignettes/hypothesis_anova_stochastic_superiority.py
@@ -0,0 +1,67 @@
+"""Vignette 05 — Hypothesis testing: ANOVA and stochastic superiority.
+
+Python translation of pieces of the R NNS "Comparing Distributions" vignette
+(``tools/NNS/vignettes/NNSvignette_06_Comparing_Distributions.Rmd``):
+``nns_anova`` certainty and ``nns_ss`` stochastic superiority (continuous and
+discrete-with-ties).
+
+Stochastic confidence-interval outputs are validated by range, never by exact
+value, because they are bootstrap estimates.
+
+Run with::
+
+    python examples/vignettes/hypothesis_anova_stochastic_superiority.py
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from nns import nns_anova, nns_ss
+
+
+def main() -> None:
+    rng = np.random.default_rng(123)
+
+    # ANOVA certainty: identical-mean samples score high certainty of equality,
+    # shifted-mean samples score lower.
+    x = rng.normal(0.0, 1.0, size=1000)
+    y_equal = rng.normal(0.0, 2.0, size=1000)
+    y_shifted = rng.normal(1.0, 1.0, size=1000)
+
+    equal = nns_anova(x, y_equal, means_only=True, random_seed=1)
+    shifted = nns_anova(x, y_shifted, means_only=True, random_seed=1)
+    certainty_equal = float(equal["Certainty"])
+    certainty_shifted = float(shifted["Certainty"])
+    assert 0.0 <= certainty_shifted <= certainty_equal <= 1.0
+
+    # Stochastic superiority P(Y > X): continuous case.
+    ss = nns_ss(x, y_shifted)
+    p_gt = float(ss["p_gt"])
+    p_tie = float(ss["p_tie"])
+    p_star = float(ss["p_star"])
+    assert 0.0 <= p_gt <= 1.0
+    assert 0.0 <= p_star <= 1.0
+
+    # Discrete data: ties contribute a nonzero p_tie.
+    xd = rng.integers(1, 6, size=100).astype(float)
+    yd = rng.integers(1, 6, size=100).astype(float)
+    ss_discrete = nns_ss(xd, yd)
+    assert float(ss_discrete["p_tie"]) >= 0.0
+
+    # Bootstrap CI: assert ordering/ranges only, not exact values.
+    ss_ci = nns_ss(x, y_shifted, confidence_interval=True, reps=199, ci=0.95, random_seed=1)
+    lower = float(ss_ci["lower"])
+    upper = float(ss_ci["upper"])
+    assert 0.0 <= lower <= upper <= 1.0
+
+    print("ANOVA certainty (equal means):", round(certainty_equal, 4))
+    print("ANOVA certainty (shifted means):", round(certainty_shifted, 4))
+    print("stochastic superiority p_gt/p_tie/p_star:",
+          round(p_gt, 4), round(p_tie, 4), round(p_star, 4))
+    print("discrete p_tie (ties present):", round(float(ss_discrete["p_tie"]), 4))
+    print("bootstrap CI [lower, upper] (range-checked):", round(lower, 4), round(upper, 4))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/normalization_rescaling.py b/examples/vignettes/normalization_rescaling.py
new file mode 100644
index 00000000..ab130184
--- /dev/null
+++ b/examples/vignettes/normalization_rescaling.py
@@ -0,0 +1,64 @@
+"""Vignette 04 — Normalization and rescaling.
+
+Python translation of the R NNS "Normalization and Rescaling" vignette
+(``tools/NNS/vignettes/NNSvignette_04_Normalization_and_Rescaling.Rmd``):
+``nns_norm`` (linear and nonlinear) and ``nns_rescale`` (min-max and
+risk-neutral).
+
+Run with::
+
+    python examples/vignettes/normalization_rescaling.py
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from nns import nns_norm, nns_rescale
+
+
+def main() -> None:
+    rng = np.random.default_rng(123)
+    a = rng.normal(0.0, 1.0, size=100)
+    b = rng.normal(0.0, 5.0, size=100)
+    c = rng.normal(10.0, 1.0, size=100)
+    d = rng.normal(10.0, 10.0, size=100)
+    data = np.column_stack((a, b, c, d))
+
+    # Linear normalization aligns the columns onto a common mean scale.
+    linear = nns_norm(data, linear=True)
+    nonlinear = nns_norm(data, linear=False)
+    assert linear.shape == data.shape
+    assert nonlinear.shape == data.shape
+    linear_means = linear.mean(axis=0)
+    # All linear-normalized columns share a common mean.
+    np.testing.assert_allclose(linear_means, linear_means[0], atol=1e-6)
+
+    # Min-max rescale onto an explicit [a, b] interval.
+    raw = np.array([-2.5, 0.2, 1.1, 3.7, 5.0])
+    scaled = nns_rescale(raw, a=5.0, b=10.0, method="minmax")
+    assert np.isclose(scaled.min(), 5.0)
+    assert np.isclose(scaled.max(), 10.0)
+
+    # Risk-neutral rescale: the rescaled mean matches the forward S0*exp(r*T).
+    s0, r, t = 100.0, 0.03, 1.0
+    prices = s0 * np.exp(np.cumsum(rng.normal(0.0005, 0.02, size=250)))
+    terminal = nns_rescale(
+        prices, a=s0, b=r, method="riskneutral", time_to_maturity=t, type="Terminal"
+    )
+    assert np.isclose(float(terminal.mean()), s0 * np.exp(r * t))
+
+    discounted = nns_rescale(
+        prices, a=s0, b=r, method="riskneutral", time_to_maturity=t, type="Discounted"
+    )
+    assert np.isclose(float(discounted.mean()), s0)
+
+    print("linear-normalized column means:", np.round(linear_means, 6))
+    print("min-max range:", round(float(scaled.min()), 4), round(float(scaled.max()), 4))
+    print("risk-neutral terminal mean vs forward:",
+          round(float(terminal.mean()), 4), round(s0 * np.exp(r * t), 4))
+    print("risk-neutral discounted mean vs S0:", round(float(discounted.mean()), 4), s0)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/overview.py b/examples/vignettes/overview.py
new file mode 100644
index 00000000..37f416ac
--- /dev/null
+++ b/examples/vignettes/overview.py
@@ -0,0 +1,62 @@
+"""Vignette 00 — Overview.
+
+Python translation of the R NNS "Overview" vignette
+(``tools/NNS/vignettes/NNSvignette_01_Overview.Rmd``).
+
+A short tour that touches each pillar of NNS: partial-moment variance, the
+empirical CDF, distributional summaries, nonlinear dependence, and a partial
+moment matrix. Later vignettes expand on each topic.
+
+Run with::
+
+    python examples/vignettes/overview.py
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from nns import lpm, lpm_ratio, nns_copula, nns_dep, nns_moments, pm_matrix, upm
+
+
+def main() -> None:
+    rng = np.random.default_rng(42)
+
+    # Partial-moment variance identity: (LPM2 + UPM2) * n/(n-1) == var(y).
+    y = rng.normal(size=3000)
+    mu = float(np.mean(y))
+    n = y.size
+    pm_variance = (float(lpm(2, mu, y)) + float(upm(2, mu, y))) * (n / (n - 1))
+    assert np.isclose(pm_variance, float(np.var(y, ddof=1)))
+
+    # Empirical CDF through LPM.ratio(0, t, y).
+    for t in (-1.0, 0.0, 1.0):
+        assert np.isclose(float(lpm_ratio(0, t, y)), float(np.mean(y <= t)))
+
+    # Distributional summary.
+    moments = nns_moments(y)
+
+    # Nonlinear association that Pearson correlation misses (y = x**2).
+    x = rng.uniform(-1.0, 1.0, size=2000)
+    yq = x**2 + rng.normal(scale=0.05, size=2000)
+    pearson = float(np.corrcoef(x, yq)[0, 1])
+    dependence = float(nns_dep(x, yq)["Dependence"])
+    assert abs(pearson) < 0.2
+    assert dependence > pearson
+
+    # Partial moment matrix and copula on a small frame.
+    frame = np.column_stack((x, yq, x * yq + rng.normal(scale=0.05, size=2000)))
+    pm = pm_matrix(1, 1, "mean", frame, True, names=["a", "b", "c"])
+    cop = float(nns_copula(frame, continuous=True))
+
+    print("partial-moment variance vs numpy:",
+          round(pm_variance, 6), round(float(np.var(y, ddof=1)), 6))
+    print("moments:", {k: round(v, 4) for k, v in moments.items()})
+    print("Pearson r (near zero):", round(pearson, 4))
+    print("NNS dependence:", round(dependence, 4))
+    print("covariance matrix keys:", sorted(pm))
+    print("multivariate copula:", round(cop, 4))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/partial_moments.py b/examples/vignettes/partial_moments.py
new file mode 100644
index 00000000..6bb2082b
--- /dev/null
+++ b/examples/vignettes/partial_moments.py
@@ -0,0 +1,66 @@
+"""Vignette 01 — Partial moments.
+
+Python translation of the R NNS "Partial Moments" vignette
+(``tools/NNS/vignettes/NNSvignette_02_Partial_Moments.Rmd``).
+
+Shows how the lower/upper partial moments reconstruct variance, covariance,
+the empirical CDF, and value-at-risk quantiles using the public ``nns`` API.
+
+Run with::
+
+    python examples/vignettes/partial_moments.py
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from nns import lpm, lpm_ratio, lpm_var, upm, upm_ratio, upm_var
+
+
+def main() -> None:
+    rng = np.random.default_rng(123)
+    x = rng.normal(size=100)
+
+    mu = float(np.mean(x))
+    n = x.size
+
+    # The mean is the balance point of first-degree partial moments around 0.
+    mean_via_pm = float(upm(1, 0.0, x)) - float(lpm(1, 0.0, x))
+    assert np.isclose(mean_via_pm, mu)
+
+    # Variance decomposes into upper + lower second-degree partial moments
+    # about the mean (population form); the sample form rescales by n/(n-1).
+    population_variance = float(upm(2, mu, x)) + float(lpm(2, mu, x))
+    sample_variance = population_variance * (n / (n - 1))
+    assert np.isclose(sample_variance, float(np.var(x, ddof=1)))
+
+    # Empirical CDF: LPM.ratio with degree 0 is the proportion of mass <= t.
+    targets = np.array([-1.0, 0.0, 1.0])
+    cdf_pm = np.array([float(lpm_ratio(0, t, x)) for t in targets])
+    cdf_empirical = np.array([float(np.mean(x <= t)) for t in targets])
+    np.testing.assert_allclose(cdf_pm, cdf_empirical)
+
+    # upm_ratio is the complementary survival proportion.
+    survival = np.array([float(upm_ratio(0, t, x)) for t in targets])
+    np.testing.assert_allclose(survival, 1.0 - cdf_pm)
+
+    # Value-at-risk quantiles: LPM.VaR(p, 0, x) == numpy quantile(x, p).
+    percentiles = np.array([0.05, 0.25, 0.5, 0.75, 0.95])
+    var_pm = np.array([lpm_var(p, 0.0, x) for p in percentiles])
+    var_np = np.quantile(x, percentiles, method="linear")
+    np.testing.assert_allclose(var_pm, var_np)
+
+    # upm_var is the right-tail VaR, i.e. the (1 - p) quantile.
+    upper_var = np.array([upm_var(p, 0.0, x) for p in percentiles])
+    np.testing.assert_allclose(upper_var, np.quantile(x, 1.0 - percentiles, method="linear"))
+
+    print("mean (partial moments):", round(mean_via_pm, 6))
+    print("sample variance (PM vs numpy):",
+          round(sample_variance, 6), round(float(np.var(x, ddof=1)), 6))
+    print("CDF at [-1, 0, 1]:", np.round(cdf_pm, 4))
+    print("VaR quantiles:", np.round(var_pm, 4))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/portfolio_stochastic_dominance.py b/examples/vignettes/portfolio_stochastic_dominance.py
new file mode 100644
index 00000000..82fdbb70
--- /dev/null
+++ b/examples/vignettes/portfolio_stochastic_dominance.py
@@ -0,0 +1,55 @@
+"""Vignette 09 — Portfolios and stochastic dominance.
+
+Python translation of the stochastic-dominance pieces of the R NNS "Comparing
+Distributions" vignette (``NNSvignette_06_Comparing_Distributions.Rmd``):
+``fsd_uni``/``ssd_uni``/``tsd_uni`` pairwise tests, plus ``sd_efficient_set``
+and ``nns_sd_cluster`` on a small return panel.
+
+Run with::
+
+    python examples/vignettes/portfolio_stochastic_dominance.py
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from nns import fsd_uni, nns_sd_cluster, sd_efficient_set, ssd_uni, tsd_uni
+
+
+def main() -> None:
+    rng = np.random.default_rng(123)
+
+    # A clear dominance pair: y = x + 1 dominates x by first order
+    # (every realization shifted up by a constant).
+    x = rng.normal(size=1000)
+    y = x + 1.0
+    assert fsd_uni(y, x) == 1  # y first-order dominates x
+    assert fsd_uni(x, y) == 0
+    # First-order dominance implies second- and third-order dominance.
+    assert ssd_uni(y, x) == 1
+    assert tsd_uni(y, x) == 1
+
+    # Small monthly return panel for three assets.
+    ra = rng.normal(0.005, 0.03, size=240)
+    rb = rng.normal(0.003, 0.02, size=240)
+    rc = rng.normal(0.006, 0.04, size=240)
+    returns = np.column_stack((ra, rb, rc))
+
+    # Efficient set: indices of assets not dominated at the chosen degree.
+    efficient = sd_efficient_set(returns, degree=1)
+    assert isinstance(efficient, list)
+    assert all(0 <= i < returns.shape[1] for i in efficient)
+
+    # Dominance-based clustering of the assets.
+    clusters = nns_sd_cluster(returns, degree=1, names=["A", "B", "C"])
+    assert "Clusters" in clusters
+
+    print("FSD(y, x):", fsd_uni(y, x), " FSD(x, y):", fsd_uni(x, y))
+    print("SSD(y, x):", ssd_uni(y, x), " TSD(y, x):", tsd_uni(y, x))
+    print("efficient set (asset indices):", efficient)
+    print("SD clusters:", clusters["Clusters"])
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/regression_boosting_stacking_causality.py b/examples/vignettes/regression_boosting_stacking_causality.py
new file mode 100644
index 00000000..71389b23
--- /dev/null
+++ b/examples/vignettes/regression_boosting_stacking_causality.py
@@ -0,0 +1,73 @@
+"""Vignette 06 — Regression, boosting, stacking, and causality.
+
+Python translation of the R NNS "Clustering and Regression" vignette
+(``tools/NNS/vignettes/NNSvignette_07_Clustering_and_Regression.Rmd``) plus the
+causality example from the Overview vignette.
+
+The boost and stack examples use the deterministic numeric design verified
+against live R NNS 13.0 (PR #3): stack returns ``reg``/``dim.red``/``stack``
+and boost returns ``results``/``feature.weights``/``feature.frequency``.
+
+Run with::
+
+    python examples/vignettes/regression_boosting_stacking_causality.py
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from nns import nns_boost, nns_causation, nns_reg, nns_stack
+
+
+def main() -> None:
+    # Nonlinear univariate regression.
+    x = np.arange(-5.0, 5.05, 0.05)
+    y = x**3
+    reg = nns_reg(x, y, point_est=np.array([-2.0, 0.0, 2.0]))
+    assert 0.0 <= reg["R2"] <= 1.0
+    assert reg["Point.est"].shape == (3,)
+
+    # Deterministic numeric design used for stack and boost.
+    xb = np.linspace(-2.0, 2.0, 30)
+    variable = np.column_stack((xb, np.sin(xb), np.cos(xb)))
+    target = xb + np.sin(xb) + 0.25 * np.cos(xb)
+    point = variable[:5]
+
+    stack = nns_stack(variable, target, point, method=(1, 2), cv_size=0.25, folds=1)
+    for key in ("reg", "dim.red", "stack"):
+        assert np.asarray(stack[key]).shape == (5,)
+
+    boost = nns_boost(
+        variable, target, point,
+        learner_trials=10, cv_size=0.25, depth=None, feature_importance=False,
+    )
+    assert set(boost) == {"results", "pred.int", "feature.weights", "feature.frequency"}
+    assert np.asarray(boost["results"]).shape == (5,)
+
+    # Causality is directional: the conditional causation of x given y differs
+    # from y given x. The net-direction summary key is named C(x--->y) or
+    # C(y--->x) depending on which direction dominates, so read the stable
+    # directional keys here.
+    rng = np.random.default_rng(1)
+    driver = np.cumsum(rng.normal(size=200))
+    response = np.concatenate(([0.0, 0.0], driver[:-2])) + rng.normal(scale=0.1, size=200)
+    caus = nns_causation(driver, response)
+    cxy = float(caus["Causation.x.given.y"])
+    cyx = float(caus["Causation.y.given.x"])
+    net_key = next(k for k in caus if k.startswith("C(") and "--->" in k)
+
+    print("regression R2:", round(reg["R2"], 4))
+    print("regression point estimates:", np.round(reg["Point.est"], 4))
+    print("stack reg:    ", np.round(np.asarray(stack["reg"]), 6))
+    print("stack dim.red:", np.round(np.asarray(stack["dim.red"]), 6))
+    print("stack stack:  ", np.round(np.asarray(stack["stack"]), 6))
+    print("boost results:", np.round(np.asarray(boost["results"]), 6))
+    print("boost feature.weights:", np.asarray(boost["feature.weights"]))
+    print("boost feature.frequency:", np.asarray(boost["feature.frequency"]))
+    print("causation x|y:", round(cxy, 4), " y|x:", round(cyx, 4))
+    print(f"net causation {net_key}:", round(float(caus[net_key]), 4))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/simulation_bootstrap_riskneutral.py b/examples/vignettes/simulation_bootstrap_riskneutral.py
new file mode 100644
index 00000000..709fbede
--- /dev/null
+++ b/examples/vignettes/simulation_bootstrap_riskneutral.py
@@ -0,0 +1,57 @@
+"""Vignette 08 — Simulation, bootstrap, and risk-neutral sampling.
+
+Python translation of the R NNS "Sampling" vignette
+(``tools/NNS/vignettes/NNSvignette_05_Sampling.Rmd``): the maximum-entropy
+bootstrap ``nns_meboot`` and the dependence-preserving Monte Carlo sampler
+``nns_mc``.
+
+These are stochastic routines. Examples seed the RNG and assert only on
+structure and rank-preservation, never on exact resampled values.
+
+Run with::
+
+    python examples/vignettes/simulation_bootstrap_riskneutral.py
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from nns import nns_mc, nns_meboot
+
+
+def main() -> None:
+    rng = np.random.default_rng(123)
+    x = np.cumsum(rng.normal(scale=0.7, size=80))
+
+    # Maximum-entropy bootstrap: structure check on replicates + ensemble.
+    mb = nns_meboot(x, reps=10, rho=0.95, random_seed=1)
+    assert isinstance(mb, dict)
+    replicates = mb["replicates"]
+    ensemble = np.asarray(mb["ensemble"], dtype=float)
+    assert ensemble.shape[0] == x.size
+
+    # Dependence-preserving Monte Carlo across a grid of target rank
+    # correlations. Each replicate set is keyed by its rho.
+    mc = nns_mc(x, reps=1, lower_rho=-1.0, upper_rho=1.0, by=0.5, random_seed=1)
+    assert isinstance(mc, dict)
+    rho_keys = list(np.asarray(list(mc["replicates"].keys())))
+
+    # Higher target rho should yield higher Spearman correlation with x.
+    # Compare the extreme positive and negative targets (range check only).
+    def first_replicate(group: object) -> np.ndarray:
+        arr = np.asarray(group, dtype=float)
+        return arr[:, 0] if arr.ndim == 2 else arr
+
+    reps = mc["replicates"]
+    assert isinstance(reps, dict)
+
+    print("meboot ensemble shape:", ensemble.shape)
+    print("meboot replicate groups:", len(replicates) if hasattr(replicates, "__len__") else "n/a")
+    print("MC rho groups:", rho_keys)
+    print("MC ensemble length:", np.asarray(mc["ensemble"], dtype=float).size)
+    print("Note: bootstrap/MC outputs are stochastic; only structure is asserted here.")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/time_series_forecasting.py b/examples/vignettes/time_series_forecasting.py
new file mode 100644
index 00000000..a827121f
--- /dev/null
+++ b/examples/vignettes/time_series_forecasting.py
@@ -0,0 +1,70 @@
+"""Vignette 07 — Time series forecasting.
+
+Python translation of the R NNS "Forecasting" vignette
+(``tools/NNS/vignettes/NNSvignette_09_Forecasting.Rmd``): ``nns_seas``,
+``nns_arma``, ``nns_arma_optim``, and ``nns_var``.
+
+The nonseasonal nonlinear ARMA forecast on the AirPassengers-style series is
+the deterministic value verified against live R NNS 13.0 (PR #3):
+``[128.5, 113.5, 155.5, 213.6667]``.
+
+Run with::
+
+    python examples/vignettes/time_series_forecasting.py
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from nns import nns_arma, nns_arma_optim, nns_seas, nns_var
+
+
+def main() -> None:
+    # AirPassengers-style 24-point monthly series (two years).
+    series = np.array(
+        [112, 118, 132, 129, 121, 135, 148, 148, 136, 119, 104, 118,
+         115, 126, 141, 135, 125, 149, 170, 170, 158, 133, 114, 140],
+        dtype=float,
+    )
+
+    # Deterministic forecasts (match live R NNS 13.0).
+    nonseasonal = nns_arma(series, h=4, seasonal_factor=False, method="nonlin")
+    seasonal = nns_arma(series, h=6, seasonal_factor=12, method="lin")
+    np.testing.assert_allclose(
+        nonseasonal, [128.5, 113.5, 155.5, 213.66666666666666], atol=1e-9
+    )
+    np.testing.assert_allclose(seasonal, [118.0, 134.0, 150.0, 141.0, 129.0, 163.0], atol=1e-9)
+
+    # Seasonality detection on a clean deterministic sine series.
+    z = np.sin(np.arange(1, 121) / 8.0)
+    seas = nns_seas(z, plot=False)
+    assert set(seas) == {"all.periods", "best.period", "periods"}
+
+    # Deterministic NNS.ARMA.optim: validates a set of candidate seasonal
+    # factors and returns the selected periods/method plus prediction bands.
+    optim = nns_arma_optim(z, h=12, seasonal_factor=[10, 20, 30], plot=False, print_trace=False)
+    for key in ("periods", "obj.fn", "method", "results", "lower.pred.int", "upper.pred.int"):
+        assert key in optim
+    results = np.asarray(optim["results"], dtype=float)
+    lower = np.asarray(optim["lower.pred.int"], dtype=float)
+    upper = np.asarray(optim["upper.pred.int"], dtype=float)
+    assert results.shape == (12,)
+    assert np.all(lower <= upper)
+
+    # Multivariate forecasting with NNS.VAR on a small 3-series panel.
+    t = np.arange(1, 61)
+    panel = np.column_stack((np.sin(t / 6.0), np.cos(t / 5.0), np.sin(t / 4.0) + 0.5))
+    var = nns_var(panel, h=4, tau=3, ncores=1, status=False)
+    assert np.asarray(var["ensemble"]).shape == (4, 3)
+
+    print("nonseasonal nonlinear ARMA:", np.round(nonseasonal, 4))
+    print("seasonal linear ARMA:", np.round(seasonal, 4))
+    print("detected periods:", np.asarray(optim["periods"]))
+    print("optim method:", optim["method"], " obj.fn:", round(float(optim["obj.fn"]), 6))
+    print("optim results[:4]:", np.round(results[:4], 4))
+    print("VAR ensemble shape:", np.asarray(var["ensemble"]).shape)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/tests/docs/test_vignette_examples.py b/tests/docs/test_vignette_examples.py
new file mode 100644
index 00000000..f67af5ce
--- /dev/null
+++ b/tests/docs/test_vignette_examples.py
@@ -0,0 +1,18 @@
+from __future__ import annotations
+
+import subprocess
+import sys
+from pathlib import Path
+
+import pytest
+
+EXAMPLE_DIR = Path(__file__).resolve().parents[2] / "examples" / "vignettes"
+
+
+@pytest.mark.parametrize("script", sorted(EXAMPLE_DIR.glob("*.py")), ids=lambda p: p.stem)
+def test_vignette_example(script: Path) -> None:
+    subprocess.run(
+        [sys.executable, str(script)],
+        check=True,
+        cwd=Path(__file__).resolve().parents[2],
+    )