diff --git a/.github/workflows/ci-python.yml b/.github/workflows/ci-python.yml
index dc87723..20daf3d 100644
--- a/.github/workflows/ci-python.yml
+++ b/.github/workflows/ci-python.yml
@@ -125,3 +125,43 @@ jobs:
         with:
           name: benchmark-vs-talib
           path: benchmark_vs_talib.json
+
+  perf-smoke:
+    name: Performance smoke and contracts
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v6
+
+      - name: Set up Python 3.12
+        uses: actions/setup-python@v6
+        with:
+          python-version: "3.12"
+
+      - name: Install maturin and perf dependencies
+        run: |
+          pip install maturin numpy pytest
+
+      - name: Build and install ferro_ta
+        run: |
+          maturin build --release --out dist
+          pip install dist/*.whl
+
+      - name: Generate reproducible perf artifacts
+        run: |
+          python benchmarks/run_perf_contract.py \
+            --output-dir perf-contract \
+            --skip-talib \
+            --batch-samples 20000 \
+            --batch-series 32 \
+            --streaming-bars 20000 \
+            --price-bars 20000 \
+            --iv-bars 50000
+
+      - name: Enforce hotspot regression policy
+        run: python benchmarks/check_hotspot_regression.py --input perf-contract/runtime_hotspots.json
+
+      - name: Upload perf artifacts
+        uses: actions/upload-artifact@v7
+        with:
+          name: perf-contract
+          path: perf-contract/
diff --git a/.github/workflows/ci-wasm.yml b/.github/workflows/ci-wasm.yml
index 91d493f..d1636ef 100644
--- a/.github/workflows/ci-wasm.yml
+++ b/.github/workflows/ci-wasm.yml
@@ -13,6 +13,11 @@ jobs:
     steps:
       - uses: actions/checkout@v6
 
+      - name: Set up Node.js
+        uses: actions/setup-node@v4
+        with:
+          node-version: "20"
+
       - name: Install Rust (stable)
         uses: dtolnay/rust-toolchain@v1
         with:
@@ -30,8 +35,18 @@ jobs:
         working-directory: wasm
         run: wasm-pack build --target nodejs --out-dir pkg
 
+      - name: Benchmark WASM package
+        working-directory: wasm
+        run: node bench.js --json ../wasm_benchmark.json
+
       - name: Upload WASM package artifact
         uses: actions/upload-artifact@v7
         with:
           name: wasm-pkg
           path: wasm/pkg/
+
+      - name: Upload WASM benchmark artifact
+        uses: actions/upload-artifact@v7
+        with:
+          name: wasm-benchmark
+          path: wasm_benchmark.json
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c819211..597a0c0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,25 @@ and the project uses [Semantic Versioning](https://semver.org/).
 
 ## [Unreleased]
 
+## [1.0.2] — 2026-03-24
+
+### Performance
+
+- Optimized rolling statistical kernels (`CORREL`, `BETA`, `LINEARREG*`, `TSF`)
+  with incremental window math and matching warmup semantics.
+- Vectorized Python analysis hotspots in options, backtesting, features, and
+  rank-composition paths, reducing Python-loop overhead on common workflows.
+- Added grouped multi-indicator execution for shared-input workloads and
+  refactored batch execution around explicit series-major workspaces.
+
+### Added
+
+- Reproducible perf-contract artifacts for single-series, batch, streaming,
+  SIMD, TA-Lib comparison, and WASM benchmark runs.
+- Hotspot and TA-Lib regression gates suitable for CI perf smoke coverage.
+- Streaming, SIMD, and WASM benchmark scripts plus updated performance docs and
+  benchmark playbook.
+
 ## [1.0.1] — 2026-03-24
 
 ### Added
@@ -274,6 +293,7 @@ and the project uses [Semantic Versioning](https://semver.org/).
 
 ---
 
-[Unreleased]: https://github.com/pratikbhadane24/ferro-ta/compare/v1.0.1...HEAD
+[Unreleased]: https://github.com/pratikbhadane24/ferro-ta/compare/v1.0.2...HEAD
+[1.0.2]: https://github.com/pratikbhadane24/ferro-ta/compare/v1.0.1...v1.0.2
 [1.0.1]: https://github.com/pratikbhadane24/ferro-ta/compare/v1.0.0...v1.0.1
 [1.0.0]: https://github.com/pratikbhadane24/ferro-ta/releases/tag/v1.0.0
diff --git a/Cargo.lock b/Cargo.lock
index ff595fa..5340f3d 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -207,7 +207,7 @@ checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719"
 
 [[package]]
 name = "ferro_ta"
-version = "1.0.1"
+version = "1.0.2"
 dependencies = [
  "criterion",
  "ferro_ta_core",
@@ -222,7 +222,7 @@ dependencies = [
 
 [[package]]
 name = "ferro_ta_core"
-version = "1.0.1"
+version = "1.0.2"
 dependencies = [
  "criterion",
  "wide",
diff --git a/Cargo.toml b/Cargo.toml
index e394178..83e017d 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -5,7 +5,7 @@ resolver = "2"
 
 [package]
 name = "ferro_ta"
-version = "1.0.1"
+version = "1.0.2"
 edition = "2021"
 license = "MIT"
 publish = false
@@ -23,7 +23,7 @@ ndarray = "0.16"
 rayon = "1.10"
 log = "0.4"
 pyo3-log = "0.12"
-ferro_ta_core = { path = "crates/ferro_ta_core", version = "1.0.1" }
+ferro_ta_core = { path = "crates/ferro_ta_core", version = "1.0.2" }
 
 [dev-dependencies]
 criterion = { version = "0.8", features = ["html_reports"] }
diff --git a/benchmarks/README.md b/benchmarks/README.md
index b55a55c..c534d7e 100644
--- a/benchmarks/README.md
+++ b/benchmarks/README.md
@@ -40,6 +40,30 @@ from benchmarks.data_generator import SMALL, MEDIUM, LARGE
 - **Machine info:** Stored in `benchmarks/results.json` (`machine_info`, `commit_info`) for reproducibility.
 - **Libraries:** Only libraries present in the environment are benchmarked; missing ones are skipped.
 
+## Reproducible Perf Artifacts
+
+Use the perf-contract runner when you want a compact set of machine-readable
+artifacts for single-series latency, batch throughput, streaming throughput,
+and hotspot attribution in one directory:
+
+```bash
+uv run python benchmarks/run_perf_contract.py --output-dir benchmarks/artifacts/latest --skip-talib
+```
+
+That command writes:
+
+- `indicator_latency.json` — canonical-fixture timings for the benchmark suite indicators
+- `batch.json` — 2-D batch throughput plus grouped multi-indicator timings
+- `streaming.json` — streaming update throughput vs batch baselines
+- `runtime_hotspots.json` — ranked hotspot report with reference speedups
+- `manifest.json` — runtime/git metadata plus hashes for the generated artifacts
+
+For CI or local guardrails, validate the hotspot report with:
+
+```bash
+uv run python benchmarks/check_hotspot_regression.py --input benchmarks/artifacts/latest/runtime_hotspots.json
+```
+
 ---
 
 ## Speed comparison (100k bars, median µs — lower is better)
@@ -141,10 +165,34 @@ uv run python benchmarks/bench_vs_talib.py --sizes 10000 100000 --json benchmark
 
 # Optional regression check used in CI
 uv run python benchmarks/check_vs_talib_regression.py --input benchmark_vs_talib.json
+
+# Batch throughput + grouped multi-indicator calls
+uv run python benchmarks/bench_batch.py --samples 100000 --series 100 --json batch_benchmark.json
+
+# Streaming update throughput vs batch baselines
+uv run python benchmarks/bench_streaming.py --bars 100000 --json streaming_benchmark.json
+
+# Ranked hotspot attribution against bundled reference implementations
+uv run python benchmarks/profile_runtime_hotspots.py --json runtime_hotspots.json
+
+# Portable vs SIMD-enabled build comparison
+uv run python benchmarks/bench_simd.py --json simd_benchmark.json
+
+# One-shot perf artifact bundle
+uv run python benchmarks/run_perf_contract.py --output-dir benchmarks/artifacts/latest
 ```
 
 Without `uv`: use `pytest` and `python` from the same environment where `ferro_ta` and optional libs (e.g. `talib`, `pandas_ta`, `ta`, `tulipy`, `finta`) are installed.
 
+### WASM
+
+From the `wasm/` directory:
+
+```bash
+wasm-pack build --target nodejs --out-dir pkg
+node bench.js --json ../wasm_benchmark.json
+```
+
 ---
 
 ## Indicator coverage
diff --git a/benchmarks/artifacts/latest/batch.json b/benchmarks/artifacts/latest/batch.json
new file mode 100644
index 0000000..89a46a4
--- /dev/null
+++ b/benchmarks/artifacts/latest/batch.json
@@ -0,0 +1,71 @@
+{
+  "metadata": {
+    "suite": "batch",
+    "runtime": {
+      "generated_at_utc": "2026-03-23T20:25:58.345834+00:00",
+      "python_version": "3.13.5",
+      "platform": "macOS-26.3.1-arm64-arm-64bit-Mach-O",
+      "machine": "arm64",
+      "processor": "arm"
+    },
+    "git": {
+      "commit": "9011250f992119170242cf17a67834c67b91bcdb",
+      "dirty": true,
+      "branch": "feat/performace-1.0.2"
+    },
+    "dataset": {
+      "n_samples": 100000,
+      "n_series": 100,
+      "total_bars": 10000000,
+      "seed": 42
+    }
+  },
+  "results": [
+    {
+      "indicator": "SMA",
+      "parallel_ms": 37.86,
+      "sequential_ms": 43.5625,
+      "loop_ms": 17.8136,
+      "parallel_speedup_vs_loop": 0.4705,
+      "sequential_speedup_vs_loop": 0.4089
+    },
+    {
+      "indicator": "RSI",
+      "parallel_ms": 40.9229,
+      "sequential_ms": 79.5345,
+      "loop_ms": 53.3368,
+      "parallel_speedup_vs_loop": 1.3033,
+      "sequential_speedup_vs_loop": 0.6706
+    },
+    {
+      "indicator": "ATR",
+      "parallel_ms": 91.76,
+      "sequential_ms": 130.1404,
+      "loop_ms": 99.5885,
+      "parallel_speedup_vs_loop": 1.0853,
+      "sequential_speedup_vs_loop": 0.7652
+    },
+    {
+      "indicator": "ADX",
+      "parallel_ms": 100.1362,
+      "sequential_ms": 149.3412,
+      "loop_ms": 125.3319,
+      "parallel_speedup_vs_loop": 1.2516,
+      "sequential_speedup_vs_loop": 0.8392
+    }
+  ],
+  "grouped_results": [
+    {
+      "case": "close_bundle_3",
+      "grouped_ms": 0.652,
+      "separate_ms": 0.9124,
+      "speedup_vs_separate": 1.3994
+    },
+    {
+      "case": "hlc_bundle_3",
+      "grouped_ms": 1.4784,
+      "separate_ms": 3.3724,
+      "speedup_vs_separate": 2.2811
+    }
+  ]
+}
\ No newline at end of file
diff --git a/benchmarks/artifacts/latest/benchmark_vs_talib.json b/benchmarks/artifacts/latest/benchmark_vs_talib.json
new file mode 100644
index 0000000..116cf5e
--- /dev/null
+++ b/benchmarks/artifacts/latest/benchmark_vs_talib.json
@@ -0,0 +1,262 @@
+{
+  "schema_version": 1,
+  "command": "python benchmarks/bench_vs_talib.py",
+  "n_warmup": 1,
+  "n_runs": 7,
+  "sizes": [
+    10000,
+    100000
+  ],
+  "talib_available": true,
+  "runtime": {
+    "generated_at_utc": "2026-03-23T20:26:40.738132+00:00",
+    "python_version": "3.13.5",
+    "platform": "macOS-26.3.1-arm64-arm-64bit-Mach-O",
+    "machine": "arm64"
+  },
+  "git": {
+    "commit": "9011250f992119170242cf17a67834c67b91bcdb",
+    "dirty": true
+  },
+  "summary": {
+    "total_rows": 24,
+    "by_size": [
+      {
+        "size": 10000,
+        "rows": 12,
+        "wins": 8,
+        "win_rate": 0.6666666666666666,
+        "median_speedup": 1.0546,
+        "min_speedup": 0.7174,
+        "max_speedup": 2.0135
+      },
+      {
+        "size": 100000,
+        "rows": 12,
+        "wins": 7,
+        "win_rate": 0.5833333333333334,
+        "median_speedup": 1.0896,
+        "min_speedup": 0.4965,
+        "max_speedup": 3.6029
+      }
+    ]
+  },
+  "results": [
+    {
+      "indicator": "SMA",
+      "size": 10000,
+      "ferro_ta_ms": 0.0093,
+      "talib_ms": 0.0167,
+      "speedup": 1.7937,
+      "ferro_ta_m_bars_s": 1076.19,
+      "talib_m_bars_s": 599.99
+    },
+    {
+      "indicator": "SMA",
+      "size": 100000,
+      "ferro_ta_ms": 0.0765,
+      "talib_ms": 0.1346,
+      "speedup": 1.7589,
+      "ferro_ta_m_bars_s": 1306.49,
+      "talib_m_bars_s": 742.8
+    },
+    {
+      "indicator": "EMA",
+      "size": 10000,
+      "ferro_ta_ms": 0.02,
+      "talib_ms": 0.0214,
+      "speedup": 1.0687,
+      "ferro_ta_m_bars_s": 500.0,
+      "talib_m_bars_s": 467.84
+    },
+    {
+      "indicator": "EMA",
+      "size": 100000,
+      "ferro_ta_ms": 0.1998,
+      "talib_ms": 0.1883,
+      "speedup": 0.9425,
+      "ferro_ta_m_bars_s": 500.42,
+      "talib_m_bars_s": 530.97
+    },
+    {
+      "indicator": "RSI",
+      "size": 10000,
+      "ferro_ta_ms": 0.0475,
+      "talib_ms": 0.048,
+      "speedup": 1.0088,
+      "ferro_ta_m_bars_s": 210.34,
+      "talib_m_bars_s": 208.52
+    },
+    {
+      "indicator": "RSI",
+      "size": 100000,
+      "ferro_ta_ms": 0.4804,
+      "talib_ms": 0.4635,
+      "speedup": 0.9647,
+      "ferro_ta_m_bars_s": 208.15,
+      "talib_m_bars_s": 215.77
+    },
+    {
+      "indicator": "BBANDS",
+      "size": 10000,
+      "ferro_ta_ms": 0.0215,
+      "talib_ms": 0.0433,
+      "speedup": 2.0135,
+      "ferro_ta_m_bars_s": 465.12,
+      "talib_m_bars_s": 230.99
+    },
+    {
+      "indicator": "BBANDS",
+      "size": 100000,
+      "ferro_ta_ms": 0.1686,
+      "talib_ms": 0.4275,
+      "speedup": 2.535,
+      "ferro_ta_m_bars_s": 593.03,
+      "talib_m_bars_s": 233.94
+    },
+    {
+      "indicator": "MACD",
+      "size": 10000,
+      "ferro_ta_ms": 0.0495,
+      "talib_ms": 0.065,
+      "speedup": 1.3134,
+      "ferro_ta_m_bars_s": 202.19,
+      "talib_m_bars_s": 153.95
+    },
+    {
+      "indicator": "MACD",
+      "size": 100000,
+      "ferro_ta_ms": 0.4459,
+      "talib_ms": 0.6334,
+      "speedup": 1.4205,
+      "ferro_ta_m_bars_s": 224.28,
+      "talib_m_bars_s": 157.88
+    },
+    {
+      "indicator": "ATR",
+      "size": 10000,
+      "ferro_ta_ms": 0.0483,
+      "talib_ms": 0.0502,
+      "speedup": 1.0405,
+      "ferro_ta_m_bars_s": 207.07,
+      "talib_m_bars_s": 199.0
+    },
+    {
+      "indicator": "ATR",
+      "size": 100000,
+      "ferro_ta_ms": 0.4705,
+      "talib_ms": 0.4788,
+      "speedup": 1.0174,
+      "ferro_ta_m_bars_s": 212.52,
+      "talib_m_bars_s": 208.88
+    },
+    {
+      "indicator": "STOCH",
+      "size": 10000,
+      "ferro_ta_ms": 0.0915,
+      "talib_ms": 0.066,
+      "speedup": 0.721,
+      "ferro_ta_m_bars_s": 109.24,
+      "talib_m_bars_s": 151.52
+    },
+    {
+      "indicator": "STOCH",
+      "size": 100000,
+      "ferro_ta_ms": 1.5285,
+      "talib_ms": 0.7589,
+      "speedup": 0.4965,
+      "ferro_ta_m_bars_s": 65.42,
+      "talib_m_bars_s": 131.77
+    },
+    {
+      "indicator": "ADX",
+      "size": 10000,
+      "ferro_ta_ms": 0.0692,
+      "talib_ms": 0.0522,
+      "speedup": 0.7538,
+      "ferro_ta_m_bars_s": 144.49,
+      "talib_m_bars_s": 191.7
+    },
+    {
+      "indicator": "ADX",
+      "size": 100000,
+      "ferro_ta_ms": 0.6209,
+      "talib_ms": 0.5731,
+      "speedup": 0.923,
+      "ferro_ta_m_bars_s": 161.05,
+      "talib_m_bars_s": 174.49
+    },
+    {
+      "indicator": "CCI",
+      "size": 10000,
+      "ferro_ta_ms": 0.0731,
+      "talib_ms": 0.0829,
+      "speedup": 1.1333,
+      "ferro_ta_m_bars_s": 136.75,
+      "talib_m_bars_s": 120.66
+    },
+    {
+      "indicator": "CCI",
+      "size": 100000,
+      "ferro_ta_ms": 0.7028,
+      "talib_ms": 0.8164,
+      "speedup": 1.1617,
+      "ferro_ta_m_bars_s": 142.3,
+      "talib_m_bars_s": 122.49
+    },
+    {
+      "indicator": "OBV",
+      "size": 10000,
+      "ferro_ta_ms": 0.0156,
+      "talib_ms": 0.0112,
+      "speedup": 0.7174,
+      "ferro_ta_m_bars_s": 640.0,
+      "talib_m_bars_s": 892.14
+    },
+    {
+      "indicator": "OBV",
+      "size": 100000,
+      "ferro_ta_ms": 0.2953,
+      "talib_ms": 0.2808,
+      "speedup": 0.9509,
+      "ferro_ta_m_bars_s": 338.6,
+      "talib_m_bars_s": 356.08
+    },
+    {
+      "indicator": "MFI",
+      "size": 10000,
+      "ferro_ta_ms": 0.0239,
+      "talib_ms": 0.0203,
+      "speedup": 0.8516,
+      "ferro_ta_m_bars_s": 418.85,
+      "talib_m_bars_s": 491.81
+    },
+    {
+      "indicator": "MFI",
+      "size": 100000,
+      "ferro_ta_ms": 0.1721,
+      "talib_ms": 0.62,
+      "speedup": 3.6029,
+      "ferro_ta_m_bars_s": 581.11,
+      "talib_m_bars_s": 161.29
+    },
+    {
+      "indicator": "WMA",
+      "size": 10000,
+      "ferro_ta_ms": 0.0105,
+      "talib_ms": 0.0203,
+      "speedup": 1.9363,
+      "ferro_ta_m_bars_s": 956.21,
+      "talib_m_bars_s": 493.83
+    },
+    {
+      "indicator": "WMA",
+      "size": 100000,
+      "ferro_ta_ms": 0.0869,
+      "talib_ms": 0.1868,
+      "speedup": 2.1506,
+      "ferro_ta_m_bars_s": 1151.08,
+      "talib_m_bars_s": 535.24
+    }
+  ]
+}
\ No newline at end of file
diff --git a/benchmarks/artifacts/latest/indicator_latency.json b/benchmarks/artifacts/latest/indicator_latency.json
new file mode 100644
index 0000000..c931fb1
--- /dev/null
+++ b/benchmarks/artifacts/latest/indicator_latency.json
@@ -0,0 +1,163 @@
+{
+  "metadata": {
+    "suite": "indicator_latency",
+    "runtime": {
+      "generated_at_utc": "2026-03-23T20:25:52.160357+00:00",
+      "python_version": "3.13.5",
+      "platform": "macOS-26.3.1-arm64-arm-64bit-Mach-O",
+      "machine": "arm64",
+      "processor": "arm"
+    },
+    "git": {
+      "commit": "9011250f992119170242cf17a67834c67b91bcdb",
+      "dirty": true,
+      "branch": "feat/performace-1.0.2"
+    },
+    "fixtures": [
+      {
+        "path": "/Users/pratikbhadane/Work/Projects/ferro-ta/benchmarks/fixtures/canonical_ohlcv.npz",
+        "size_bytes": 75586,
+        "sha256": "60192f8349fb06cd59ef7f70fd77aa8280399e819d7cc5eed3ca95cf5ee1a89c"
+      }
+    ],
+    "dataset": {
+      "fixture": "/Users/pratikbhadane/Work/Projects/ferro-ta/benchmarks/fixtures/canonical_ohlcv.npz",
+      "bars": 2000,
+      "rounds": 5
+    }
+  },
+  "results": [
+    {
+      "name": "VAR_20",
+      "inputs": "close",
+      "kwargs": {
+        "timeperiod": 20
+      },
+      "elapsed_ms": 0.0229
+    },
+    {
+      "name": "STOCH",
+      "inputs": "hlc",
+      "kwargs": {},
+      "elapsed_ms": 0.0201
+    },
+    {
+      "name": "WILLR_14",
+      "inputs": "hlc",
+      "kwargs": {
+        "timeperiod": 14
+      },
+      "elapsed_ms": 0.0201
+    },
+    {
+      "name": "CCI_14",
+      "inputs": "hlc",
+      "kwargs": {
+        "timeperiod": 14
+      },
+      "elapsed_ms": 0.0163
+    },
+    {
+      "name": "ADX_14",
+      "inputs": "hlc",
+      "kwargs": {
+        "timeperiod": 14
+      },
+      "elapsed_ms": 0.015
+    },
+    {
+      "name": "MACD",
+      "inputs": "close",
+      "kwargs": {},
+      "elapsed_ms": 0.0135
+    },
+    {
+      "name": "ATR_14",
+      "inputs": "hlc",
+      "kwargs": {
+        "timeperiod": 14
+      },
+      "elapsed_ms": 0.0105
+    },
+    {
+      "name": "RSI_14",
+      "inputs": "close",
+      "kwargs": {
+        "timeperiod": 14
+      },
+      "elapsed_ms": 0.0098
+    },
+    {
+      "name": "STDDEV_20",
+      "inputs": "close",
+      "kwargs": {
+        "timeperiod": 20
+      },
+      "elapsed_ms": 0.009
+    },
+    {
+      "name": "BETA_5",
+      "inputs": "pair_hl",
+      "kwargs": {
+        "timeperiod": 5
+      },
+      "elapsed_ms": 0.0083
+    },
+    {
+      "name": "CORREL_30",
+      "inputs": "pair_hl",
+      "kwargs": {
+        "timeperiod": 30
+      },
+      "elapsed_ms": 0.0076
+    },
+    {
+      "name": "BBANDS_20",
+      "inputs": "close",
+      "kwargs": {
+        "timeperiod": 20
+      },
+      "elapsed_ms": 0.0055
+    },
+    {
+      "name": "LINEARREG_14",
+      "inputs": "close",
+      "kwargs": {
+        "timeperiod": 14
+      },
+      "elapsed_ms": 0.0055
+    },
+    {
+      "name": "TSF_14",
+      "inputs": "close",
+      "kwargs": {
+        "timeperiod": 14
+      },
+      "elapsed_ms": 0.0054
+    },
+    {
+      "name": "LINEARREG_SLOPE_14",
+      "inputs": "close",
+      "kwargs": {
+        "timeperiod": 14
+      },
+      "elapsed_ms": 0.0052
+    },
+    {
+      "name": "EMA_20",
+      "inputs": "close",
+      "kwargs": {
+        "timeperiod": 20
+      },
+      "elapsed_ms": 0.005
+    },
+    {
+      "name": "SMA_20",
+      "inputs": "close",
+      "kwargs": {
+        "timeperiod": 20
+      },
+      "elapsed_ms": 0.0026
+    }
+  ]
+}
\ No newline at end of file
diff --git a/benchmarks/artifacts/latest/manifest.json b/benchmarks/artifacts/latest/manifest.json
new file mode 100644
index 0000000..a9325c5
--- /dev/null
+++ b/benchmarks/artifacts/latest/manifest.json
@@ -0,0 +1,62 @@
+{
+  "metadata": {
+    "suite": "perf_contract",
+    "runtime": {
+      "generated_at_utc": "2026-03-23T20:26:40.776130+00:00",
+      "python_version": "3.13.5",
+      "platform": "macOS-26.3.1-arm64-arm-64bit-Mach-O",
+      "machine": "arm64",
+      "processor": "arm"
+    },
+    "git": {
+      "commit": "9011250f992119170242cf17a67834c67b91bcdb",
+      "dirty": true,
+      "branch": "feat/performace-1.0.2"
+    },
+    "fixtures": [
+      {
+        "path": "/Users/pratikbhadane/Work/Projects/ferro-ta/benchmarks/fixtures/canonical_ohlcv.npz",
+        "size_bytes": 75586,
+        "sha256": "60192f8349fb06cd59ef7f70fd77aa8280399e819d7cc5eed3ca95cf5ee1a89c"
+      }
+    ],
+    "output_dir": "benchmarks/artifacts/latest"
+  },
+  "artifacts": {
+    "indicator_latency": {
+      "path": "benchmarks/artifacts/latest/indicator_latency.json",
+      "size_bytes": 3217,
+      "sha256": "43b88a50a4d7f91e30ff8e57dbf859ae5e76ecabaaf05cfcb7d8db67df920f7f"
+    },
+    "batch": {
+      "path": "benchmarks/artifacts/latest/batch.json",
+      "size_bytes": 1701,
+      "sha256": "bc900c885c48ec1903ea4870ca1de8cb9f33c609d2cefa4688cbdd18cb977f11"
+    },
+    "streaming": {
+      "path": "benchmarks/artifacts/latest/streaming.json",
+      "size_bytes": 1944,
+      "sha256": "925ba1be66d0d499daa81dfc148b03ac325ad685ce0c71e28ca1fc6927f16415"
+    },
+    "runtime_hotspots": {
+      "path": "benchmarks/artifacts/latest/runtime_hotspots.json",
+      "size_bytes": 2366,
+      "sha256": "920553b14b545f211b119c099ec59885de8b9e8056271cb2d2ac34c0c69b0906"
+    },
+    "simd": {
+      "path": "benchmarks/artifacts/latest/simd.json",
+      "size_bytes": 7700,
+      "sha256": "d48943a5dfcf4f8d8d2ca42f0004f02f9fc894de7477791b686231da665e3335"
+    },
+    "benchmark_vs_talib": {
+      "path": "benchmarks/artifacts/latest/benchmark_vs_talib.json",
+      "size_bytes": 5923,
+      "sha256": "8a4e847517f1334255353982a5266c0323bf433a1eb78dafeff808d5ad3bf7f0"
+    },
+    "wasm": {
+      "path": "benchmarks/artifacts/latest/wasm.json",
+      "size_bytes": 935,
+      "sha256": "f31fd871990c44e24a2259d618ae40a52866d20b95aa6047af3d38b9371c2ab7"
+    }
+  }
+}
\ No newline at end of file
diff --git a/benchmarks/artifacts/latest/runtime_hotspots.json b/benchmarks/artifacts/latest/runtime_hotspots.json
new file mode 100644
index 0000000..6c32fe8
--- /dev/null
+++ b/benchmarks/artifacts/latest/runtime_hotspots.json
@@ -0,0 +1,96 @@
+{
+  "metadata": {
+    "suite": "runtime_hotspots",
+    "runtime": {
+      "generated_at_utc": "2026-03-23T20:26:02.236710+00:00",
+      "python_version": "3.13.5",
+      "platform": "macOS-26.3.1-arm64-arm-64bit-Mach-O",
+      "machine": "arm64",
+      "processor": "arm"
+    },
+    "git": {
+      "commit": "9011250f992119170242cf17a67834c67b91bcdb",
+      "dirty": true,
+      "branch": "feat/performace-1.0.2"
+    },
+    "dataset": {
+      "price_bars": 20000,
+      "iv_bars": 50000,
+      "window": 252
+    }
+  },
+  "results": [
+    {
+      "category": "python_analysis",
+      "name": "iv_zscore",
+      "fast_ms": 35.984,
+      "reference_ms": 944.5804,
+      "speedup_vs_reference": 26.25,
+      "share_of_suite_pct": 77.27
+    },
+    {
+      "category": "python_analysis",
+      "name": "iv_percentile",
+      "fast_ms": 7.6624,
+      "reference_ms": 81.581,
+      "speedup_vs_reference": 10.6469,
+      "share_of_suite_pct": 16.45
+    },
+    {
+      "category": "python_analysis",
+      "name": "iv_rank",
+      "fast_ms": 2.2905,
+      "reference_ms": 198.2937,
+      "speedup_vs_reference": 86.5738,
+      "share_of_suite_pct": 4.92
+    },
+    {
+      "category": "ffi_grouping",
+      "name": "feature_matrix",
+      "fast_ms": 0.2872,
+      "reference_ms": 0.2377,
+      "speedup_vs_reference": 0.8275,
+      "share_of_suite_pct": 0.62
+    },
+    {
+      "category": "ffi_grouping",
+      "name": "compute_many_close",
+      "fast_ms": 0.1448,
+      "reference_ms": 0.1505,
+      "speedup_vs_reference": 1.0391,
+      "share_of_suite_pct": 0.31
+    },
+    {
+      "category": "rust_kernel",
+      "name": "BETA",
+      "fast_ms": 0.0637,
+      "reference_ms": 164.1752,
+      "speedup_vs_reference": 2575.2975,
+      "share_of_suite_pct": 0.14
+    },
+    {
+      "category": "rust_kernel",
+      "name": "CORREL",
+      "fast_ms": 0.0553,
+      "reference_ms": 159.6473,
+      "speedup_vs_reference": 2885.1573,
+      "share_of_suite_pct": 0.12
+    },
+    {
+      "category": "rust_kernel",
+      "name": "LINEARREG",
+      "fast_ms": 0.0415,
+      "reference_ms": 47.4665,
+      "speedup_vs_reference": 1143.77,
+      "share_of_suite_pct": 0.09
+    },
+    {
+      "category": "rust_kernel",
+      "name": "TSF",
+      "fast_ms": 0.0414,
+      "reference_ms": 47.921,
+      "speedup_vs_reference": 1157.036,
+      "share_of_suite_pct": 0.09
+    }
+  ]
+}
\ No newline at end of file
diff --git a/benchmarks/artifacts/latest/simd.json b/benchmarks/artifacts/latest/simd.json
new file mode 100644
index 0000000..01d8d65
--- /dev/null
+++ b/benchmarks/artifacts/latest/simd.json
@@ -0,0 +1,285 @@
+{
+  "metadata": {
+    "suite": "simd",
+    "runtime": {
+      "generated_at_utc": "2026-03-23T20:26:40.566511+00:00",
+      "python_version": "3.13.5",
+      "platform": "macOS-26.3.1-arm64-arm-64bit-Mach-O",
+      "machine": "arm64",
+      "processor": "arm"
+    },
+    "git": {
+      "commit": "9011250f992119170242cf17a67834c67b91bcdb",
+      "dirty": true,
+      "branch": "feat/performace-1.0.2"
+    },
+    "dataset": {
+      "price_bars": 20000,
+      "iv_bars": 50000,
+      "window": 252
+    },
+    "variants": [
+      "portable_release",
+      "simd_release"
+    ]
+  },
+  "results": [
+    {
+      "name": "BETA",
+      "category": "rust_kernel",
+      "portable_ms": 0.0635,
+      "simd_ms": 0.0636,
+      "speedup_simd_vs_portable": 0.9984
+    },
+    {
+      "name": "TSF",
+      "category": "rust_kernel",
+      "portable_ms": 0.0415,
+      "simd_ms": 0.0417,
+      "speedup_simd_vs_portable": 0.9952
+    },
+    {
+      "name": "compute_many_close",
+      "category": "ffi_grouping",
+      "portable_ms": 0.1548,
+      "simd_ms": 0.1572,
+      "speedup_simd_vs_portable": 0.9847
+    },
+    {
+      "name": "iv_zscore",
+      "category": "python_analysis",
+      "portable_ms": 36.0643,
+      "simd_ms": 37.2041,
+      "speedup_simd_vs_portable": 0.9694
+    },
+    {
+      "name": "feature_matrix",
+      "category": "ffi_grouping",
+      "portable_ms": 0.2556,
+      "simd_ms": 0.2667,
+      "speedup_simd_vs_portable": 0.9584
+    },
+    {
+      "name": "iv_percentile",
+      "category": "python_analysis",
+      "portable_ms": 7.7548,
+      "simd_ms": 8.1565,
+      "speedup_simd_vs_portable": 0.9508
+    },
+    {
+      "name": "LINEARREG",
+      "category": "rust_kernel",
+      "portable_ms": 0.0416,
+      "simd_ms": 0.0443,
+      "speedup_simd_vs_portable": 0.9391
+    },
+    {
+      "name": "iv_rank",
+      "category": "python_analysis",
+      "portable_ms": 2.2813,
+      "simd_ms": 2.4386,
+      "speedup_simd_vs_portable": 0.9355
+    },
+    {
+      "name": "CORREL",
+      "category": "rust_kernel",
+      "portable_ms": 0.0552,
+      "simd_ms": 0.0633,
+      "speedup_simd_vs_portable": 0.872
+    }
+  ],
+  "reports": {
+    "portable_release": {
+      "metadata": {
+        "suite": "runtime_hotspots",
+        "runtime": {
+          "generated_at_utc": "2026-03-23T20:26:06.920513+00:00",
+          "python_version": "3.13.5",
+          "platform": "macOS-26.3.1-arm64-arm-64bit-Mach-O",
+          "machine": "arm64",
+          "processor": "arm"
+        },
+        "git": {
+          "commit": "9011250f992119170242cf17a67834c67b91bcdb",
+          "dirty": true,
+          "branch": "feat/performace-1.0.2"
+        },
+        "dataset": {
+          "price_bars": 20000,
+          "iv_bars": 50000,
+          "window": 252
+        }
+      },
+      "results": [
+        {
+          "category": "python_analysis",
+          "name": "iv_zscore",
+          "fast_ms": 36.0643,
+          "reference_ms": 908.5403,
+          "speedup_vs_reference": 25.1922,
+          "share_of_suite_pct": 77.2
+        },
+        {
+          "category": "python_analysis",
+          "name": "iv_percentile",
+          "fast_ms": 7.7548,
+          "reference_ms": 82.2352,
+          "speedup_vs_reference": 10.6045,
+          "share_of_suite_pct": 16.6
+        },
+        {
+          "category": "python_analysis",
+          "name": "iv_rank",
+          "fast_ms": 2.2813,
+          "reference_ms": 202.5375,
+          "speedup_vs_reference": 88.7803,
+          "share_of_suite_pct": 4.88
+        },
+        {
+          "category": "ffi_grouping",
+          "name": "feature_matrix",
+          "fast_ms": 0.2556,
+          "reference_ms": 0.2252,
+          "speedup_vs_reference": 0.8812,
+          "share_of_suite_pct": 0.55
+        },
+        {
+          "category": "ffi_grouping",
+          "name": "compute_many_close",
+          "fast_ms": 0.1548,
+          "reference_ms": 0.1508,
+          "speedup_vs_reference": 0.9742,
+          "share_of_suite_pct": 0.33
+        },
+        {
+          "category": "rust_kernel",
+          "name": "BETA",
+          "fast_ms": 0.0635,
+          "reference_ms": 162.8972,
+          "speedup_vs_reference": 2563.6148,
+          "share_of_suite_pct": 0.14
+        },
+        {
+          "category": "rust_kernel",
+          "name": "CORREL",
+          "fast_ms": 0.0552,
+          "reference_ms": 163.1357,
+          "speedup_vs_reference": 2952.6826,
+          "share_of_suite_pct": 0.12
+        },
+        {
+          "category": "rust_kernel",
+          "name": "LINEARREG",
+          "fast_ms": 0.0416,
+          "reference_ms": 48.0097,
+          "speedup_vs_reference": 1153.3863,
+          "share_of_suite_pct": 0.09
+        },
+        {
+          "category": "rust_kernel",
+          "name": "TSF",
+          "fast_ms": 0.0415,
+          "reference_ms": 47.9395,
+          "speedup_vs_reference": 1155.1696,
+          "share_of_suite_pct": 0.09
+        }
+      ]
+    },
+    "simd_release": {
+      "metadata": {
+        "suite": "runtime_hotspots",
+        "runtime": {
+          "generated_at_utc": "2026-03-23T20:26:25.789478+00:00",
+          "python_version": "3.13.5",
+          "platform": "macOS-26.3.1-arm64-arm-64bit-Mach-O",
+          "machine": "arm64",
+          "processor": "arm"
+        },
+        "git": {
+          "commit": "9011250f992119170242cf17a67834c67b91bcdb",
+          "dirty": true,
+          "branch": "feat/performace-1.0.2"
+        },
+        "dataset": {
+          "price_bars": 20000,
+          "iv_bars": 50000,
+          "window": 252
+        }
+      },
+      "results": [
+        {
+          "category": "python_analysis",
+          "name": "iv_zscore",
+          "fast_ms": 37.2041,
+          "reference_ms": 930.7842,
+          "speedup_vs_reference": 25.0183,
+          "share_of_suite_pct": 76.81
+        },
+        {
+          "category": "python_analysis",
+          "name": "iv_percentile",
+          "fast_ms": 8.1565,
+          "reference_ms": 88.3639,
+          "speedup_vs_reference": 10.8336,
+          "share_of_suite_pct": 16.84
+        },
+        {
+          "category": "python_analysis",
+          "name": "iv_rank",
+          "fast_ms": 2.4386,
+          "reference_ms": 221.0389,
+          "speedup_vs_reference": 90.6424,
+          "share_of_suite_pct": 5.03
+        },
+        {
+          "category": "ffi_grouping",
+          "name": "feature_matrix",
+          "fast_ms": 0.2667,
+          "reference_ms": 0.2436,
+          "speedup_vs_reference": 0.9134,
+          "share_of_suite_pct": 0.55
+        },
+        {
+          "category": "ffi_grouping",
+          "name": "compute_many_close",
+          "fast_ms": 0.1572,
+          "reference_ms": 0.1593,
+          "speedup_vs_reference": 1.0135,
+          "share_of_suite_pct": 0.32
+        },
+        {
+          "category": "rust_kernel",
+          "name": "BETA",
+          "fast_ms": 0.0636,
+          "reference_ms": 172.9198,
+          "speedup_vs_reference": 2717.7961,
+          "share_of_suite_pct": 0.13
+        },
+        {
+          "category": "rust_kernel",
+          "name": "CORREL",
+          "fast_ms": 0.0633,
+          "reference_ms": 170.0262,
+          "speedup_vs_reference": 2686.3776,
+          "share_of_suite_pct": 0.13
+        },
+        {
+          "category": "rust_kernel",
+          "name": "LINEARREG",
+          "fast_ms": 0.0443,
+          "reference_ms": 50.5614,
+          "speedup_vs_reference": 1141.5474,
+          "share_of_suite_pct": 0.09
+        },
+        {
+          "category": "rust_kernel",
+          "name": "TSF",
+          "fast_ms": 0.0417,
+          "reference_ms": 50.9599,
+          "speedup_vs_reference": 1221.8259,
+          "share_of_suite_pct": 0.09
+        }
+      ]
+    }
+  }
+}
\ No newline at end of file
diff --git a/benchmarks/artifacts/latest/streaming.json b/benchmarks/artifacts/latest/streaming.json
new file mode 100644
index 0000000..c635596
--- /dev/null
+++ b/benchmarks/artifacts/latest/streaming.json
@@ -0,0 +1,73 @@
+{
+  "metadata": {
+    "suite": "streaming",
+    "runtime": {
+      "generated_at_utc": "2026-03-23T20:25:58.628657+00:00",
+      "python_version": "3.13.5",
+      "platform": "macOS-26.3.1-arm64-arm-64bit-Mach-O",
+      "machine": "arm64",
+      "processor": "arm"
+    },
+    "git": {
+      "commit": "9011250f992119170242cf17a67834c67b91bcdb",
+      "dirty": true,
+      "branch": "feat/performace-1.0.2"
+    },
+    "dataset": {
+      "n_bars": 100000,
+      "seed": 2026
+    }
+  },
+  "results": [
+    {
+      "indicator": "StreamingSMA",
+      "inputs": "close",
+      "stream_total_ms": 4.5729,
+      "batch_total_ms": 0.0685,
+      "stream_ns_per_update": 45.73,
+      "batch_ns_per_bar": 0.68,
+      "updates_per_second": 21867879.95,
+      "stream_over_batch_ratio": 66.7989
+    },
+    {
+      "indicator": "StreamingEMA",
+      "inputs": "close",
+      "stream_total_ms": 4.535,
+      "batch_total_ms": 0.1969,
+      "stream_ns_per_update": 45.35,
+      "batch_ns_per_bar": 1.97,
+      "updates_per_second": 22050716.65,
+      "stream_over_batch_ratio": 23.0301
+    },
+    {
+      "indicator": "StreamingRSI",
+      "inputs": "close",
+      "stream_total_ms": 4.6421,
+      "batch_total_ms": 0.4597,
+      "stream_ns_per_update": 46.42,
+      "batch_ns_per_bar": 4.6,
+      "updates_per_second": 21541858.52,
+      "stream_over_batch_ratio": 10.098
+    },
+    {
+      "indicator": "StreamingATR",
+      "inputs": "hlc",
+      "stream_total_ms": 10.2098,
+      "batch_total_ms": 0.4599,
+      "stream_ns_per_update": 102.1,
+      "batch_ns_per_bar": 4.6,
+      "updates_per_second": 9794518.83,
+      "stream_over_batch_ratio": 22.2012
+    },
+    {
+      "indicator": "StreamingVWAP",
+      "inputs": "hlcv",
+      "stream_total_ms": 12.5109,
+      "batch_total_ms": 0.1027,
+      "stream_ns_per_update": 125.11,
+      "batch_ns_per_bar": 1.03,
+      "updates_per_second": 7993046.05,
+      "stream_over_batch_ratio": 121.7603
+    }
+  ]
+}
\ No newline at end of file
diff --git a/benchmarks/artifacts/latest/wasm.json b/benchmarks/artifacts/latest/wasm.json
new file mode 100644
index 0000000..32a796a
--- /dev/null
+++ b/benchmarks/artifacts/latest/wasm.json
@@ -0,0 +1,46 @@
+{
+  "metadata": {
+    "suite": "wasm",
+    "runtime": {
+      "generated_at_utc": "2026-03-23T20:15:04.885Z",
+      "node_version": "v25.8.1",
+      "platform": "darwin",
+      "arch": "arm64"
+    },
+    "dataset": {
+      "bars": 100000
+    }
+  },
+  "results": [
+    {
+      "indicator": "SMA",
+      "elapsed_ms": 0.1702,
+      "ns_per_bar": 1.7,
+      "million_bars_per_second": 587.52
+    },
+    {
+      "indicator": "EMA",
+      "elapsed_ms": 0.2923,
+      "ns_per_bar": 2.92,
+      "million_bars_per_second": 342.08
+    },
+    {
+      "indicator": "RSI",
+      "elapsed_ms": 0.5962,
+      "ns_per_bar": 5.96,
+      "million_bars_per_second": 167.73
+    },
+    {
+      "indicator": "ATR",
+      "elapsed_ms": 0.642,
+      "ns_per_bar": 6.42,
+      "million_bars_per_second": 155.75
+    },
+    {
+      "indicator": "BBANDS",
+      "elapsed_ms": 1.7411,
+      "ns_per_bar": 17.41,
+      "million_bars_per_second": 57.43
+    }
+  ]
+}
diff --git a/benchmarks/bench_batch.py b/benchmarks/bench_batch.py
index d2335a2..31d9b65 100644
--- a/benchmarks/bench_batch.py
+++ b/benchmarks/bench_batch.py
@@ -1,65 +1,221 @@
+from __future__ import annotations
+
+import argparse
+import json
 import time
+from pathlib import Path
+from typing import Any
+
 import numpy as np
+
 import ferro_ta
 
-def _time_fn(fn, *args, **kwargs):
-    times = []
-    # Warmup
+try:
+    from benchmarks.metadata import benchmark_metadata
+except ModuleNotFoundError:  # pragma: no cover - script execution fallback
+    from metadata import benchmark_metadata
+
+
+def _time_fn(fn, *args, rounds: int = 5, **kwargs) -> float:
     fn(*args, **kwargs)
-    for _ in range(5):
+    times: list[float] = []
+    for _ in range(rounds):
         t0 = time.perf_counter()
         fn(*args, **kwargs)
         times.append(time.perf_counter() - t0)
     return min(times)
 
-def main():
-    n_samples = 100_000
-    n_series = 100
-    print(f"Batch Benchmark: {n_samples} bars, {n_series} series (Total: {n_samples*n_series/1e6:.1f} M bars)")
-    
-    np.random.seed(42)
-    # contiguous array in row-major
-    close2d = np.random.uniform(100.0, 200.0, (n_samples, n_series))
-    h2d = close2d + np.random.uniform(0.1, 2.0, (n_samples, n_series))
-    l2d = close2d - np.random.uniform(0.1, 2.0, (n_samples, n_series))
-    
-    print("-" * 50)
-    print(f"{'Indicator':<15} {'Batch (ms)':>12} {'Loop (ms)':>12} {'Speedup':>10}")
-    print("-" * 50)
 
-    # 1. SMA
-    kwargs = {"timeperiod": 14}
-    def loop_sma(arr):
-        for j in range(arr.shape[1]):
-            ferro_ta.SMA(arr[:, j], **kwargs)
-            
-    t_batch_sma = _time_fn(ferro_ta.batch.batch_sma, close2d, **kwargs)
-    t_loop_sma = _time_fn(loop_sma, close2d)
-    print(f"SMA             {t_batch_sma*1000:12.1f} {t_loop_sma*1000:12.1f} {t_loop_sma/t_batch_sma:9.1f}x")
+def run_batch_benchmark(
+    *,
+    n_samples: int = 100_000,
+    n_series: int = 100,
+    seed: int = 42,
+) -> dict[str, Any]:
+    rng = np.random.default_rng(seed)
+    close2d = rng.uniform(100.0, 200.0, (n_samples, n_series))
+    high2d = close2d + rng.uniform(0.1, 2.0, (n_samples, n_series))
+    low2d = close2d - rng.uniform(0.1, 2.0, (n_samples, n_series))
+    close1d = close2d[:, 0]
+    high1d = high2d[:, 0]
+    low1d = low2d[:, 0]
 
-    # 2. RSI
-    def loop_rsi(arr):
-        for j in range(arr.shape[1]):
-            ferro_ta.RSI(arr[:, j], **kwargs)
-    t_batch_rsi = _time_fn(ferro_ta.batch.batch_rsi, close2d, **kwargs)
-    t_loop_rsi = _time_fn(loop_rsi, close2d)
-    print(f"RSI             {t_batch_rsi*1000:12.1f} {t_loop_rsi*1000:12.1f} {t_loop_rsi/t_batch_rsi:9.1f}x")
+    batch_rows: list[dict[str, Any]] = []
+    grouped_rows: list[dict[str, Any]] = []
 
-    # 3. ATR
-    def loop_atr(h, l, c):
-        for j in range(h.shape[1]):
-            ferro_ta.ATR(h[:, j], l[:, j], c[:, j], **kwargs)
-    t_batch_atr = _time_fn(ferro_ta.batch.batch_atr, h2d, l2d, close2d, **kwargs)
-    t_loop_atr = _time_fn(loop_atr, h2d, l2d, close2d)
-    print(f"ATR             {t_batch_atr*1000:12.1f} {t_loop_atr*1000:12.1f} {t_loop_atr/t_batch_atr:9.1f}x")
+    indicators = [
+        (
+            "SMA",
+            lambda: ferro_ta.batch.batch_sma(close2d, timeperiod=14, parallel=True),
+            lambda: ferro_ta.batch.batch_sma(close2d, timeperiod=14, parallel=False),
+            lambda: [ferro_ta.SMA(close2d[:, j], timeperiod=14) for j in range(n_series)],
+        ),
+        (
+            "RSI",
+            lambda: ferro_ta.batch.batch_rsi(close2d, timeperiod=14, parallel=True),
+            lambda: ferro_ta.batch.batch_rsi(close2d, timeperiod=14, parallel=False),
+            lambda: [ferro_ta.RSI(close2d[:, j], timeperiod=14) for j in range(n_series)],
+        ),
+        (
+            "ATR",
+            lambda: ferro_ta.batch.batch_atr(
+                high2d, low2d, close2d, timeperiod=14, parallel=True
+            ),
+            lambda: ferro_ta.batch.batch_atr(
+                high2d, low2d, close2d, timeperiod=14, parallel=False
+            ),
+            lambda: [
+                ferro_ta.ATR(high2d[:, j], low2d[:, j], close2d[:, j], timeperiod=14)
+                for j in range(n_series)
+            ],
+        ),
+        (
+            "ADX",
+            lambda: ferro_ta.batch.batch_adx(
+                high2d, low2d, close2d, timeperiod=14, parallel=True
+            ),
+            lambda: ferro_ta.batch.batch_adx(
+                high2d, low2d, close2d, timeperiod=14, parallel=False
+            ),
+            lambda: [
+                ferro_ta.ADX(high2d[:, j], low2d[:, j], close2d[:, j], timeperiod=14)
+                for j in range(n_series)
+            ],
+        ),
+    ]
 
-    # 4. ADX
-    def loop_adx(h, l, c):
-        for j in range(h.shape[1]):
-            ferro_ta.ADX(h[:, j], l[:, j], c[:, j], **kwargs)
-    t_batch_adx = _time_fn(ferro_ta.batch.batch_adx, h2d, l2d, close2d, **kwargs)
-    t_loop_adx = _time_fn(loop_adx, h2d, l2d, close2d)
-    print(f"ADX             {t_batch_adx*1000:12.1f} {t_loop_adx*1000:12.1f} {t_loop_adx/t_batch_adx:9.1f}x")
+    for name, parallel_fn, sequential_fn, loop_fn in indicators:
+        batch_parallel_s = _time_fn(parallel_fn)
+        batch_sequential_s = _time_fn(sequential_fn)
+        loop_s = _time_fn(loop_fn)
+        batch_rows.append(
+            {
+                "indicator": name,
+                "parallel_ms": round(batch_parallel_s * 1000, 4),
+                "sequential_ms": round(batch_sequential_s * 1000, 4),
+                "loop_ms": round(loop_s * 1000, 4),
+                "parallel_speedup_vs_loop": round(loop_s / batch_parallel_s, 4),
+                "sequential_speedup_vs_loop": round(loop_s / batch_sequential_s, 4),
+            }
+        )
 
-if __name__ == '__main__':
-    main()
+    grouped_cases = [
+        (
+            "close_bundle_3",
+            lambda: ferro_ta.batch.compute_many(
+                [
+                    ("SMA", {"timeperiod": 10}),
+                    ("EMA", {"timeperiod": 12}),
+                    ("RSI", {"timeperiod": 14}),
+                ],
+                close=close1d,
+            ),
+            lambda: (
+                ferro_ta.SMA(close1d, timeperiod=10),
+                ferro_ta.EMA(close1d, timeperiod=12),
+                ferro_ta.RSI(close1d, timeperiod=14),
+            ),
+        ),
+        (
+            "hlc_bundle_3",
+            lambda: ferro_ta.batch.compute_many(
+                [
+                    ("ATR", {"timeperiod": 14}),
+                    ("ADX", {"timeperiod": 14}),
+                    ("CCI", {"timeperiod": 14}),
+                ],
+                close=close1d,
+                high=high1d,
+                low=low1d,
+            ),
+            lambda: (
+                ferro_ta.ATR(high1d, low1d, close1d, timeperiod=14),
+                ferro_ta.ADX(high1d, low1d, close1d, timeperiod=14),
+                ferro_ta.CCI(high1d, low1d, close1d, timeperiod=14),
+            ),
+        ),
+    ]
+
+    for name, grouped_fn, separate_fn in grouped_cases:
+        grouped_s = _time_fn(grouped_fn)
+        separate_s = _time_fn(separate_fn)
+        grouped_rows.append(
+            {
+                "case": name,
+                "grouped_ms": round(grouped_s * 1000, 4),
+                "separate_ms": round(separate_s * 1000, 4),
+                "speedup_vs_separate": round(separate_s / grouped_s, 4),
+            }
+        )
+
+    return {
+        "metadata": benchmark_metadata(
+            "batch",
+            extra={
+                "dataset": {
+                    "n_samples": n_samples,
+                    "n_series": n_series,
+                    "total_bars": n_samples * n_series,
+                    "seed": seed,
+                }
+            },
+        ),
+        "results": batch_rows,
+        "grouped_results": grouped_rows,
+    }
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(description="Benchmark batch indicator execution.")
+    parser.add_argument("--samples", type=int, default=100_000)
+    parser.add_argument("--series", type=int, default=100)
+    parser.add_argument("--seed", type=int, default=42)
+    parser.add_argument("--json", dest="json_path")
+    args = parser.parse_args()
+
+    payload = run_batch_benchmark(
+        n_samples=args.samples,
+        n_series=args.series,
+        seed=args.seed,
+    )
+
+    dataset = payload["metadata"]["dataset"]
+    print(
+        "Batch Benchmark: "
+        f"{dataset['n_samples']} bars, {dataset['n_series']} series "
+        f"(Total: {dataset['total_bars'] / 1e6:.1f} M bars)"
+    )
+    print("-" * 74)
+    print(
+        f"{'Indicator':<12} {'Parallel (ms)':>14} {'Sequential (ms)':>16} "
+        f"{'Loop (ms)':>12} {'P speedup':>10}"
+    )
+    print("-" * 74)
+    for row in payload["results"]:
+        print(
+            f"{row['indicator']:<12} {row['parallel_ms']:14.1f} "
+            f"{row['sequential_ms']:16.1f} {row['loop_ms']:12.1f} "
+            f"{row['parallel_speedup_vs_loop']:10.2f}x"
+        )
+
+    if payload["grouped_results"]:
+        print("\nGrouped Multi-Indicator Calls")
+        print("-" * 64)
+        print(f"{'Case':<18} {'Grouped (ms)':>14} {'Separate (ms)':>16} {'Speedup':>12}")
+        print("-" * 64)
+        for row in payload["grouped_results"]:
+            print(
+                f"{row['case']:<18} {row['grouped_ms']:14.1f} "
+                f"{row['separate_ms']:16.1f} {row['speedup_vs_separate']:12.2f}x"
+            )
+
+    if args.json_path:
+        json_path = Path(args.json_path)
+        json_path.write_text(json.dumps(payload, indent=2), encoding="utf-8")
+        print(f"\nWrote JSON results to {json_path}")
+
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/benchmarks/bench_simd.py b/benchmarks/bench_simd.py
new file mode 100644
index 0000000..82a27de
--- /dev/null
+++ b/benchmarks/bench_simd.py
@@ -0,0 +1,158 @@
+from __future__ import annotations
+
+import argparse
+import json
+import subprocess
+import sys
+import tempfile
+from pathlib import Path
+from typing import Any
+
+try:
+    from benchmarks.metadata import benchmark_metadata
+except ModuleNotFoundError:  # pragma: no cover - script execution fallback
+    from metadata import benchmark_metadata
+
+ROOT = Path(__file__).resolve().parents[1]
+
+
+def _run(cmd: list[str], *, cwd: Path = ROOT) -> None:
+    subprocess.run(cmd, cwd=cwd, check=True)
+
+
+def _profile_variant(
+    *,
+    label: str,
+    maturin_args: list[str],
+    price_bars: int,
+    iv_bars: int,
+    window: int,
+) -> dict[str, Any]:
+    _run([sys.executable, "-m", "maturin", "develop", "--release", *maturin_args])
+    with tempfile.TemporaryDirectory(prefix=f"ferro_ta_{label}_") as tmp_dir:
+        json_path = Path(tmp_dir) / "runtime_hotspots.json"
+        _run(
+            [
+                sys.executable,
+                "benchmarks/profile_runtime_hotspots.py",
+                "--price-bars",
+                str(price_bars),
+                "--iv-bars",
+                str(iv_bars),
+                "--window",
+                str(window),
+                "--json",
+                str(json_path),
+            ]
+        )
+        payload = json.loads(json_path.read_text(encoding="utf-8"))
+    return payload
+
+
+def run_simd_benchmark(
+    *,
+    price_bars: int = 20_000,
+    iv_bars: int = 50_000,
+    window: int = 252,
+) -> dict[str, Any]:
+    variants = [
+        ("portable_release", []),
+        ("simd_release", ["--features", "simd"]),
+    ]
+    reports = {
+        label: _profile_variant(
+            label=label,
+            maturin_args=args,
+            price_bars=price_bars,
+            iv_bars=iv_bars,
+            window=window,
+        )
+        for label, args in variants
+    }
+
+    portable_rows = {
+        row["name"]: row for row in reports["portable_release"]["results"]
+    }
+    simd_rows = {row["name"]: row for row in reports["simd_release"]["results"]}
+
+    comparison: list[dict[str, Any]] = []
+    for name in sorted(portable_rows):
+        portable = portable_rows[name]
+        simd = simd_rows.get(name)
+        if simd is None:
+            continue
+        portable_ms = float(portable["fast_ms"])
+        simd_ms = float(simd["fast_ms"])
+        comparison.append(
+            {
+                "name": name,
+                "category": portable["category"],
+                "portable_ms": round(portable_ms, 4),
+                "simd_ms": round(simd_ms, 4),
+                "speedup_simd_vs_portable": round(
+                    portable_ms / simd_ms if simd_ms > 0.0 else float("inf"), 4
+                ),
+            }
+        )
+
+    comparison.sort(
+        key=lambda row: float(row["speedup_simd_vs_portable"]), reverse=True
+    )
+
+    # Restore the default portable editable build so the workspace ends in the
+    # distributable configuration.
+    _run([sys.executable, "-m", "maturin", "develop", "--release"])
+
+    return {
+        "metadata": benchmark_metadata(
+            "simd",
+            extra={
+                "dataset": {
+                    "price_bars": price_bars,
+                    "iv_bars": iv_bars,
+                    "window": window,
+                },
+                "variants": [label for label, _ in variants],
+            },
+        ),
+        "results": comparison,
+        "reports": reports,
+    }
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(
+        description="Benchmark portable vs SIMD-enabled ferro-ta builds."
+    )
+    parser.add_argument("--price-bars", type=int, default=20_000)
+    parser.add_argument("--iv-bars", type=int, default=50_000)
+    parser.add_argument("--window", type=int, default=252)
+    parser.add_argument("--json", dest="json_path")
+    args = parser.parse_args()
+
+    payload = run_simd_benchmark(
+        price_bars=args.price_bars,
+        iv_bars=args.iv_bars,
+        window=args.window,
+    )
+
+    print(
+        f"{'Case':<20} {'Portable (ms)':>14} {'SIMD (ms)':>12} {'SIMD speedup':>14}"
+    )
+    print("-" * 64)
+    for row in payload["results"]:
+        print(
+            f"{row['name']:<20} {row['portable_ms']:14.4f} "
+            f"{row['simd_ms']:12.4f} {row['speedup_simd_vs_portable']:14.2f}x"
+        )
+
+    if args.json_path:
+        path = Path(args.json_path)
+        path.write_text(json.dumps(payload, indent=2), encoding="utf-8")
+        print(f"\nWrote JSON results to {path}")
+
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/benchmarks/bench_streaming.py b/benchmarks/bench_streaming.py
new file mode 100644
index 0000000..10a071d
--- /dev/null
+++ b/benchmarks/bench_streaming.py
@@ -0,0 +1,189 @@
+from __future__ import annotations
+
+import argparse
+import json
+import time
+from collections.abc import Callable
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+import ferro_ta as ft
+
+try:
+    from benchmarks.metadata import benchmark_metadata
+except ModuleNotFoundError:  # pragma: no cover - script execution fallback
+    from metadata import benchmark_metadata
+
+
+def _time_min(fn: Callable[[], object], rounds: int = 5) -> float:
+    fn()
+    samples: list[float] = []
+    for _ in range(rounds):
+        t0 = time.perf_counter()
+        fn()
+        samples.append(time.perf_counter() - t0)
+    return min(samples)
+
+
+def _stream_close(close: np.ndarray, factory: Callable[[], Any]) -> float:
+    streamer = factory()
+    last = np.nan
+    for value in close:
+        last = streamer.update(float(value))
+    return float(last) if not np.isnan(last) else np.nan
+
+
+def _stream_hlc(
+    high: np.ndarray,
+    low: np.ndarray,
+    close: np.ndarray,
+    factory: Callable[[], Any],
+) -> float:
+    streamer = factory()
+    last = np.nan
+    for high_value, low_value, close_value in zip(high, low, close):
+        last = streamer.update(float(high_value), float(low_value), float(close_value))
+    return float(last) if not np.isnan(last) else np.nan
+
+
+def _stream_hlcv(
+    high: np.ndarray,
+    low: np.ndarray,
+    close: np.ndarray,
+    volume: np.ndarray,
+    factory: Callable[[], Any],
+) -> float:
+    streamer = factory()
+    last = np.nan
+    for high_value, low_value, close_value, volume_value in zip(high, low, close, volume):
+        last = streamer.update(
+            float(high_value),
+            float(low_value),
+            float(close_value),
+            float(volume_value),
+        )
+    return float(last) if not np.isnan(last) else np.nan
+
+
+def run_streaming_benchmark(
+    *,
+    n_bars: int = 100_000,
+    seed: int = 2026,
+) -> dict[str, Any]:
+    rng = np.random.default_rng(seed)
+    close = 100.0 + np.cumsum(rng.normal(0.0, 1.0, n_bars)).astype(np.float64)
+    high = close + rng.uniform(0.1, 2.0, n_bars)
+    low = close - rng.uniform(0.1, 2.0, n_bars)
+    volume = rng.uniform(1_000.0, 100_000.0, n_bars)
+
+    cases = [
+        (
+            "StreamingSMA",
+            "close",
+            lambda: _stream_close(close, lambda: ft.StreamingSMA(period=20)),
+            lambda: ft.SMA(close, timeperiod=20),
+        ),
+        (
+            "StreamingEMA",
+            "close",
+            lambda: _stream_close(close, lambda: ft.StreamingEMA(period=20)),
+            lambda: ft.EMA(close, timeperiod=20),
+        ),
+        (
+            "StreamingRSI",
+            "close",
+            lambda: _stream_close(close, lambda: ft.StreamingRSI(period=14)),
+            lambda: ft.RSI(close, timeperiod=14),
+        ),
+        (
+            "StreamingATR",
+            "hlc",
+            lambda: _stream_hlc(
+                high,
+                low,
+                close,
+                lambda: ft.StreamingATR(period=14),
+            ),
+            lambda: ft.ATR(high, low, close, timeperiod=14),
+        ),
+        (
+            "StreamingVWAP",
+            "hlcv",
+            lambda: _stream_hlcv(
+                high,
+                low,
+                close,
+                volume,
+                lambda: ft.StreamingVWAP(),
+            ),
+            lambda: ft.VWAP(high, low, close, volume),
+        ),
+    ]
+
+    rows: list[dict[str, Any]] = []
+    for name, input_kind, stream_fn, batch_fn in cases:
+        stream_s = _time_min(stream_fn)
+        batch_s = _time_min(batch_fn)
+        rows.append(
+            {
+                "indicator": name,
+                "inputs": input_kind,
+                "stream_total_ms": round(stream_s * 1000.0, 4),
+                "batch_total_ms": round(batch_s * 1000.0, 4),
+                "stream_ns_per_update": round(stream_s * 1e9 / n_bars, 2),
+                "batch_ns_per_bar": round(batch_s * 1e9 / n_bars, 2),
+                "updates_per_second": round(n_bars / stream_s, 2),
+                "stream_over_batch_ratio": round(stream_s / batch_s, 4),
+            }
+        )
+
+    return {
+        "metadata": benchmark_metadata(
+            "streaming",
+            extra={
+                "dataset": {
+                    "n_bars": n_bars,
+                    "seed": seed,
+                }
+            },
+        ),
+        "results": rows,
+    }
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(description="Benchmark streaming indicator execution.")
+    parser.add_argument("--bars", type=int, default=100_000)
+    parser.add_argument("--seed", type=int, default=2026)
+    parser.add_argument("--json", dest="json_path")
+    args = parser.parse_args()
+
+    payload = run_streaming_benchmark(n_bars=args.bars, seed=args.seed)
+
+    dataset = payload["metadata"]["dataset"]
+    print(f"Streaming Benchmark: {dataset['n_bars']} bars")
+    print("-" * 86)
+    print(
+        f"{'Indicator':<16} {'Stream (ms)':>12} {'Batch (ms)':>12} "
+        f"{'ns/update':>12} {'upd/s':>12} {'ratio':>10}"
+    )
+    print("-" * 86)
+    for row in payload["results"]:
+        print(
+            f"{row['indicator']:<16} {row['stream_total_ms']:12.2f} "
+            f"{row['batch_total_ms']:12.2f} {row['stream_ns_per_update']:12.2f} "
+            f"{row['updates_per_second']:12.1f} {row['stream_over_batch_ratio']:10.2f}"
+        )
+
+    if args.json_path:
+        path = Path(args.json_path)
+        path.write_text(json.dumps(payload, indent=2), encoding="utf-8")
+        print(f"\nWrote JSON results to {path}")
+
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/benchmarks/check_hotspot_regression.py b/benchmarks/check_hotspot_regression.py
new file mode 100644
index 0000000..4902baa
--- /dev/null
+++ b/benchmarks/check_hotspot_regression.py
@@ -0,0 +1,107 @@
+#!/usr/bin/env python3
+"""
+Validate hotspot benchmark JSON against conservative speedup floors.
+
+This gate is intentionally lightweight: it checks that the optimized paths
+remain faster than their bundled reference implementations and that all
+expected cases were present in the report.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+from pathlib import Path
+
+
+def _parse_threshold_items(items: list[str]) -> dict[str, float]:
+    thresholds: dict[str, float] = {}
+    for item in items:
+        if "=" not in item:
+            raise ValueError(f"Invalid threshold '{item}', expected NAME=VALUE")
+        name, value_s = item.split("=", 1)
+        thresholds[name] = float(value_s)
+    return thresholds
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(
+        description="Check hotspot benchmark JSON against regression thresholds."
+    )
+    parser.add_argument(
+        "--input",
+        default="runtime_hotspots.json",
+        help="Path to JSON produced by benchmarks/profile_runtime_hotspots.py",
+    )
+    parser.add_argument(
+        "--min-speedup",
+        action="append",
+        default=[
+            "CORREL=2.0",
+            "BETA=2.0",
+            "LINEARREG=2.0",
+            "TSF=2.0",
+            "iv_rank=1.1",
+            "iv_percentile=1.1",
+            "iv_zscore=1.05",
+            "compute_many_close=0.85",
+            "feature_matrix=0.80",
+        ],
+        help="Required minimum speedup per named case, e.g. CORREL=5.0 (repeatable)",
+    )
+    parser.add_argument(
+        "--min-cases",
+        type=int,
+        default=9,
+        help="Minimum number of benchmark rows expected in the report",
+    )
+    args = parser.parse_args()
+
+    path = Path(args.input)
+    if not path.exists():
+        print(f"ERROR: hotspot benchmark file not found: {path}")
+        return 1
+
+    payload = json.loads(path.read_text(encoding="utf-8"))
+    rows = payload.get("results", [])
+    if len(rows) < args.min_cases:
+        print(
+            f"ERROR: hotspot report contains {len(rows)} rows, expected at least {args.min_cases}"
+        )
+        return 1
+
+    thresholds = _parse_threshold_items(args.min_speedup)
+    rows_by_name = {str(row.get("name")): row for row in rows}
+    failures: list[str] = []
+
+    for name, floor in thresholds.items():
+        row = rows_by_name.get(name)
+        if row is None:
+            failures.append(f"missing row for {name}")
+            continue
+
+        speedup = float(row.get("speedup_vs_reference", 0.0))
+        fast_ms = float(row.get("fast_ms", 0.0))
+        reference_ms = float(row.get("reference_ms", 0.0))
+        print(
+            f"{name}: fast_ms={fast_ms:.4f}, reference_ms={reference_ms:.4f}, "
+            f"speedup={speedup:.4f}"
+        )
+
+        if fast_ms <= 0.0 or reference_ms <= 0.0:
+            failures.append(f"{name} has non-positive timing values")
+        if speedup < floor:
+            failures.append(f"{name} speedup {speedup:.4f} < floor {floor:.4f}")
+
+    if failures:
+        print("FAILED hotspot regression policy:")
+        for failure in failures:
+            print(f" - {failure}")
+        return 1
+
+    print("PASS hotspot regression policy.")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/benchmarks/metadata.py b/benchmarks/metadata.py
new file mode 100644
index 0000000..cecb48d
--- /dev/null
+++ b/benchmarks/metadata.py
@@ -0,0 +1,79 @@
+from __future__ import annotations
+
+import hashlib
+import platform
+import subprocess
+import sys
+from datetime import datetime, timezone
+from pathlib import Path
+from typing import Any
+
+
+def git_info() -> dict[str, Any]:
+    """Best-effort git metadata for reproducible benchmark artifacts."""
+    try:
+        commit = subprocess.check_output(
+            ["git", "rev-parse", "HEAD"], text=True, stderr=subprocess.DEVNULL
+        ).strip()
+    except Exception:
+        commit = None
+
+    try:
+        dirty = bool(
+            subprocess.check_output(
+                ["git", "status", "--porcelain"],
+                text=True,
+                stderr=subprocess.DEVNULL,
+            ).strip()
+        )
+    except Exception:
+        dirty = None
+
+    try:
+        branch = subprocess.check_output(
+            ["git", "rev-parse", "--abbrev-ref", "HEAD"],
+            text=True,
+            stderr=subprocess.DEVNULL,
+        ).strip()
+    except Exception:
+        branch = None
+
+    return {"commit": commit, "dirty": dirty, "branch": branch}
+
+
+def runtime_info() -> dict[str, Any]:
+    return {
+        "generated_at_utc": datetime.now(timezone.utc).isoformat(),
+        "python_version": sys.version.split()[0],
+        "platform": platform.platform(),
+        "machine": platform.machine(),
+        "processor": platform.processor() or None,
+    }
+
+
+def file_info(path: str | Path) -> dict[str, Any]:
+    file_path = Path(path)
+    data = file_path.read_bytes()
+    return {
+        "path": str(file_path),
+        "size_bytes": file_path.stat().st_size,
+        "sha256": hashlib.sha256(data).hexdigest(),
+    }
+
+
+def benchmark_metadata(
+    suite: str,
+    *,
+    fixtures: list[str | Path] | None = None,
+    extra: dict[str, Any] | None = None,
+) -> dict[str, Any]:
+    metadata: dict[str, Any] = {
+        "suite": suite,
+        "runtime": runtime_info(),
+        "git": git_info(),
+    }
+    if fixtures:
+        metadata["fixtures"] = [file_info(path) for path in fixtures]
+    if extra:
+        metadata.update(extra)
+    return metadata
diff --git a/benchmarks/profile_runtime_hotspots.py b/benchmarks/profile_runtime_hotspots.py
new file mode 100644
index 0000000..ba92797
--- /dev/null
+++ b/benchmarks/profile_runtime_hotspots.py
@@ -0,0 +1,269 @@
+from __future__ import annotations
+
+import argparse
+import json
+import time
+from collections.abc import Callable
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+import ferro_ta as ft
+from ferro_ta.analysis.features import feature_matrix
+from ferro_ta.analysis.options import iv_percentile, iv_rank, iv_zscore
+from ferro_ta.data.batch import compute_many
+
+try:
+    from benchmarks.metadata import benchmark_metadata
+except ModuleNotFoundError:  # pragma: no cover - script execution fallback
+    from metadata import benchmark_metadata
+
+
+def _time_min(fn: Callable[[], object], rounds: int = 5) -> float:
+    fn()
+    samples: list[float] = []
+    for _ in range(rounds):
+        t0 = time.perf_counter()
+        fn()
+        samples.append(time.perf_counter() - t0)
+    return min(samples) * 1000.0
+
+
+def _naive_correl(x: np.ndarray, y: np.ndarray, window: int) -> np.ndarray:
+    out = np.full(len(x), np.nan, dtype=np.float64)
+    for end in range(window - 1, len(x)):
+        x_window = x[end + 1 - window : end + 1]
+        y_window = y[end + 1 - window : end + 1]
+        mean_x = float(np.sum(x_window)) / window
+        mean_y = float(np.sum(y_window)) / window
+        cov = float(np.sum((x_window - mean_x) * (y_window - mean_y)))
+        std_x = float(np.sqrt(np.sum((x_window - mean_x) ** 2)))
+        std_y = float(np.sqrt(np.sum((y_window - mean_y) ** 2)))
+        denom = std_x * std_y
+        out[end] = cov / denom if denom != 0.0 else np.nan
+    return out
+
+
+def _naive_beta(x: np.ndarray, y: np.ndarray, window: int) -> np.ndarray:
+    out = np.full(len(x), np.nan, dtype=np.float64)
+    for end in range(window, len(x)):
+        start = end - window
+        rx = np.array(
+            [x[idx + 1] / x[idx] - 1.0 if x[idx] != 0.0 else np.nan for idx in range(start, end)],
+            dtype=np.float64,
+        )
+        ry = np.array(
+            [y[idx + 1] / y[idx] - 1.0 if y[idx] != 0.0 else np.nan for idx in range(start, end)],
+            dtype=np.float64,
+        )
+        mean_x = float(np.sum(rx)) / window
+        mean_y = float(np.sum(ry)) / window
+        cov = float(np.sum((rx - mean_x) * (ry - mean_y))) / window
+        var_x = float(np.sum((rx - mean_x) ** 2)) / window
+        out[end] = cov / var_x if var_x != 0.0 else np.nan
+    return out
+
+
+def _naive_linearreg(series: np.ndarray, timeperiod: int, x_value: float) -> np.ndarray:
+    out = np.full(len(series), np.nan, dtype=np.float64)
+    xs = np.arange(timeperiod, dtype=np.float64)
+    sum_x = float(np.sum(xs))
+    sum_x2 = float(np.sum(xs * xs))
+    for end in range(timeperiod - 1, len(series)):
+        window = series[end + 1 - timeperiod : end + 1]
+        sum_y = float(np.sum(window))
+        sum_xy = float(np.sum(xs * window))
+        denom = timeperiod * sum_x2 - sum_x * sum_x
+        slope = (timeperiod * sum_xy - sum_x * sum_y) / denom if denom != 0.0 else 0.0
+        intercept = (sum_y - slope * sum_x) / timeperiod
+        out[end] = intercept + slope * x_value
+    return out
+
+
+def _old_iv_rank(iv: np.ndarray, window: int) -> np.ndarray:
+    out = np.full(len(iv), np.nan, dtype=np.float64)
+    for idx in range(window - 1, len(iv)):
+        win = iv[idx - window + 1 : idx + 1]
+        lower = float(np.nanmin(win))
+        upper = float(np.nanmax(win))
+        out[idx] = 0.0 if upper == lower else (iv[idx] - lower) / (upper - lower)
+    return out
+
+
+def _old_iv_percentile(iv: np.ndarray, window: int) -> np.ndarray:
+    out = np.full(len(iv), np.nan, dtype=np.float64)
+    for idx in range(window - 1, len(iv)):
+        win = iv[idx - window + 1 : idx + 1]
+        out[idx] = float(np.sum(win <= iv[idx])) / window
+    return out
+
+
+def _old_iv_zscore(iv: np.ndarray, window: int) -> np.ndarray:
+    out = np.full(len(iv), np.nan, dtype=np.float64)
+    for idx in range(window - 1, len(iv)):
+        win = iv[idx - window + 1 : idx + 1]
+        mean = float(np.nanmean(win))
+        std = float(np.nanstd(win, ddof=0))
+        out[idx] = np.nan if std == 0.0 else (iv[idx] - mean) / std
+    return out
+
+
+def build_hotspot_report(
+    *,
+    price_bars: int = 20_000,
+    iv_bars: int = 50_000,
+    window: int = 252,
+) -> dict[str, Any]:
+    rng = np.random.default_rng(2026)
+    close = 100 + np.cumsum(rng.normal(0, 1, price_bars)).astype(np.float64)
+    high = close + rng.uniform(0.1, 2.0, price_bars)
+    low = close - rng.uniform(0.1, 2.0, price_bars)
+    iv = rng.uniform(10.0, 40.0, iv_bars).astype(np.float64)
+    ohlcv = {"close": close, "high": high, "low": low, "volume": np.full(price_bars, 1000.0)}
+
+    rows = [
+        (
+            "rust_kernel",
+            "CORREL",
+            lambda: ft.CORREL(high, low, timeperiod=30),
+            lambda: _naive_correl(high, low, 30),
+        ),
+        (
+            "rust_kernel",
+            "BETA",
+            lambda: ft.BETA(high, low, timeperiod=5),
+            lambda: _naive_beta(high, low, 5),
+        ),
+        (
+            "rust_kernel",
+            "LINEARREG",
+            lambda: ft.LINEARREG(close, timeperiod=14),
+            lambda: _naive_linearreg(close, 14, 13.0),
+        ),
+        (
+            "rust_kernel",
+            "TSF",
+            lambda: ft.TSF(close, timeperiod=14),
+            lambda: _naive_linearreg(close, 14, 14.0),
+        ),
+        (
+            "python_analysis",
+            "iv_rank",
+            lambda: iv_rank(iv, window),
+            lambda: _old_iv_rank(iv, window),
+        ),
+        (
+            "python_analysis",
+            "iv_percentile",
+            lambda: iv_percentile(iv, window),
+            lambda: _old_iv_percentile(iv, window),
+        ),
+        (
+            "python_analysis",
+            "iv_zscore",
+            lambda: iv_zscore(iv, window),
+            lambda: _old_iv_zscore(iv, window),
+        ),
+        (
+            "ffi_grouping",
+            "compute_many_close",
+            lambda: compute_many(
+                [
+                    ("SMA", {"timeperiod": 10}),
+                    ("EMA", {"timeperiod": 12}),
+                    ("RSI", {"timeperiod": 14}),
+                ],
+                close=close,
+            ),
+            lambda: (
+                ft.SMA(close, timeperiod=10),
+                ft.EMA(close, timeperiod=12),
+                ft.RSI(close, timeperiod=14),
+            ),
+        ),
+        (
+            "ffi_grouping",
+            "feature_matrix",
+            lambda: feature_matrix(
+                ohlcv,
+                [
+                    ("SMA", {"timeperiod": 10}),
+                    ("ATR", {"timeperiod": 14}),
+                    ("ADX", {"timeperiod": 14}),
+                ],
+            ),
+            lambda: {
+                "SMA": ft.SMA(close, timeperiod=10),
+                "ATR": ft.ATR(high, low, close, timeperiod=14),
+                "ADX": ft.ADX(high, low, close, timeperiod=14),
+            },
+        ),
+    ]
+
+    results: list[dict[str, Any]] = []
+    for category, name, fast_fn, reference_fn in rows:
+        fast_ms = _time_min(fast_fn)
+        reference_ms = _time_min(reference_fn, rounds=1)
+        results.append(
+            {
+                "category": category,
+                "name": name,
+                "fast_ms": round(fast_ms, 4),
+                "reference_ms": round(reference_ms, 4),
+                "speedup_vs_reference": round(reference_ms / fast_ms, 4),
+            }
+        )
+
+    results.sort(key=lambda row: row["fast_ms"], reverse=True)
+    total_fast_ms = sum(float(row["fast_ms"]) for row in results) or 1.0
+    for row in results:
+        row["share_of_suite_pct"] = round(float(row["fast_ms"]) / total_fast_ms * 100.0, 2)
+
+    return {
+        "metadata": benchmark_metadata(
+            "runtime_hotspots",
+            extra={
+                "dataset": {
+                    "price_bars": price_bars,
+                    "iv_bars": iv_bars,
+                    "window": window,
+                }
+            },
+        ),
+        "results": results,
+    }
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(description="Profile ferro-ta runtime hotspots.")
+    parser.add_argument("--price-bars", type=int, default=20_000)
+    parser.add_argument("--iv-bars", type=int, default=50_000)
+    parser.add_argument("--window", type=int, default=252)
+    parser.add_argument("--json", dest="json_path")
+    args = parser.parse_args()
+
+    payload = build_hotspot_report(
+        price_bars=args.price_bars,
+        iv_bars=args.iv_bars,
+        window=args.window,
+    )
+
+    print(f"{'Category':<16} {'Case':<18} {'Fast (ms)':>10} {'Ref (ms)':>10} {'Speedup':>10}")
+    print("-" * 70)
+    for row in payload["results"]:
+        print(
+            f"{row['category']:<16} {row['name']:<18} {row['fast_ms']:10.2f} "
+            f"{row['reference_ms']:10.2f} {row['speedup_vs_reference']:10.2f}x"
+        )
+
+    if args.json_path:
+        path = Path(args.json_path)
+        path.write_text(json.dumps(payload, indent=2), encoding="utf-8")
+        print(f"\nWrote JSON results to {path}")
+
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/benchmarks/run_perf_contract.py b/benchmarks/run_perf_contract.py
new file mode 100644
index 0000000..682fbd7
--- /dev/null
+++ b/benchmarks/run_perf_contract.py
@@ -0,0 +1,211 @@
+from __future__ import annotations
+
+import argparse
+import json
+import time
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+try:
+    from benchmarks.bench_batch import run_batch_benchmark
+    from benchmarks.bench_simd import run_simd_benchmark
+    from benchmarks.bench_streaming import run_streaming_benchmark
+    from benchmarks.bench_vs_talib import run_comparison
+    from benchmarks.metadata import benchmark_metadata, file_info
+    from benchmarks.profile_runtime_hotspots import build_hotspot_report
+    from benchmarks.test_benchmark_suite import (
+        FIXTURE_PATH,
+        INDICATOR_SUITE,
+        _run_indicator,
+    )
+except ModuleNotFoundError:  # pragma: no cover - script execution fallback
+    from bench_batch import run_batch_benchmark
+    from bench_simd import run_simd_benchmark
+    from bench_streaming import run_streaming_benchmark
+    from bench_vs_talib import run_comparison
+    from metadata import benchmark_metadata, file_info
+    from profile_runtime_hotspots import build_hotspot_report
+    from test_benchmark_suite import FIXTURE_PATH, INDICATOR_SUITE, _run_indicator
+
+
+def _time_min(fn, rounds: int = 5) -> float:
+    fn()
+    samples: list[float] = []
+    for _ in range(rounds):
+        t0 = time.perf_counter()
+        fn()
+        samples.append(time.perf_counter() - t0)
+    return min(samples) * 1000.0
+
+
+def build_indicator_latency_report(*, rounds: int = 5) -> dict[str, Any]:
+    if not FIXTURE_PATH.exists():
+        raise FileNotFoundError(
+            f"Canonical fixture not found: {FIXTURE_PATH}. "
+            "Run benchmarks/fixtures/generate_canonical.py first."
+        )
+
+    fixture = np.load(FIXTURE_PATH)
+    ohlcv = {key: fixture[key] for key in fixture.files}
+
+    rows: list[dict[str, Any]] = []
+    for entry in INDICATOR_SUITE:
+        elapsed_ms = _time_min(lambda entry=entry: _run_indicator(entry, ohlcv), rounds=rounds)
+        rows.append(
+            {
+                "name": entry["name"],
+                "inputs": entry["inputs"],
+                "kwargs": entry["kwargs"],
+                "elapsed_ms": round(elapsed_ms, 4),
+            }
+        )
+
+    rows.sort(key=lambda row: float(row["elapsed_ms"]), reverse=True)
+    return {
+        "metadata": benchmark_metadata(
+            "indicator_latency",
+            fixtures=[FIXTURE_PATH],
+            extra={
+                "dataset": {
+                    "fixture": str(FIXTURE_PATH),
+                    "bars": len(ohlcv["close"]),
+                    "rounds": rounds,
+                }
+            },
+        ),
+        "results": rows,
+    }
+
+
+def _write_json(path: Path, payload: dict[str, Any]) -> None:
+    path.parent.mkdir(parents=True, exist_ok=True)
+    path.write_text(json.dumps(payload, indent=2), encoding="utf-8")
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(
+        description="Generate reproducible performance baseline artifacts."
+    )
+    parser.add_argument(
+        "--output-dir",
+        default="benchmarks/artifacts/latest",
+        help="Directory where benchmark JSON artifacts are written",
+    )
+    parser.add_argument("--indicator-rounds", type=int, default=5)
+    parser.add_argument("--batch-samples", type=int, default=100_000)
+    parser.add_argument("--batch-series", type=int, default=100)
+    parser.add_argument("--batch-seed", type=int, default=42)
+    parser.add_argument("--streaming-bars", type=int, default=100_000)
+    parser.add_argument("--streaming-seed", type=int, default=2026)
+    parser.add_argument("--price-bars", type=int, default=20_000)
+    parser.add_argument("--iv-bars", type=int, default=50_000)
+    parser.add_argument("--window", type=int, default=252)
+    parser.add_argument(
+        "--skip-simd",
+        action="store_true",
+        help="Skip portable-vs-SIMD comparison",
+    )
+    parser.add_argument(
+        "--talib-sizes",
+        type=int,
+        nargs="+",
+        default=[10_000, 100_000],
+        help="Bar counts used for the TA-Lib comparison suite",
+    )
+    parser.add_argument(
+        "--skip-talib",
+        action="store_true",
+        help="Skip the TA-Lib comparison artifact",
+    )
+    args = parser.parse_args()
+
+    output_dir = Path(args.output_dir)
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    artifacts: dict[str, str] = {}
+
+    indicator_path = output_dir / "indicator_latency.json"
+    _write_json(
+        indicator_path,
+        build_indicator_latency_report(rounds=args.indicator_rounds),
+    )
+    artifacts["indicator_latency"] = str(indicator_path)
+
+    batch_path = output_dir / "batch.json"
+    _write_json(
+        batch_path,
+        run_batch_benchmark(
+            n_samples=args.batch_samples,
+            n_series=args.batch_series,
+            seed=args.batch_seed,
+        ),
+    )
+    artifacts["batch"] = str(batch_path)
+
+    streaming_path = output_dir / "streaming.json"
+    _write_json(
+        streaming_path,
+        run_streaming_benchmark(
+            n_bars=args.streaming_bars,
+            seed=args.streaming_seed,
+        ),
+    )
+    artifacts["streaming"] = str(streaming_path)
+
+    hotspot_path = output_dir / "runtime_hotspots.json"
+    _write_json(
+        hotspot_path,
+        build_hotspot_report(
+            price_bars=args.price_bars,
+            iv_bars=args.iv_bars,
+            window=args.window,
+        ),
+    )
+    artifacts["runtime_hotspots"] = str(hotspot_path)
+
+    if not args.skip_simd:
+        simd_path = output_dir / "simd.json"
+        _write_json(
+            simd_path,
+            run_simd_benchmark(
+                price_bars=args.price_bars,
+                iv_bars=args.iv_bars,
+                window=args.window,
+            ),
+        )
+        artifacts["simd"] = str(simd_path)
+
+    if not args.skip_talib:
+        talib_path = output_dir / "benchmark_vs_talib.json"
+        run_comparison(args.talib_sizes, str(talib_path))
+        artifacts["benchmark_vs_talib"] = str(talib_path)
+
+    wasm_path = output_dir / "wasm.json"
+    if wasm_path.exists():
+        artifacts["wasm"] = str(wasm_path)
+
+    manifest = {
+        "metadata": benchmark_metadata(
+            "perf_contract",
+            fixtures=[FIXTURE_PATH],
+            extra={"output_dir": str(output_dir)},
+        ),
+        "artifacts": {
+            name: file_info(path)
+            for name, path in artifacts.items()
+        },
+    }
+    manifest_path = output_dir / "manifest.json"
+    _write_json(manifest_path, manifest)
+
+    print(f"Generated performance contract artifacts in {output_dir}")
+    for name, path in artifacts.items():
+        print(f" - {name}: {path}")
+    print(f" - manifest: {manifest_path}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/benchmarks/test_benchmark_suite.py b/benchmarks/test_benchmark_suite.py
index d95a7a9..05c669a 100644
--- a/benchmarks/test_benchmark_suite.py
+++ b/benchmarks/test_benchmark_suite.py
@@ -32,7 +32,8 @@ from __future__ import annotations
 
 import pathlib
 import time
-from typing import Any, Callable, Dict, List
+from collections.abc import Callable
+from typing import Any
 
 import numpy as np
 import pytest
@@ -46,7 +47,7 @@ BASELINE_PATH = pathlib.Path(__file__).parent / "baselines.npz"
 
 
 @pytest.fixture(scope="session")
-def ohlcv() -> Dict[str, np.ndarray]:
+def ohlcv() -> dict[str, np.ndarray]:
     """Load canonical OHLCV fixture."""
     if not FIXTURE_PATH.exists():
         pytest.skip(f"Canonical fixture not found: {FIXTURE_PATH}")
@@ -60,7 +61,7 @@ def ohlcv() -> Dict[str, np.ndarray]:
 
 # Each entry: (name, callable, kwargs)
 # The callable receives (close,) or (high, low, close,) based on 'inputs' key.
-INDICATOR_SUITE: List[Dict[str, Any]] = [
+INDICATOR_SUITE: list[dict[str, Any]] = [
     {
         "name": "SMA_20",
         "inputs": "close",
@@ -131,6 +132,20 @@ INDICATOR_SUITE: List[Dict[str, Any]] = [
         "fn_name": "LINEARREG",
         "kwargs": {"timeperiod": 14},
     },
+    {
+        "name": "LINEARREG_SLOPE_14",
+        "inputs": "close",
+        "fn": None,
+        "fn_name": "LINEARREG_SLOPE",
+        "kwargs": {"timeperiod": 14},
+    },
+    {
+        "name": "TSF_14",
+        "inputs": "close",
+        "fn": None,
+        "fn_name": "TSF",
+        "kwargs": {"timeperiod": 14},
+    },
     {
         "name": "VAR_20",
         "inputs": "close",
@@ -138,6 +153,20 @@ INDICATOR_SUITE: List[Dict[str, Any]] = [
         "fn_name": "VAR",
         "kwargs": {"timeperiod": 20},
     },
+    {
+        "name": "CORREL_30",
+        "inputs": "pair_hl",
+        "fn": None,
+        "fn_name": "CORREL",
+        "kwargs": {"timeperiod": 30},
+    },
+    {
+        "name": "BETA_5",
+        "inputs": "pair_hl",
+        "fn": None,
+        "fn_name": "BETA",
+        "kwargs": {"timeperiod": 5},
+    },
     {
         "name": "CCI_14",
         "inputs": "hlc",
@@ -161,12 +190,14 @@ def _load_fn(fn_name: str) -> Callable[..., Any]:
     return getattr(ft, fn_name)
 
 
-def _run_indicator(entry: Dict[str, Any], data: Dict[str, np.ndarray]) -> np.ndarray:
+def _run_indicator(entry: dict[str, Any], data: dict[str, np.ndarray]) -> np.ndarray:
     fn = _load_fn(entry["fn_name"])
     if entry["inputs"] == "close":
         result = fn(data["close"], **entry["kwargs"])
-    else:  # hlc
+    elif entry["inputs"] == "hlc":
         result = fn(data["high"], data["low"], data["close"], **entry["kwargs"])
+    else:  # pair_hl
+        result = fn(data["high"], data["low"], **entry["kwargs"])
     if isinstance(result, tuple):
         result = result[0]
     return np.asarray(result, dtype=np.float64)
@@ -184,7 +215,7 @@ class TestNumericalRegression:
         "entry", INDICATOR_SUITE, ids=[e["name"] for e in INDICATOR_SUITE]
     )
     def test_output_shape(
-        self, entry: Dict[str, Any], ohlcv: Dict[str, np.ndarray]
+        self, entry: dict[str, Any], ohlcv: dict[str, np.ndarray]
     ) -> None:
         """Indicator output length must equal input length."""
         out = _run_indicator(entry, ohlcv)
@@ -196,7 +227,7 @@ class TestNumericalRegression:
         "entry", INDICATOR_SUITE, ids=[e["name"] for e in INDICATOR_SUITE]
     )
     def test_warmup_is_nan(
-        self, entry: Dict[str, Any], ohlcv: Dict[str, np.ndarray]
+        self, entry: dict[str, Any], ohlcv: dict[str, np.ndarray]
     ) -> None:
         """First bar must be NaN (warm-up)."""
         out = _run_indicator(entry, ohlcv)
@@ -205,7 +236,7 @@ class TestNumericalRegression:
     @pytest.mark.parametrize(
         "entry", INDICATOR_SUITE, ids=[e["name"] for e in INDICATOR_SUITE]
     )
-    def test_no_inf(self, entry: Dict[str, Any], ohlcv: Dict[str, np.ndarray]) -> None:
+    def test_no_inf(self, entry: dict[str, Any], ohlcv: dict[str, np.ndarray]) -> None:
         """Output must not contain infinities."""
         out = _run_indicator(entry, ohlcv)
         assert not np.any(np.isinf(out)), f"{entry['name']}: output contains Inf"
@@ -214,7 +245,7 @@ class TestNumericalRegression:
         "entry", INDICATOR_SUITE, ids=[e["name"] for e in INDICATOR_SUITE]
     )
     def test_last_values_stable(
-        self, entry: Dict[str, Any], ohlcv: Dict[str, np.ndarray]
+        self, entry: dict[str, Any], ohlcv: dict[str, np.ndarray]
     ) -> None:
         """Last 10 non-NaN values must be finite and stable (no sudden jumps)."""
         out = _run_indicator(entry, ohlcv)
@@ -230,7 +261,7 @@ class TestNumericalRegression:
         "entry", INDICATOR_SUITE, ids=[e["name"] for e in INDICATOR_SUITE]
     )
     def test_regression_vs_baseline(
-        self, entry: Dict[str, Any], ohlcv: Dict[str, np.ndarray]
+        self, entry: dict[str, Any], ohlcv: dict[str, np.ndarray]
     ) -> None:
         """Compare last 10 values to stored baselines."""
         baselines = np.load(BASELINE_PATH)
@@ -265,8 +296,8 @@ class TestPerformance:
     )
     def test_timing(
         self,
-        entry: Dict[str, Any],
-        ohlcv: Dict[str, np.ndarray],
+        entry: dict[str, Any],
+        ohlcv: dict[str, np.ndarray],
         request: pytest.FixtureRequest,
     ) -> None:
         """Time the indicator on the canonical dataset."""
@@ -302,7 +333,7 @@ class TestPerformance:
 # ---------------------------------------------------------------------------
 
 
-def update_baselines(ohlcv_data: Dict[str, np.ndarray]) -> None:
+def update_baselines(ohlcv_data: dict[str, np.ndarray]) -> None:
     """Write current indicator outputs and timings to baselines.npz.
 
     Call this after intentional changes to update the stored baselines::
@@ -314,7 +345,7 @@ def update_baselines(ohlcv_data: Dict[str, np.ndarray]) -> None:
         update_baselines(data)
         "
     """
-    store: Dict[str, np.ndarray] = {}
+    store: dict[str, np.ndarray] = {}
     for entry in INDICATOR_SUITE:
         out = _run_indicator(entry, ohlcv_data)
         valid = out[~np.isnan(out)]
diff --git a/crates/ferro_ta_core/Cargo.toml b/crates/ferro_ta_core/Cargo.toml
index 093efd5..b432a32 100644
--- a/crates/ferro_ta_core/Cargo.toml
+++ b/crates/ferro_ta_core/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "ferro_ta_core"
-version = "1.0.1"
+version = "1.0.2"
 edition = "2021"
 description = "Pure Rust core indicator library — no PyO3, no numpy dependency"
 license = "MIT"
diff --git a/crates/ferro_ta_core/README.md b/crates/ferro_ta_core/README.md
index e4859ce..5a49414 100644
--- a/crates/ferro_ta_core/README.md
+++ b/crates/ferro_ta_core/README.md
@@ -13,7 +13,7 @@ PyO3, NumPy, or Python runtime dependency, which makes it a good fit for:
 
 ```toml
 [dependencies]
-ferro_ta_core = "1.0.1"
+ferro_ta_core = "1.0.2"
 ```
 
 ## Design
diff --git a/docs/performance.md b/docs/performance.md
index ffd32aa..4fddcf9 100644
--- a/docs/performance.md
+++ b/docs/performance.md
@@ -14,6 +14,7 @@ practical advice on how to get the best speed from the library.
 | polars users                          | Pass `pl.Series`; result is `pl.Series`| Small overhead for type conversion     |
 | Raw Rust access (expert)              | `from ferro_ta._ferro_ta import sma`     | Bypasses all Python wrappers           |
 | Multiple series at once               | `batch_sma`, `batch_ema`, `batch_rsi`  | One Python call for all columns        |
+| Many indicators on same arrays        | `compute_many`                         | Amortizes Python→Rust overhead         |
 
 **Recorded baseline and roadmap:** Performance roadmap and trade-offs are tracked
 in [PERFORMANCE_ROADMAP.md](../PERFORMANCE_ROADMAP.md). For reproducible benchmark
@@ -41,10 +42,12 @@ fast.  The bottlenecks for most users are in the Python wrapping layer:
    np.asarray(result))`), which avoids the O(n) `.tolist()` conversion of
    earlier versions.
 
-4. **Batch** — `batch_sma`/`batch_ema`/`batch_rsi` use Rust-side batch functions
-   for 2-D input (single GIL release for all columns).  The generic
-   `batch_apply` runs any indicator in a Python loop over columns; use the
-   dedicated batch functions when available.
+4. **Batch / grouped execution** — `batch_sma`/`batch_ema`/`batch_rsi` use
+   Rust-side batch functions for 2-D input (single GIL release for all
+   columns). `compute_many(...)` groups supported 1-D indicator bundles into
+   one Rust call, which helps most on medium-to-large workloads. The generic
+   `batch_apply` still runs a Python loop over columns; use it only when there
+   is no dedicated fast path.
 
 ---
 
@@ -156,6 +159,25 @@ For 2-D input, `batch_sma`/`batch_ema`/`batch_rsi` use Rust-side batch
 functions (single GIL release for all columns).  Use `batch_apply` for other
 indicators that do not have a dedicated Rust batch implementation.
 
+When you have several indicators over the same 1-D arrays, use `compute_many`:
+
+```python
+from ferro_ta.batch import compute_many
+
+results = compute_many(
+    [
+        ("SMA", {"timeperiod": 10}),
+        ("EMA", {"timeperiod": 12}),
+        ("RSI", {"timeperiod": 14}),
+    ],
+    close=close,
+)
+```
+
+Supported grouped paths currently cover common close-only indicators plus a
+small HLC bundle (`ATR`, `NATR`, `ADX`, `ADXR`, `CCI`, `WILLR`). Unsupported
+parameter shapes fall back to the normal registry path automatically.
+
 ---
 
 ## Streaming (Bar-by-Bar)
@@ -208,6 +230,53 @@ wrapper with validation and `_to_f64`; all computation runs in the extension.
 5. **Profile before optimising.** Use `cProfile` or `py-spy` to find the
    actual bottleneck before assuming a particular layer is slow.
 
+6. **Use the perf-contract scripts for evidence.** `benchmarks/run_perf_contract.py`
+   and `benchmarks/profile_runtime_hotspots.py` record timings with git/runtime
+   metadata so you can compare apples to apples across machines and commits.
+
+## Benchmark Tooling
+
+The benchmark suite now includes a small set of machine-readable scripts for
+performance work beyond the full pytest benchmark table:
+
+- `python benchmarks/bench_batch.py --json batch_benchmark.json`
+- `python benchmarks/bench_streaming.py --json streaming_benchmark.json`
+- `python benchmarks/profile_runtime_hotspots.py --json runtime_hotspots.json`
+- `python benchmarks/bench_simd.py --json simd_benchmark.json`
+- `python benchmarks/run_perf_contract.py --output-dir benchmarks/artifacts/latest`
+- `python benchmarks/check_hotspot_regression.py --input runtime_hotspots.json`
+
+The WASM bindings also ship with a Node benchmark:
+
+- `cd wasm && wasm-pack build --target nodejs --out-dir pkg`
+- `node bench.js --json ../wasm_benchmark.json`
+
+## SIMD And Build Flags
+
+Distributable wheels should stay on the portable release profile:
+
+- `cargo`/`maturin` release build
+- `lto = true`
+- `codegen-units = 1`
+- no architecture-specific `target-cpu=native` in shipped artifacts
+
+For local source builds, there are two opt-in tuning levers:
+
+```bash
+# Portable SIMD-enabled local build
+uv run maturin develop --release --features simd
+
+# Maximum local tuning for your current machine only
+RUSTFLAGS="-C target-cpu=native" uv run maturin develop --release --features simd
+```
+
+Policy:
+
+- Ship portable wheels with the default release settings.
+- Use `--features simd` for measured local/source wins.
+- Reserve `target-cpu=native` for developer workstations or private deploys,
+  because those binaries are not portable across CPU families.
+
 ---
 
 ## Performance Improvements (implemented)
@@ -247,17 +316,19 @@ bottlenecks are fixed or deferred.
   (unlike `_to_f64` for 1-D); could avoid a potential copy.
 
 **Options** (`python/ferro_ta/options.py`):
-- `iv_rank`, `iv_percentile`, `iv_zscore` use Python loops over windows
-  (O(n) iterations with per-window NumPy). Could move to Rust or vectorize.
-  See also `docs/options-volatility.md`.
+- `iv_rank`, `iv_percentile`, and `iv_zscore` are vectorized now, but
+  `iv_percentile`/`iv_zscore` still spend meaningful time in NumPy window
+  materialization on very long series.
 
 **Features** (`python/ferro_ta/features.py`):
-- With `nan_policy="fill"` and no pandas, a Python loop fills NaN per column.
-- Indicators are run in a Python loop (one call per indicator); no bulk API.
+- `nan_policy="fill"` is vectorized now.
+- `feature_matrix(...)` uses `compute_many(...)`, but grouped HLC bundles are
+  still only near parity on medium workloads and are best on larger arrays.
 
 **Signals** (`python/ferro_ta/signals.py`):
-- `compose(..., method="rank")` uses a list comprehension over columns (one
-  Python round-trip per column). Could add a Rust batch rank for 2-D input.
+- `compose(..., method="rank")` now uses a one-call Rust rank-composition
+  path, but its gains are moderate rather than dramatic. Keep measuring before
+  treating it as a major optimization lever.
 
 **Other**:
 - **dsl.py**: Some code paths use Python loops over bars.
diff --git a/pyproject.toml b/pyproject.toml
index 5ad0309..556cd98 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "maturin"
 
 [project]
 name = "ferro-ta"
-version = "1.0.1"
+version = "1.0.2"
 description = "A fast Technical Analysis library — TA-Lib alternative powered by Rust and PyO3"
 readme = "README.md"
 license = { text = "MIT" }
diff --git a/python/ferro_ta/__init__.py b/python/ferro_ta/__init__.py
index 0f042ea..8d32d3f 100644
--- a/python/ferro_ta/__init__.py
+++ b/python/ferro_ta/__init__.py
@@ -525,6 +525,7 @@ from ferro_ta.data.batch import (  # noqa: F401, E402
     batch_ema,
     batch_rsi,
     batch_sma,
+    compute_many,
 )
 from ferro_ta.data.chunked import (  # noqa: F401, E402
     chunk_apply,
diff --git a/python/ferro_ta/__init__.pyi b/python/ferro_ta/__init__.pyi
index 705d896..46fa809 100644
--- a/python/ferro_ta/__init__.pyi
+++ b/python/ferro_ta/__init__.pyi
@@ -710,6 +710,7 @@ from ferro_ta.batch import batch_apply as batch_apply
 from ferro_ta.batch import batch_ema as batch_ema
 from ferro_ta.batch import batch_rsi as batch_rsi
 from ferro_ta.batch import batch_sma as batch_sma
+from ferro_ta.batch import compute_many as compute_many
 
 # ---------------------------------------------------------------------------
 # Exception hierarchy (re-exported from ferro_ta.exceptions)
diff --git a/python/ferro_ta/analysis/backtest.py b/python/ferro_ta/analysis/backtest.py
index 6c0fc14..af4ddc3 100644
--- a/python/ferro_ta/analysis/backtest.py
+++ b/python/ferro_ta/analysis/backtest.py
@@ -360,10 +360,10 @@ def backtest(
     bar_returns[1:] = np.diff(c) / c[:-1]
 
     strategy_returns = positions * bar_returns
+    position_changed = np.concatenate([[False], positions[1:] != positions[:-1]])
 
     # Slippage: on each position change, reduce return by slippage_bps/10000 (one-way)
     if slippage_bps > 0:
-        position_changed = np.concatenate([[False], positions[1:] != positions[:-1]])
         strategy_returns = strategy_returns.copy()
         strategy_returns[position_changed] -= slippage_bps / 10_000.0
 
@@ -371,13 +371,18 @@ def backtest(
     if commission_per_trade <= 0:
         equity = np.cumprod(1.0 + strategy_returns)
     else:
-        equity = np.empty(len(c), dtype=np.float64)
-        equity[0] = 1.0
-        position_changed = np.concatenate([[False], positions[1:] != positions[:-1]])
-        for i in range(1, len(c)):
-            equity[i] = equity[i - 1] * (1.0 + strategy_returns[i])
-            if position_changed[i]:
-                equity[i] -= commission_per_trade
+        gross_equity = np.cumprod(1.0 + strategy_returns)
+        if np.any(gross_equity == 0.0):
+            equity = np.empty(len(c), dtype=np.float64)
+            equity[0] = 1.0
+            for i in range(1, len(c)):
+                equity[i] = equity[i - 1] * (1.0 + strategy_returns[i])
+                if position_changed[i]:
+                    equity[i] -= commission_per_trade
+        else:
+            commissions = position_changed.astype(np.float64) * commission_per_trade
+            discounted_commissions = np.cumsum(commissions / gross_equity)
+            equity = gross_equity * (1.0 - discounted_commissions)
 
     return BacktestResult(
         signals=signals,
diff --git a/python/ferro_ta/analysis/features.py b/python/ferro_ta/analysis/features.py
index be3d5ba..19e8be7 100644
--- a/python/ferro_ta/analysis/features.py
+++ b/python/ferro_ta/analysis/features.py
@@ -24,12 +24,23 @@ import numpy as np
 from numpy.typing import NDArray
 
 from ferro_ta._utils import _to_f64
-from ferro_ta.core.registry import run as _registry_run
+from ferro_ta.data.batch import compute_many
 
 __all__ = [
     "feature_matrix",
 ]
 
+
+def _forward_fill_nan(arr: NDArray[np.float64]) -> NDArray[np.float64]:
+    mask = np.isnan(arr)
+    if not mask.any():
+        return arr
+
+    last_valid = np.where(~mask, np.arange(len(arr)), 0)
+    np.maximum.accumulate(last_valid, out=last_valid)
+    return arr[last_valid]
+
+
 # ---------------------------------------------------------------------------
 # feature_matrix
 # ---------------------------------------------------------------------------
@@ -117,67 +128,23 @@ def feature_matrix(
     n = len(close)
     columns: dict[str, NDArray[np.float64]] = {}
 
-    # --- Indicators needing HLCV ---
-    _multi_input = {
-        "ATR",
-        "NATR",
-        "TRANGE",
-        "ADX",
-        "ADXR",
-        "PLUS_DI",
-        "MINUS_DI",
-        "PLUS_DM",
-        "MINUS_DM",
-        "DX",
-        "AROON",
-        "AROONOSC",
-        "CCI",
-        "MFI",
-        "STOCH",
-        "STOCHF",
-        "STOCHRSI",
-        "WILLR",
-        "AD",
-        "ADOSC",
-        "OBV",
-        "VWAP",
-        "DONCHIAN",
-        "ICHIMOKU",
-    }
+    results = compute_many(
+        indicators,
+        close=close,
+        high=high if high is not None else None,
+        low=low if low is not None else None,
+        volume=volume if volume is not None else None,
+    )
 
-    def _call_indicator(name: str, kwargs: dict[str, Any]) -> Any:
-        # Try with close only first; if that fails try with hlcv
-        try:
-            return _registry_run(name, close, **kwargs)
-        except (TypeError, Exception):
-            pass
-        # Build appropriate positional args from available arrays
-        if name in _multi_input and high is not None and low is not None:
-            try:
-                return _registry_run(name, high, low, close, **kwargs)
-            except Exception:
-                pass
-            if volume is not None:
-                try:
-                    return _registry_run(name, high, low, close, volume, **kwargs)
-                except Exception:
-                    pass
-        raise ValueError(
-            f"Cannot call indicator '{name}': insufficient data columns or incompatible parameters."
-        )
-
-    for spec in indicators:
+    for spec, result in zip(indicators, results):
         if isinstance(spec, str):
             name = spec
-            kwargs: dict[str, Any] = {}
             out_key: Optional[Any] = None
         elif len(spec) == 2:
-            name, kwargs = spec  # type: ignore[misc]
+            name, _ = spec  # type: ignore[misc]
             out_key = None
         else:
-            name, kwargs, out_key = spec  # type: ignore[misc]
-
-        result = _call_indicator(name, kwargs)
+            name, _, out_key = spec  # type: ignore[misc]
 
         if isinstance(result, tuple):
             if out_key is not None:
@@ -215,8 +182,6 @@ def feature_matrix(
                 mask &= ~np.isnan(arr)
             return {k: v[mask] for k, v in columns.items()}
         elif nan_policy == "fill":
-            for k, arr in columns.items():
-                for i in range(1, len(arr)):
-                    if np.isnan(arr[i]):
-                        arr[i] = arr[i - 1]
+            for key, arr in columns.items():
+                columns[key] = _forward_fill_nan(arr)
         return columns
diff --git a/python/ferro_ta/analysis/options.py b/python/ferro_ta/analysis/options.py
index 2b14620..44327e9 100644
--- a/python/ferro_ta/analysis/options.py
+++ b/python/ferro_ta/analysis/options.py
@@ -45,6 +45,7 @@ iv_zscore(iv_series, window)
 from __future__ import annotations
 
 import numpy as np
+from numpy.lib.stride_tricks import sliding_window_view
 from numpy.typing import ArrayLike, NDArray
 
 from ferro_ta.core.exceptions import FerroTAInputError, FerroTAValueError
@@ -103,15 +104,15 @@ def iv_rank(
     arr = _validate_iv(np.asarray(iv_series, dtype=np.float64), window)
     n = len(arr)
     out = np.full(n, np.nan, dtype=np.float64)
+    if window > n:
+        return out
 
-    for i in range(window - 1, n):
-        window_slice = arr[i - window + 1 : i + 1]
-        lo = float(np.nanmin(window_slice))
-        hi = float(np.nanmax(window_slice))
-        if hi == lo:
-            out[i] = 0.0
-        else:
-            out[i] = (arr[i] - lo) / (hi - lo)
+    windows = sliding_window_view(arr, window_shape=window)
+    lower = np.nanmin(windows, axis=1)
+    upper = np.nanmax(windows, axis=1)
+    current = arr[window - 1 :]
+    spread = upper - lower
+    out[window - 1 :] = np.where(spread == 0.0, 0.0, (current - lower) / spread)
 
     return out
 
@@ -149,11 +150,12 @@ def iv_percentile(
     arr = _validate_iv(np.asarray(iv_series, dtype=np.float64), window)
     n = len(arr)
     out = np.full(n, np.nan, dtype=np.float64)
+    if window > n:
+        return out
 
-    for i in range(window - 1, n):
-        window_slice = arr[i - window + 1 : i + 1]
-        current = arr[i]
-        out[i] = float(np.sum(window_slice <= current)) / window
+    windows = sliding_window_view(arr, window_shape=window)
+    current = arr[window - 1 :, None]
+    out[window - 1 :] = np.sum(windows <= current, axis=1, dtype=np.int64) / window
 
     return out
 
@@ -192,14 +194,13 @@ def iv_zscore(
     arr = _validate_iv(np.asarray(iv_series, dtype=np.float64), window)
     n = len(arr)
     out = np.full(n, np.nan, dtype=np.float64)
+    if window > n:
+        return out
 
-    for i in range(window - 1, n):
-        window_slice = arr[i - window + 1 : i + 1]
-        mu = float(np.nanmean(window_slice))
-        sigma = float(np.nanstd(window_slice, ddof=0))
-        if sigma == 0.0:
-            out[i] = np.nan
-        else:
-            out[i] = (arr[i] - mu) / sigma
+    windows = sliding_window_view(arr, window_shape=window)
+    mean = np.nanmean(windows, axis=1)
+    std = np.nanstd(windows, axis=1, ddof=0)
+    current = arr[window - 1 :]
+    out[window - 1 :] = np.where(std == 0.0, np.nan, (current - mean) / std)
 
     return out
diff --git a/python/ferro_ta/analysis/signals.py b/python/ferro_ta/analysis/signals.py
index 32e8c4f..5920e46 100644
--- a/python/ferro_ta/analysis/signals.py
+++ b/python/ferro_ta/analysis/signals.py
@@ -33,6 +33,7 @@ import numpy as np
 from numpy.typing import ArrayLike, NDArray
 
 from ferro_ta._ferro_ta import bottom_n_indices as _rust_bottom_n
+from ferro_ta._ferro_ta import compose_rank as _rust_compose_rank
 from ferro_ta._ferro_ta import compose_weighted as _rust_compose_weighted
 from ferro_ta._ferro_ta import rank_series as _rust_rank_series
 from ferro_ta._ferro_ta import top_n_indices as _rust_top_n
@@ -131,13 +132,7 @@ def compose(
         w = np.full(n_sigs, 1.0 / n_sigs)
         return _rust_compose_weighted(arr, w)
     elif method == "rank":
-        # Replace each column with its rank, then sum (ensure contiguous slices)
-        ranked = np.column_stack(
-            [_rust_rank_series(np.ascontiguousarray(arr[:, j])) for j in range(n_sigs)]
-        )
-        ranked = np.ascontiguousarray(ranked)
-        w = np.full(n_sigs, 1.0)
-        return _rust_compose_weighted(ranked, w)
+        return _rust_compose_rank(arr)
     else:
         # weighted (default)
         if weights is None:
diff --git a/python/ferro_ta/data/batch.py b/python/ferro_ta/data/batch.py
index 1c709c3..f31e92b 100644
--- a/python/ferro_ta/data/batch.py
+++ b/python/ferro_ta/data/batch.py
@@ -52,6 +52,13 @@ from ferro_ta._ferro_ta import (
 from ferro_ta._ferro_ta import (
     batch_stoch as _rust_batch_stoch,
 )
+from ferro_ta._ferro_ta import (
+    run_close_indicators as _rust_run_close_indicators,
+)
+from ferro_ta._ferro_ta import (
+    run_hlc_indicators as _rust_run_hlc_indicators,
+)
+from ferro_ta.core.registry import run as _registry_run
 from ferro_ta.indicators.momentum import RSI
 from ferro_ta.indicators.overlap import EMA, SMA
 
@@ -60,8 +67,154 @@ __all__ = [
     "batch_ema",
     "batch_rsi",
     "batch_apply",
+    "compute_many",
 ]
 
+_CLOSE_FASTPATH_DEFAULTS: dict[str, int] = {
+    "SMA": 30,
+    "EMA": 30,
+    "RSI": 14,
+    "STDDEV": 5,
+    "VAR": 5,
+    "LINEARREG": 14,
+    "LINEARREG_SLOPE": 14,
+    "LINEARREG_INTERCEPT": 14,
+    "LINEARREG_ANGLE": 14,
+    "TSF": 14,
+}
+
+_HLC_FASTPATH_DEFAULTS: dict[str, int] = {
+    "ATR": 14,
+    "NATR": 14,
+    "ADX": 14,
+    "ADXR": 14,
+    "CCI": 14,
+    "WILLR": 14,
+}
+
+
+def _normalize_indicator_spec(
+    spec: str | tuple[str, dict[str, object]] | tuple[str, dict[str, object], object],
+) -> tuple[str, dict[str, object], object | None]:
+    if isinstance(spec, str):
+        return spec, {}, None
+    if len(spec) == 2:
+        name, kwargs = spec
+        return name, kwargs, None
+    name, kwargs, out_key = spec
+    return name, kwargs, out_key
+
+
+def _extract_timeperiod(
+    name: str, kwargs: dict[str, object], defaults: dict[str, int]
+) -> int | None:
+    if name not in defaults:
+        return None
+    extra_keys = set(kwargs) - {"timeperiod"}
+    if extra_keys:
+        return None
+    raw_value = kwargs.get("timeperiod", defaults[name])
+    if not isinstance(raw_value, int):
+        return None
+    return raw_value
+
+
+def compute_many(
+    indicators: list[str | tuple[str, dict[str, object]] | tuple[str, dict[str, object], object]],
+    *,
+    close: ArrayLike,
+    high: ArrayLike | None = None,
+    low: ArrayLike | None = None,
+    volume: ArrayLike | None = None,
+    parallel: bool = True,
+) -> list[object]:
+    """Compute multiple indicators over the same arrays with grouped Rust calls.
+
+    Supported single-output indicators are grouped into one Rust boundary crossing
+    per input-shape family (`close` only or `high/low/close`). Unsupported specs
+    fall back to the regular registry path, preserving behavior.
+    """
+
+    close_arr = np.ascontiguousarray(close, dtype=np.float64)
+    high_arr = None if high is None else np.ascontiguousarray(high, dtype=np.float64)
+    low_arr = None if low is None else np.ascontiguousarray(low, dtype=np.float64)
+    volume_arr = None if volume is None else np.ascontiguousarray(volume, dtype=np.float64)
+
+    normalized = [_normalize_indicator_spec(spec) for spec in indicators]
+    results: list[object | None] = [None] * len(normalized)
+
+    close_indices: list[int] = []
+    close_names: list[str] = []
+    close_periods: list[int] = []
+
+    hlc_indices: list[int] = []
+    hlc_names: list[str] = []
+    hlc_periods: list[int] = []
+
+    for idx, (name, kwargs, out_key) in enumerate(normalized):
+        if out_key is None:
+            close_period = _extract_timeperiod(name, kwargs, _CLOSE_FASTPATH_DEFAULTS)
+            if close_period is not None:
+                close_indices.append(idx)
+                close_names.append(name)
+                close_periods.append(close_period)
+                continue
+
+            hlc_period = _extract_timeperiod(name, kwargs, _HLC_FASTPATH_DEFAULTS)
+            if (
+                hlc_period is not None
+                and high_arr is not None
+                and low_arr is not None
+            ):
+                hlc_indices.append(idx)
+                hlc_names.append(name)
+                hlc_periods.append(hlc_period)
+                continue
+
+    if close_names:
+        grouped = _rust_run_close_indicators(
+            close_arr, close_names, close_periods, parallel
+        )
+        for idx, value in zip(close_indices, grouped):
+            results[idx] = np.asarray(value, dtype=np.float64)
+
+    if hlc_names and high_arr is not None and low_arr is not None:
+        grouped = _rust_run_hlc_indicators(
+            high_arr, low_arr, close_arr, hlc_names, hlc_periods, parallel
+        )
+        for idx, value in zip(hlc_indices, grouped):
+            results[idx] = np.asarray(value, dtype=np.float64)
+
+    for idx, (name, kwargs, _) in enumerate(normalized):
+        if results[idx] is not None:
+            continue
+        try:
+            results[idx] = _registry_run(name, close_arr, **kwargs)
+            continue
+        except (TypeError, Exception):
+            pass
+
+        if high_arr is not None and low_arr is not None:
+            try:
+                results[idx] = _registry_run(name, high_arr, low_arr, close_arr, **kwargs)
+                continue
+            except Exception:
+                pass
+            if volume_arr is not None:
+                try:
+                    results[idx] = _registry_run(
+                        name, high_arr, low_arr, close_arr, volume_arr, **kwargs
+                    )
+                    continue
+                except Exception:
+                    pass
+
+        raise ValueError(
+            f"Cannot call indicator '{name}': insufficient data columns or incompatible parameters."
+        )
+
+    return [result for result in results]
+
 
 def batch_apply(
     data: ArrayLike,
diff --git a/src/batch/mod.rs b/src/batch/mod.rs
index dbec78b..db813a1 100644
--- a/src/batch/mod.rs
+++ b/src/batch/mod.rs
@@ -9,11 +9,273 @@
 //! the sequential path (`parallel = false`) may be faster due to thread-pool
 //! overhead.
 
-use ndarray::Array2;
-use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2};
+use ndarray::{Array2, ArrayView2};
+use numpy::{IntoPyArray, PyArray1, PyArray2, PyReadonlyArray1, PyReadonlyArray2};
 use pyo3::exceptions::PyValueError;
 use pyo3::prelude::*;
 use rayon::prelude::*;
+use ta::indicators::{Maximum, Minimum};
+use ta::Next;
+
+fn transpose_to_series_major(data: ArrayView2<'_, f64>) -> Array2<f64> {
+    let (n_samples, n_series) = data.dim();
+    Array2::from_shape_vec((n_series, n_samples), data.t().iter().copied().collect())
+        .expect("shape matches transposed data")
+}
+
+fn validate_same_shape(
+    expected: (usize, usize),
+    actual: (usize, usize),
+    name: &str,
+) -> PyResult<()> {
+    if actual == expected {
+        Ok(())
+    } else {
+        Err(PyValueError::new_err(format!(
+            "{name} must have shape {:?}, got {:?}",
+            expected, actual
+        )))
+    }
+}
+
+fn finish_single_output<'py>(
+    py: Python<'py>,
+    n_samples: usize,
+    n_series: usize,
+    col_results: Vec<Vec<f64>>,
+) -> Bound<'py, PyArray2<f64>> {
+    let mut result = Array2::<f64>::from_elem((n_samples, n_series), f64::NAN);
+    for (series_idx, values) in col_results.into_iter().enumerate() {
+        debug_assert_eq!(values.len(), n_samples);
+        for (sample_idx, value) in values.into_iter().enumerate() {
+            result[[sample_idx, series_idx]] = value;
+        }
+    }
+    result.into_pyarray(py)
+}
+
+fn finish_pair_output<'py>(
+    py: Python<'py>,
+    n_samples: usize,
+    n_series: usize,
+    col_results: Vec<(Vec<f64>, Vec<f64>)>,
+) -> (Bound<'py, PyArray2<f64>>, Bound<'py, PyArray2<f64>>) {
+    let mut result_k = Array2::<f64>::from_elem((n_samples, n_series), f64::NAN);
+    let mut result_d = Array2::<f64>::from_elem((n_samples, n_series), f64::NAN);
+
+    for (series_idx, (k_values, d_values)) in col_results.into_iter().enumerate() {
+        debug_assert_eq!(k_values.len(), n_samples);
+        debug_assert_eq!(d_values.len(), n_samples);
+        for (sample_idx, value) in k_values.into_iter().enumerate() {
+            result_k[[sample_idx, series_idx]] = value;
+        }
+        for (sample_idx, value) in d_values.into_iter().enumerate() {
+            result_d[[sample_idx, series_idx]] = value;
+        }
+    }
+
+    (result_k.into_pyarray(py), result_d.into_pyarray(py))
+}
+
+fn run_unary_batch<'py, F>(
+    py: Python<'py>,
+    data: PyReadonlyArray2<'py, f64>,
+    parallel: bool,
+    process_col: F,
+) -> Bound<'py, PyArray2<f64>>
+where
+    F: Fn(&[f64]) -> Vec<f64> + Sync,
+{
+    let arr = data.as_array();
+    let (n_samples, n_series) = arr.dim();
+    let series_major = transpose_to_series_major(arr);
+
+    let col_results: Vec<Vec<f64>> = py.allow_threads(|| {
+        let run = |series_idx: usize| {
+            let column_row = series_major.row(series_idx);
+            let column = column_row
+                .as_slice()
+                .expect("series-major rows are contiguous");
+            process_col(column)
+        };
+        if parallel {
+            (0..n_series).into_par_iter().map(run).collect()
+        } else {
+            (0..n_series).map(run).collect()
+        }
+    });
+
+    finish_single_output(py, n_samples, n_series, col_results)
+}
+
+fn validate_indicator_requests(names: &[String], timeperiods: &[usize]) -> PyResult<()> {
+    if names.len() != timeperiods.len() {
+        return Err(PyValueError::new_err(format!(
+            "names length ({}) must equal timeperiods length ({})",
+            names.len(),
+            timeperiods.len()
+        )));
+    }
+    for (name, &timeperiod) in names.iter().zip(timeperiods.iter()) {
+        if timeperiod == 0 {
+            return Err(PyValueError::new_err(format!(
+                "{name}: timeperiod must be >= 1"
+            )));
+        }
+    }
+    Ok(())
+}
+
+fn compute_cci(high: &[f64], low: &[f64], close: &[f64], timeperiod: usize) -> Vec<f64> {
+    let n = high.len();
+    let typical_price: Vec<f64> = high
+        .iter()
+        .zip(low.iter())
+        .zip(close.iter())
+        .map(|((&h, &l), &c)| (h + l + c) / 3.0)
+        .collect();
+
+    let mut result = vec![f64::NAN; n];
+    for end in (timeperiod - 1)..n {
+        let window = &typical_price[(end + 1 - timeperiod)..=end];
+        let mean = window.iter().sum::<f64>() / timeperiod as f64;
+        let mad = window
+            .iter()
+            .map(|&value| (value - mean).abs())
+            .sum::<f64>()
+            / timeperiod as f64;
+        result[end] = if mad != 0.0 {
+            (typical_price[end] - mean) / (0.015 * mad)
+        } else {
+            0.0
+        };
+    }
+    result
+}
+
+fn compute_willr(
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    timeperiod: usize,
+) -> PyResult<Vec<f64>> {
+    let n = high.len();
+    let mut result = vec![f64::NAN; n];
+    let mut max_ind =
+        Maximum::new(timeperiod).map_err(|err| PyValueError::new_err(err.to_string()))?;
+    let mut min_ind =
+        Minimum::new(timeperiod).map_err(|err| PyValueError::new_err(err.to_string()))?;
+
+    for (idx, ((&high_value, &low_value), &close_value)) in
+        high.iter().zip(low.iter()).zip(close.iter()).enumerate()
+    {
+        let highest = max_ind.next(high_value);
+        let lowest = min_ind.next(low_value);
+        if idx + 1 >= timeperiod {
+            let range = highest - lowest;
+            result[idx] = if range != 0.0 {
+                -100.0 * (highest - close_value) / range
+            } else {
+                -50.0
+            };
+        }
+    }
+
+    Ok(result)
+}
+
+fn compute_close_indicator(name: &str, close: &[f64], timeperiod: usize) -> PyResult<Vec<f64>> {
+    match name {
+        "SMA" => Ok(ferro_ta_core::overlap::sma(close, timeperiod)),
+        "EMA" => Ok(ferro_ta_core::overlap::ema(close, timeperiod)),
+        "RSI" => Ok(ferro_ta_core::momentum::rsi(close, timeperiod)),
+        "STDDEV" => Ok(ferro_ta_core::statistic::stddev(close, timeperiod, 1.0)),
+        "VAR" => Ok(ferro_ta_core::statistic::stddev(close, timeperiod, 1.0)
+            .into_iter()
+            .map(|value| if value.is_nan() { value } else { value * value })
+            .collect()),
+        "LINEARREG" => {
+            use crate::statistic::common::rolling_linreg_apply;
+            let last_x = (timeperiod - 1) as f64;
+            Ok(rolling_linreg_apply(
+                close,
+                timeperiod,
+                |slope: f64, intercept: f64| intercept + slope * last_x,
+            ))
+        }
+        "LINEARREG_SLOPE" => {
+            use crate::statistic::common::rolling_linreg_apply;
+            Ok(rolling_linreg_apply(
+                close,
+                timeperiod,
+                |slope: f64, _: f64| slope,
+            ))
+        }
+        "LINEARREG_INTERCEPT" => {
+            use crate::statistic::common::rolling_linreg_apply;
+            Ok(rolling_linreg_apply(
+                close,
+                timeperiod,
+                |_: f64, intercept: f64| intercept,
+            ))
+        }
+        "LINEARREG_ANGLE" => {
+            use crate::statistic::common::rolling_linreg_apply;
+            Ok(rolling_linreg_apply(
+                close,
+                timeperiod,
+                |slope: f64, _: f64| slope.atan() * 180.0 / std::f64::consts::PI,
+            ))
+        }
+        "TSF" => {
+            use crate::statistic::common::rolling_linreg_apply;
+            let forecast_x = timeperiod as f64;
+            Ok(rolling_linreg_apply(
+                close,
+                timeperiod,
+                |slope: f64, intercept: f64| intercept + slope * forecast_x,
+            ))
+        }
+        _ => Err(PyValueError::new_err(format!(
+            "unsupported close indicator for grouped execution: {name}"
+        ))),
+    }
+}
+
+fn compute_hlc_indicator(
+    name: &str,
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    timeperiod: usize,
+) -> PyResult<Vec<f64>> {
+    match name {
+        "ATR" => Ok(ferro_ta_core::volatility::atr(high, low, close, timeperiod)),
+        "NATR" => {
+            let atr = ferro_ta_core::volatility::atr(high, low, close, timeperiod);
+            Ok(atr
+                .into_iter()
+                .zip(close.iter())
+                .map(|(atr_value, &close_value)| {
+                    if atr_value.is_nan() || close_value == 0.0 {
+                        f64::NAN
+                    } else {
+                        (atr_value / close_value) * 100.0
+                    }
+                })
+                .collect())
+        }
+        "ADX" => Ok(ferro_ta_core::momentum::adx(high, low, close, timeperiod)),
+        "ADXR" => Ok(ferro_ta_core::momentum::adxr(high, low, close, timeperiod)),
+        "CCI" => Ok(compute_cci(high, low, close, timeperiod)),
+        "WILLR" => compute_willr(high, low, close, timeperiod),
+        _ => Err(PyValueError::new_err(format!(
+            "unsupported HLC indicator for grouped execution: {name}"
+        ))),
+    }
+}
+
+type IndicatorArrayList = Vec<Py<PyArray1<f64>>>;
 
 // ---------------------------------------------------------------------------
 // batch_sma
@@ -43,34 +305,13 @@ pub fn batch_sma<'py>(
     if timeperiod == 0 {
         return Err(PyValueError::new_err("timeperiod must be >= 1"));
     }
-    let arr = data.as_array();
-    let (n_samples, n_series) = arr.dim();
+    let (n_samples, n_series) = data.as_array().dim();
     log::debug!(
         "batch_sma: timeperiod={timeperiod}, shape=({n_samples}, {n_series}), parallel={parallel}"
     );
-
-    // Extract columns to owned Vecs so we can release the GIL for parallel work.
-    let columns: Vec<Vec<f64>> = (0..n_series)
-        .map(|j| (0..n_samples).map(|i| arr[[i, j]]).collect())
-        .collect();
-
-    let process_col = |col: &Vec<f64>| -> Vec<f64> { ferro_ta_core::overlap::sma(col, timeperiod) };
-
-    let col_results: Vec<Vec<f64>> = py.allow_threads(|| {
-        if parallel {
-            columns.par_iter().map(process_col).collect()
-        } else {
-            columns.iter().map(process_col).collect()
-        }
-    });
-
-    let mut result = Array2::<f64>::from_elem((n_samples, n_series), f64::NAN);
-    for (j, col_result) in col_results.iter().enumerate() {
-        for (i, &val) in col_result.iter().enumerate() {
-            result[[i, j]] = val;
-        }
-    }
-    Ok(result.into_pyarray(py))
+    Ok(run_unary_batch(py, data, parallel, |col| {
+        ferro_ta_core::overlap::sma(col, timeperiod)
+    }))
 }
 
 // ---------------------------------------------------------------------------
@@ -100,33 +341,13 @@ pub fn batch_ema<'py>(
     if timeperiod == 0 {
         return Err(PyValueError::new_err("timeperiod must be >= 1"));
     }
-    let arr = data.as_array();
-    let (n_samples, n_series) = arr.dim();
+    let (n_samples, n_series) = data.as_array().dim();
     log::debug!(
         "batch_ema: timeperiod={timeperiod}, shape=({n_samples}, {n_series}), parallel={parallel}"
     );
-
-    let columns: Vec<Vec<f64>> = (0..n_series)
-        .map(|j| (0..n_samples).map(|i| arr[[i, j]]).collect())
-        .collect();
-
-    let process_col = |col: &Vec<f64>| -> Vec<f64> { ferro_ta_core::overlap::ema(col, timeperiod) };
-
-    let col_results: Vec<Vec<f64>> = py.allow_threads(|| {
-        if parallel {
-            columns.par_iter().map(process_col).collect()
-        } else {
-            columns.iter().map(process_col).collect()
-        }
-    });
-
-    let mut result = Array2::<f64>::from_elem((n_samples, n_series), f64::NAN);
-    for (j, col_result) in col_results.iter().enumerate() {
-        for (i, &val) in col_result.iter().enumerate() {
-            result[[i, j]] = val;
-        }
-    }
-    Ok(result.into_pyarray(py))
+    Ok(run_unary_batch(py, data, parallel, |col| {
+        ferro_ta_core::overlap::ema(col, timeperiod)
+    }))
 }
 
 // ---------------------------------------------------------------------------
@@ -157,18 +378,13 @@ pub fn batch_rsi<'py>(
     if timeperiod == 0 {
         return Err(PyValueError::new_err("timeperiod must be >= 1"));
     }
-    let arr = data.as_array();
-    let (n_samples, n_series) = arr.dim();
+    let (n_samples, n_series) = data.as_array().dim();
     log::debug!(
         "batch_rsi: timeperiod={timeperiod}, shape=({n_samples}, {n_series}), parallel={parallel}"
     );
 
-    let columns: Vec<Vec<f64>> = (0..n_series)
-        .map(|j| (0..n_samples).map(|i| arr[[i, j]]).collect())
-        .collect();
-
     let period_f = timeperiod as f64;
-    let process_col = |col: &Vec<f64>| -> Vec<f64> {
+    Ok(run_unary_batch(py, data, parallel, |col| {
         let mut col_result = vec![f64::NAN; n_samples];
         if n_samples <= timeperiod {
             return col_result;
@@ -208,23 +424,7 @@ pub fn batch_rsi<'py>(
             col_result[i] = 100.0 - 100.0 / (1.0 + rs);
         }
         col_result
-    };
-
-    let col_results: Vec<Vec<f64>> = py.allow_threads(|| {
-        if parallel {
-            columns.par_iter().map(process_col).collect()
-        } else {
-            columns.iter().map(process_col).collect()
-        }
-    });
-
-    let mut result = Array2::<f64>::from_elem((n_samples, n_series), f64::NAN);
-    for (j, col_result) in col_results.iter().enumerate() {
-        for (i, &val) in col_result.iter().enumerate() {
-            result[[i, j]] = val;
-        }
-    }
-    Ok(result.into_pyarray(py))
+    }))
 }
 
 // ---------------------------------------------------------------------------
@@ -248,20 +448,28 @@ pub fn batch_atr<'py>(
     let arr_l = low.as_array();
     let arr_c = close.as_array();
     let (n_samples, n_series) = arr_h.dim();
+    validate_same_shape((n_samples, n_series), arr_l.dim(), "low")?;
+    validate_same_shape((n_samples, n_series), arr_c.dim(), "close")?;
 
-    let cols_h: Vec<Vec<f64>> = (0..n_series)
-        .map(|j| (0..n_samples).map(|i| arr_h[[i, j]]).collect())
-        .collect();
-    let cols_l: Vec<Vec<f64>> = (0..n_series)
-        .map(|j| (0..n_samples).map(|i| arr_l[[i, j]]).collect())
-        .collect();
-    let cols_c: Vec<Vec<f64>> = (0..n_series)
-        .map(|j| (0..n_samples).map(|i| arr_c[[i, j]]).collect())
-        .collect();
+    let high_by_series = transpose_to_series_major(arr_h);
+    let low_by_series = transpose_to_series_major(arr_l);
+    let close_by_series = transpose_to_series_major(arr_c);
 
     let col_results: Vec<Vec<f64>> = py.allow_threads(|| {
-        let process_col = |j: usize| -> Vec<f64> {
-            ferro_ta_core::volatility::atr(&cols_h[j], &cols_l[j], &cols_c[j], timeperiod)
+        let process_col = |series_idx: usize| -> Vec<f64> {
+            let high_row = high_by_series.row(series_idx);
+            let low_row = low_by_series.row(series_idx);
+            let close_row = close_by_series.row(series_idx);
+            let high_col = high_row
+                .as_slice()
+                .expect("series-major rows are contiguous");
+            let low_col = low_row
+                .as_slice()
+                .expect("series-major rows are contiguous");
+            let close_col = close_row
+                .as_slice()
+                .expect("series-major rows are contiguous");
+            ferro_ta_core::volatility::atr(high_col, low_col, close_col, timeperiod)
         };
         if parallel {
             (0..n_series).into_par_iter().map(process_col).collect()
@@ -269,14 +477,7 @@ pub fn batch_atr<'py>(
             (0..n_series).map(process_col).collect()
         }
     });
-
-    let mut result = Array2::<f64>::from_elem((n_samples, n_series), f64::NAN);
-    for (j, col_result) in col_results.iter().enumerate() {
-        for (i, &val) in col_result.iter().enumerate() {
-            result[[i, j]] = val;
-        }
-    }
-    Ok(result.into_pyarray(py))
+    Ok(finish_single_output(py, n_samples, n_series, col_results))
 }
 
 // ---------------------------------------------------------------------------
@@ -303,23 +504,31 @@ pub fn batch_stoch<'py>(
     let arr_l = low.as_array();
     let arr_c = close.as_array();
     let (n_samples, n_series) = arr_h.dim();
+    validate_same_shape((n_samples, n_series), arr_l.dim(), "low")?;
+    validate_same_shape((n_samples, n_series), arr_c.dim(), "close")?;
 
-    let cols_h: Vec<Vec<f64>> = (0..n_series)
-        .map(|j| (0..n_samples).map(|i| arr_h[[i, j]]).collect())
-        .collect();
-    let cols_l: Vec<Vec<f64>> = (0..n_series)
-        .map(|j| (0..n_samples).map(|i| arr_l[[i, j]]).collect())
-        .collect();
-    let cols_c: Vec<Vec<f64>> = (0..n_series)
-        .map(|j| (0..n_samples).map(|i| arr_c[[i, j]]).collect())
-        .collect();
+    let high_by_series = transpose_to_series_major(arr_h);
+    let low_by_series = transpose_to_series_major(arr_l);
+    let close_by_series = transpose_to_series_major(arr_c);
 
     let col_results: Vec<(Vec<f64>, Vec<f64>)> = py.allow_threads(|| {
-        let process_col = |j: usize| -> (Vec<f64>, Vec<f64>) {
+        let process_col = |series_idx: usize| -> (Vec<f64>, Vec<f64>) {
+            let high_row = high_by_series.row(series_idx);
+            let low_row = low_by_series.row(series_idx);
+            let close_row = close_by_series.row(series_idx);
+            let high_col = high_row
+                .as_slice()
+                .expect("series-major rows are contiguous");
+            let low_col = low_row
+                .as_slice()
+                .expect("series-major rows are contiguous");
+            let close_col = close_row
+                .as_slice()
+                .expect("series-major rows are contiguous");
             ferro_ta_core::momentum::stoch(
-                &cols_h[j],
-                &cols_l[j],
-                &cols_c[j],
+                high_col,
+                low_col,
+                close_col,
                 fastk_period,
                 slowk_period,
                 slowd_period,
@@ -331,16 +540,7 @@ pub fn batch_stoch<'py>(
             (0..n_series).map(process_col).collect()
         }
     });
-
-    let mut result_k = Array2::<f64>::from_elem((n_samples, n_series), f64::NAN);
-    let mut result_d = Array2::<f64>::from_elem((n_samples, n_series), f64::NAN);
-    for (j, (k_col, d_col)) in col_results.iter().enumerate() {
-        for i in 0..n_samples {
-            result_k[[i, j]] = k_col[i];
-            result_d[[i, j]] = d_col[i];
-        }
-    }
-    Ok((result_k.into_pyarray(py), result_d.into_pyarray(py)))
+    Ok(finish_pair_output(py, n_samples, n_series, col_results))
 }
 
 // ---------------------------------------------------------------------------
@@ -364,20 +564,28 @@ pub fn batch_adx<'py>(
     let arr_l = low.as_array();
     let arr_c = close.as_array();
     let (n_samples, n_series) = arr_h.dim();
+    validate_same_shape((n_samples, n_series), arr_l.dim(), "low")?;
+    validate_same_shape((n_samples, n_series), arr_c.dim(), "close")?;
 
-    let cols_h: Vec<Vec<f64>> = (0..n_series)
-        .map(|j| (0..n_samples).map(|i| arr_h[[i, j]]).collect())
-        .collect();
-    let cols_l: Vec<Vec<f64>> = (0..n_series)
-        .map(|j| (0..n_samples).map(|i| arr_l[[i, j]]).collect())
-        .collect();
-    let cols_c: Vec<Vec<f64>> = (0..n_series)
-        .map(|j| (0..n_samples).map(|i| arr_c[[i, j]]).collect())
-        .collect();
+    let high_by_series = transpose_to_series_major(arr_h);
+    let low_by_series = transpose_to_series_major(arr_l);
+    let close_by_series = transpose_to_series_major(arr_c);
 
     let col_results: Vec<Vec<f64>> = py.allow_threads(|| {
-        let process_col = |j: usize| -> Vec<f64> {
-            ferro_ta_core::momentum::adx(&cols_h[j], &cols_l[j], &cols_c[j], timeperiod)
+        let process_col = |series_idx: usize| -> Vec<f64> {
+            let high_row = high_by_series.row(series_idx);
+            let low_row = low_by_series.row(series_idx);
+            let close_row = close_by_series.row(series_idx);
+            let high_col = high_row
+                .as_slice()
+                .expect("series-major rows are contiguous");
+            let low_col = low_row
+                .as_slice()
+                .expect("series-major rows are contiguous");
+            let close_col = close_row
+                .as_slice()
+                .expect("series-major rows are contiguous");
+            ferro_ta_core::momentum::adx(high_col, low_col, close_col, timeperiod)
         };
         if parallel {
             (0..n_series).into_par_iter().map(process_col).collect()
@@ -385,14 +593,83 @@ pub fn batch_adx<'py>(
             (0..n_series).map(process_col).collect()
         }
     });
+    Ok(finish_single_output(py, n_samples, n_series, col_results))
+}
 
-    let mut result = Array2::<f64>::from_elem((n_samples, n_series), f64::NAN);
-    for (j, col_result) in col_results.iter().enumerate() {
-        for (i, &val) in col_result.iter().enumerate() {
-            result[[i, j]] = val;
+// ---------------------------------------------------------------------------
+// grouped 1-D execution
+// ---------------------------------------------------------------------------
+
+#[pyfunction]
+#[pyo3(signature = (close, names, timeperiods, parallel = true))]
+pub fn run_close_indicators<'py>(
+    py: Python<'py>,
+    close: PyReadonlyArray1<'py, f64>,
+    names: Vec<String>,
+    timeperiods: Vec<usize>,
+    parallel: bool,
+) -> PyResult<IndicatorArrayList> {
+    validate_indicator_requests(&names, &timeperiods)?;
+    let close_values = close.as_slice()?;
+
+    let results: Vec<PyResult<Vec<f64>>> = py.allow_threads(|| {
+        let run = |idx: usize| compute_close_indicator(&names[idx], close_values, timeperiods[idx]);
+        if parallel {
+            (0..names.len()).into_par_iter().map(run).collect()
+        } else {
+            (0..names.len()).map(run).collect()
         }
+    });
+
+    results
+        .into_iter()
+        .map(|result| result.map(|values| values.into_pyarray(py).unbind()))
+        .collect()
+}
+
+#[pyfunction]
+#[pyo3(signature = (high, low, close, names, timeperiods, parallel = true))]
+pub fn run_hlc_indicators<'py>(
+    py: Python<'py>,
+    high: PyReadonlyArray1<'py, f64>,
+    low: PyReadonlyArray1<'py, f64>,
+    close: PyReadonlyArray1<'py, f64>,
+    names: Vec<String>,
+    timeperiods: Vec<usize>,
+    parallel: bool,
+) -> PyResult<IndicatorArrayList> {
+    validate_indicator_requests(&names, &timeperiods)?;
+    let high_values = high.as_slice()?;
+    let low_values = low.as_slice()?;
+    let close_values = close.as_slice()?;
+
+    if high_values.len() != low_values.len() || high_values.len() != close_values.len() {
+        return Err(PyValueError::new_err(
+            "high, low, and close must have equal length",
+        ));
     }
-    Ok(result.into_pyarray(py))
+
+    let results: Vec<PyResult<Vec<f64>>> = py.allow_threads(|| {
+        let run = |idx: usize| {
+            compute_hlc_indicator(
+                &names[idx],
+                high_values,
+                low_values,
+                close_values,
+                timeperiods[idx],
+            )
+        };
+        if parallel {
+            (0..names.len()).into_par_iter().map(run).collect()
+        } else {
+            (0..names.len()).map(run).collect()
+        }
+    });
+
+    results
+        .into_iter()
+        .map(|result| result.map(|values| values.into_pyarray(py).unbind()))
+        .collect()
 }
 
 // ---------------------------------------------------------------------------
@@ -406,5 +683,7 @@ pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(pyo3::wrap_pyfunction!(batch_atr, m)?)?;
     m.add_function(pyo3::wrap_pyfunction!(batch_stoch, m)?)?;
     m.add_function(pyo3::wrap_pyfunction!(batch_adx, m)?)?;
+    m.add_function(pyo3::wrap_pyfunction!(run_close_indicators, m)?)?;
+    m.add_function(pyo3::wrap_pyfunction!(run_hlc_indicators, m)?)?;
     Ok(())
 }
diff --git a/src/signals/mod.rs b/src/signals/mod.rs
index a706cf8..da55cf0 100644
--- a/src/signals/mod.rs
+++ b/src/signals/mod.rs
@@ -4,10 +4,36 @@
 //! - `top_n_indices`    — indices of the N largest values in a 1-D array
 //! - `bottom_n_indices` — indices of the N smallest values in a 1-D array
 
-use numpy::{IntoPyArray, PyArray1, PyReadonlyArray1};
+use numpy::{IntoPyArray, PyArray1, PyReadonlyArray1, PyReadonlyArray2};
 use pyo3::exceptions::PyValueError;
 use pyo3::prelude::*;
 
+fn rank_values(xv: &[f64]) -> Vec<f64> {
+    let n = xv.len();
+    let mut order: Vec<usize> = (0..n).collect();
+    order.sort_by(|&a, &b| {
+        xv[a]
+            .partial_cmp(&xv[b])
+            .unwrap_or(std::cmp::Ordering::Equal)
+    });
+
+    let mut ranks = vec![0.0_f64; n];
+    let mut i = 0;
+    while i < n {
+        let val = xv[order[i]];
+        let mut j = i + 1;
+        while j < n && xv[order[j]] == val {
+            j += 1;
+        }
+        let avg_rank = (i + 1 + j) as f64 / 2.0;
+        for k in i..j {
+            ranks[order[k]] = avg_rank;
+        }
+        i = j;
+    }
+    ranks
+}
+
 // ---------------------------------------------------------------------------
 // rank_series
 // ---------------------------------------------------------------------------
@@ -33,30 +59,43 @@ pub fn rank_series<'py>(
     if n == 0 {
         return Err(PyValueError::new_err("x must be non-empty"));
     }
-    // Sort indices by value
-    let mut order: Vec<usize> = (0..n).collect();
-    order.sort_by(|&a, &b| {
-        xv[a]
-            .partial_cmp(&xv[b])
-            .unwrap_or(std::cmp::Ordering::Equal)
+    Ok(rank_values(xv).into_pyarray(py))
+}
+
+// ---------------------------------------------------------------------------
+// compose_rank
+// ---------------------------------------------------------------------------
+
+/// Compute rank-based composite scores for a 2-D signal matrix.
+///
+/// Each column is ranked independently (ascending, fractional ranks for ties),
+/// and the per-row ranks are summed across columns.
+#[pyfunction]
+pub fn compose_rank<'py>(
+    py: Python<'py>,
+    signals: PyReadonlyArray2<'py, f64>,
+) -> PyResult<Bound<'py, PyArray1<f64>>> {
+    let arr = signals.as_array();
+    let (n_bars, n_sigs) = arr.dim();
+    if n_bars == 0 || n_sigs == 0 {
+        return Err(PyValueError::new_err(
+            "signals must be a non-empty 2-D array",
+        ));
+    }
+
+    let scores = py.allow_threads(|| {
+        let mut scores = vec![0.0_f64; n_bars];
+        for sig_idx in 0..n_sigs {
+            let column: Vec<f64> = arr.column(sig_idx).iter().copied().collect();
+            let ranks = rank_values(&column);
+            for (bar_idx, rank) in ranks.into_iter().enumerate() {
+                scores[bar_idx] += rank;
+            }
+        }
+        scores
     });
 
-    let mut ranks = vec![0.0_f64; n];
-    let mut i = 0;
-    while i < n {
-        let val = xv[order[i]];
-        let mut j = i + 1;
-        while j < n && xv[order[j]] == val {
-            j += 1;
-        }
-        // Positions [i..j) all have the same value; average rank = (i+1 + j)/2
-        let avg_rank = (i + 1 + j) as f64 / 2.0;
-        for k in i..j {
-            ranks[order[k]] = avg_rank;
-        }
-        i = j;
-    }
-    Ok(ranks.into_pyarray(py))
+    Ok(scores.into_pyarray(py))
 }
 
 // ---------------------------------------------------------------------------
@@ -122,6 +161,7 @@ pub fn bottom_n_indices<'py>(
 
 pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(rank_series, m)?)?;
+    m.add_function(wrap_pyfunction!(compose_rank, m)?)?;
     m.add_function(wrap_pyfunction!(top_n_indices, m)?)?;
     m.add_function(wrap_pyfunction!(bottom_n_indices, m)?)?;
     Ok(())
diff --git a/src/statistic/beta.rs b/src/statistic/beta.rs
index f97484b..39d3097 100644
--- a/src/statistic/beta.rs
+++ b/src/statistic/beta.rs
@@ -2,6 +2,45 @@ use crate::validation;
 use numpy::{IntoPyArray, PyArray1, PyReadonlyArray1};
 use pyo3::prelude::*;
 
+fn price_return(curr: f64, prev: f64) -> f64 {
+    if prev != 0.0 {
+        curr / prev - 1.0
+    } else {
+        f64::NAN
+    }
+}
+
+fn beta_fallback(x: &[f64], y: &[f64], timeperiod: usize) -> Vec<f64> {
+    let n = x.len();
+    let mut result = vec![f64::NAN; n];
+    for end in timeperiod..n {
+        let start = end - timeperiod;
+        let mut rx = vec![0.0_f64; timeperiod];
+        let mut ry = vec![0.0_f64; timeperiod];
+        for offset in 0..timeperiod {
+            let prev = start + offset;
+            let curr = prev + 1;
+            rx[offset] = price_return(x[curr], x[prev]);
+            ry[offset] = price_return(y[curr], y[prev]);
+        }
+        let mean_x = rx.iter().sum::<f64>() / timeperiod as f64;
+        let mean_y = ry.iter().sum::<f64>() / timeperiod as f64;
+        let cov = rx
+            .iter()
+            .zip(ry.iter())
+            .map(|(&lhs, &rhs)| (lhs - mean_x) * (rhs - mean_y))
+            .sum::<f64>()
+            / timeperiod as f64;
+        let var_x = rx
+            .iter()
+            .map(|&value| (value - mean_x).powi(2))
+            .sum::<f64>()
+            / timeperiod as f64;
+        result[end] = if var_x != 0.0 { cov / var_x } else { f64::NAN };
+    }
+    result
+}
+
 /// Beta: regression of *real1* daily returns on *real0* daily returns over a
 /// rolling window of *timeperiod* return pairs.
 ///
@@ -24,39 +63,84 @@ pub fn beta<'py>(
     let y = real1.as_slice()?;
     let n = x.len();
     validation::validate_equal_length(&[(n, "real0"), (y.len(), "real1")])?;
-    let mut result = vec![f64::NAN; n];
-    // Need at least timeperiod+1 bars to compute timeperiod return pairs
-    #[allow(clippy::needless_range_loop)]
-    for i in timeperiod..n {
-        // returns from bar (i - timeperiod) to bar i  => timeperiod pairs
-        let start = i - timeperiod;
-        let mut rx = vec![0.0_f64; timeperiod];
-        let mut ry = vec![0.0_f64; timeperiod];
-        for k in 0..timeperiod {
-            let prev = start + k;
-            let curr = start + k + 1;
-            rx[k] = if x[prev] != 0.0 {
-                x[curr] / x[prev] - 1.0
-            } else {
-                f64::NAN
-            };
-            ry[k] = if y[prev] != 0.0 {
-                y[curr] / y[prev] - 1.0
-            } else {
-                f64::NAN
-            };
-        }
-        let mean_x: f64 = rx.iter().sum::<f64>() / timeperiod as f64;
-        let mean_y: f64 = ry.iter().sum::<f64>() / timeperiod as f64;
-        let cov: f64 = rx
-            .iter()
-            .zip(ry.iter())
-            .map(|(&xi, &yi)| (xi - mean_x) * (yi - mean_y))
-            .sum::<f64>()
-            / timeperiod as f64;
-        let var_x: f64 =
-            rx.iter().map(|&xi| (xi - mean_x).powi(2)).sum::<f64>() / timeperiod as f64;
-        result[i] = if var_x != 0.0 { cov / var_x } else { f64::NAN };
+
+    if x.iter().any(|value| !value.is_finite()) || y.iter().any(|value| !value.is_finite()) {
+        return Ok(beta_fallback(x, y, timeperiod).into_pyarray(py));
     }
+
+    let mut result = vec![f64::NAN; n];
+    if n <= timeperiod {
+        return Ok(result.into_pyarray(py));
+    }
+
+    let rx: Vec<f64> = x
+        .windows(2)
+        .map(|window| price_return(window[1], window[0]))
+        .collect();
+    let ry: Vec<f64> = y
+        .windows(2)
+        .map(|window| price_return(window[1], window[0]))
+        .collect();
+
+    let period = timeperiod as f64;
+    let mut invalid_pairs = 0_usize;
+    let mut sum_rx = 0.0_f64;
+    let mut sum_ry = 0.0_f64;
+    let mut sum_rx2 = 0.0_f64;
+    let mut sum_rxry = 0.0_f64;
+
+    for idx in 0..timeperiod {
+        let ret_x = rx[idx];
+        let ret_y = ry[idx];
+        if ret_x.is_finite() && ret_y.is_finite() {
+            sum_rx += ret_x;
+            sum_ry += ret_y;
+            sum_rx2 += ret_x * ret_x;
+            sum_rxry += ret_x * ret_y;
+        } else {
+            invalid_pairs += 1;
+        }
+    }
+
+    for end in timeperiod..n {
+        result[end] = if invalid_pairs == 0 {
+            let denom = period * sum_rx2 - sum_rx * sum_rx;
+            if denom != 0.0 {
+                (period * sum_rxry - sum_rx * sum_ry) / denom
+            } else {
+                f64::NAN
+            }
+        } else {
+            f64::NAN
+        };
+
+        if end + 1 < n {
+            let outgoing = end - timeperiod;
+            let incoming = end;
+
+            let outgoing_x = rx[outgoing];
+            let outgoing_y = ry[outgoing];
+            if outgoing_x.is_finite() && outgoing_y.is_finite() {
+                sum_rx -= outgoing_x;
+                sum_ry -= outgoing_y;
+                sum_rx2 -= outgoing_x * outgoing_x;
+                sum_rxry -= outgoing_x * outgoing_y;
+            } else {
+                invalid_pairs -= 1;
+            }
+
+            let incoming_x = rx[incoming];
+            let incoming_y = ry[incoming];
+            if incoming_x.is_finite() && incoming_y.is_finite() {
+                sum_rx += incoming_x;
+                sum_ry += incoming_y;
+                sum_rx2 += incoming_x * incoming_x;
+                sum_rxry += incoming_x * incoming_y;
+            } else {
+                invalid_pairs += 1;
+            }
+        }
+    }
+
     Ok(result.into_pyarray(py))
 }
diff --git a/src/statistic/common.rs b/src/statistic/common.rs
index 722ec7a..d833554 100644
--- a/src/statistic/common.rs
+++ b/src/statistic/common.rs
@@ -14,3 +14,57 @@ pub(super) fn linreg(window: &[f64]) -> (f64, f64) {
     let intercept = (sum_y - slope * sum_x) / n;
     (slope, intercept)
 }
+
+pub(crate) fn rolling_linreg_apply<F>(prices: &[f64], timeperiod: usize, mut map: F) -> Vec<f64>
+where
+    F: FnMut(f64, f64) -> f64,
+{
+    let n = prices.len();
+    let mut result = vec![f64::NAN; n];
+    if timeperiod == 0 || n < timeperiod {
+        return result;
+    }
+
+    if prices.iter().any(|value| !value.is_finite()) {
+        for end in (timeperiod - 1)..n {
+            let window = &prices[(end + 1 - timeperiod)..=end];
+            let (slope, intercept) = linreg(window);
+            result[end] = map(slope, intercept);
+        }
+        return result;
+    }
+
+    let period = timeperiod as f64;
+    let last_x = (timeperiod - 1) as f64;
+    let sum_x = last_x * period / 2.0;
+    let sum_x2 = last_x * period * (2.0 * period - 1.0) / 6.0;
+    let denom = period * sum_x2 - sum_x * sum_x;
+
+    let mut sum_y = prices[..timeperiod].iter().sum::<f64>();
+    let mut sum_xy = prices[..timeperiod]
+        .iter()
+        .enumerate()
+        .map(|(idx, &value)| idx as f64 * value)
+        .sum::<f64>();
+
+    for end in (timeperiod - 1)..n {
+        let slope = if denom != 0.0 {
+            (period * sum_xy - sum_x * sum_y) / denom
+        } else {
+            0.0
+        };
+        let intercept = (sum_y - slope * sum_x) / period;
+        result[end] = map(slope, intercept);
+
+        if end + 1 < n {
+            let outgoing = prices[end + 1 - timeperiod];
+            let incoming = prices[end + 1];
+            let prev_sum_y = sum_y;
+
+            sum_y = prev_sum_y - outgoing + incoming;
+            sum_xy = sum_xy - (prev_sum_y - outgoing) + last_x * incoming;
+        }
+    }
+
+    result
+}
diff --git a/src/statistic/correl.rs b/src/statistic/correl.rs
index f10ab4a..6131f2f 100644
--- a/src/statistic/correl.rs
+++ b/src/statistic/correl.rs
@@ -2,6 +2,35 @@ use crate::validation;
 use numpy::{IntoPyArray, PyArray1, PyReadonlyArray1};
 use pyo3::prelude::*;
 
+fn correl_fallback(x: &[f64], y: &[f64], timeperiod: usize) -> Vec<f64> {
+    let n = x.len();
+    let mut result = vec![f64::NAN; n];
+    for end in (timeperiod - 1)..n {
+        let wx = &x[(end + 1 - timeperiod)..=end];
+        let wy = &y[(end + 1 - timeperiod)..=end];
+        let mean_x = wx.iter().sum::<f64>() / timeperiod as f64;
+        let mean_y = wy.iter().sum::<f64>() / timeperiod as f64;
+        let cov = wx
+            .iter()
+            .zip(wy.iter())
+            .map(|(&xi, &yi)| (xi - mean_x) * (yi - mean_y))
+            .sum::<f64>();
+        let std_x = wx
+            .iter()
+            .map(|&xi| (xi - mean_x).powi(2))
+            .sum::<f64>()
+            .sqrt();
+        let std_y = wy
+            .iter()
+            .map(|&yi| (yi - mean_y).powi(2))
+            .sum::<f64>()
+            .sqrt();
+        let denom = std_x * std_y;
+        result[end] = if denom != 0.0 { cov / denom } else { f64::NAN };
+    }
+    result
+}
+
 /// Pearson correlation coefficient between two series over the rolling window.
 #[pyfunction]
 #[pyo3(signature = (real0, real1, timeperiod = 30))]
@@ -16,21 +45,57 @@ pub fn correl<'py>(
     let y = real1.as_slice()?;
     let n = x.len();
     validation::validate_equal_length(&[(n, "real0"), (y.len(), "real1")])?;
+
+    if x.iter().any(|value| !value.is_finite()) || y.iter().any(|value| !value.is_finite()) {
+        return Ok(correl_fallback(x, y, timeperiod).into_pyarray(py));
+    }
+
     let mut result = vec![f64::NAN; n];
-    for i in (timeperiod - 1)..n {
-        let wx = &x[(i + 1 - timeperiod)..=i];
-        let wy = &y[(i + 1 - timeperiod)..=i];
-        let mean_x: f64 = wx.iter().sum::<f64>() / timeperiod as f64;
-        let mean_y: f64 = wy.iter().sum::<f64>() / timeperiod as f64;
-        let cov: f64 = wx
-            .iter()
-            .zip(wy.iter())
-            .map(|(&xi, &yi)| (xi - mean_x) * (yi - mean_y))
-            .sum::<f64>();
-        let std_x: f64 = (wx.iter().map(|&xi| (xi - mean_x).powi(2)).sum::<f64>()).sqrt();
-        let std_y: f64 = (wy.iter().map(|&yi| (yi - mean_y).powi(2)).sum::<f64>()).sqrt();
-        let denom = std_x * std_y;
-        result[i] = if denom != 0.0 { cov / denom } else { f64::NAN };
+    if n < timeperiod {
+        return Ok(result.into_pyarray(py));
+    }
+
+    let period = timeperiod as f64;
+    let mut sum_x = x[..timeperiod].iter().sum::<f64>();
+    let mut sum_y = y[..timeperiod].iter().sum::<f64>();
+    let mut sum_x2 = x[..timeperiod]
+        .iter()
+        .map(|value| value * value)
+        .sum::<f64>();
+    let mut sum_y2 = y[..timeperiod]
+        .iter()
+        .map(|value| value * value)
+        .sum::<f64>();
+    let mut sum_xy = x[..timeperiod]
+        .iter()
+        .zip(y[..timeperiod].iter())
+        .map(|(&lhs, &rhs)| lhs * rhs)
+        .sum::<f64>();
+
+    for end in (timeperiod - 1)..n {
+        let denom_x = period * sum_x2 - sum_x * sum_x;
+        let denom_y = period * sum_y2 - sum_y * sum_y;
+        result[end] = if denom_x > 0.0 && denom_y > 0.0 {
+            (period * sum_xy - sum_x * sum_y) / (denom_x * denom_y).sqrt()
+        } else {
+            f64::NAN
+        };
+
+        if end + 1 < n {
+            let outgoing = end + 1 - timeperiod;
+            let incoming = end + 1;
+
+            let outgoing_x = x[outgoing];
+            let outgoing_y = y[outgoing];
+            let incoming_x = x[incoming];
+            let incoming_y = y[incoming];
+
+            sum_x += incoming_x - outgoing_x;
+            sum_y += incoming_y - outgoing_y;
+            sum_x2 += incoming_x * incoming_x - outgoing_x * outgoing_x;
+            sum_y2 += incoming_y * incoming_y - outgoing_y * outgoing_y;
+            sum_xy += incoming_x * incoming_y - outgoing_x * outgoing_y;
+        }
     }
     Ok(result.into_pyarray(py))
 }
diff --git a/src/statistic/linearreg.rs b/src/statistic/linearreg.rs
index ca0fae3..e22e7ec 100644
--- a/src/statistic/linearreg.rs
+++ b/src/statistic/linearreg.rs
@@ -1,4 +1,4 @@
-use super::common::linreg;
+use super::common::rolling_linreg_apply;
 use crate::validation;
 use numpy::{IntoPyArray, PyArray1, PyReadonlyArray1};
 use pyo3::prelude::*;
@@ -14,13 +14,10 @@ pub fn linearreg<'py>(
 ) -> PyResult<Bound<'py, PyArray1<f64>>> {
     validation::validate_timeperiod(timeperiod, "timeperiod", 1)?;
     let prices = close.as_slice()?;
-    let n = prices.len();
-    let mut result = vec![f64::NAN; n];
-    for i in (timeperiod - 1)..n {
-        let window = &prices[(i + 1 - timeperiod)..=i];
-        let (slope, intercept) = linreg(window);
-        result[i] = intercept + slope * (timeperiod - 1) as f64;
-    }
+    let last_x = (timeperiod - 1) as f64;
+    let result = rolling_linreg_apply(prices, timeperiod, |slope, intercept| {
+        intercept + slope * last_x
+    });
     Ok(result.into_pyarray(py))
 }
 
@@ -34,13 +31,7 @@ pub fn linearreg_slope<'py>(
 ) -> PyResult<Bound<'py, PyArray1<f64>>> {
     validation::validate_timeperiod(timeperiod, "timeperiod", 1)?;
     let prices = close.as_slice()?;
-    let n = prices.len();
-    let mut result = vec![f64::NAN; n];
-    for i in (timeperiod - 1)..n {
-        let window = &prices[(i + 1 - timeperiod)..=i];
-        let (slope, _) = linreg(window);
-        result[i] = slope;
-    }
+    let result = rolling_linreg_apply(prices, timeperiod, |slope, _| slope);
     Ok(result.into_pyarray(py))
 }
 
@@ -54,13 +45,7 @@ pub fn linearreg_intercept<'py>(
 ) -> PyResult<Bound<'py, PyArray1<f64>>> {
     validation::validate_timeperiod(timeperiod, "timeperiod", 1)?;
     let prices = close.as_slice()?;
-    let n = prices.len();
-    let mut result = vec![f64::NAN; n];
-    for i in (timeperiod - 1)..n {
-        let window = &prices[(i + 1 - timeperiod)..=i];
-        let (_, intercept) = linreg(window);
-        result[i] = intercept;
-    }
+    let result = rolling_linreg_apply(prices, timeperiod, |_, intercept| intercept);
     Ok(result.into_pyarray(py))
 }
 
@@ -74,13 +59,7 @@ pub fn linearreg_angle<'py>(
 ) -> PyResult<Bound<'py, PyArray1<f64>>> {
     validation::validate_timeperiod(timeperiod, "timeperiod", 1)?;
     let prices = close.as_slice()?;
-    let n = prices.len();
-    let mut result = vec![f64::NAN; n];
-    for i in (timeperiod - 1)..n {
-        let window = &prices[(i + 1 - timeperiod)..=i];
-        let (slope, _) = linreg(window);
-        result[i] = slope.atan() * 180.0 / PI;
-    }
+    let result = rolling_linreg_apply(prices, timeperiod, |slope, _| slope.atan() * 180.0 / PI);
     Ok(result.into_pyarray(py))
 }
 
@@ -94,13 +73,9 @@ pub fn tsf<'py>(
 ) -> PyResult<Bound<'py, PyArray1<f64>>> {
     validation::validate_timeperiod(timeperiod, "timeperiod", 1)?;
     let prices = close.as_slice()?;
-    let n = prices.len();
-    let mut result = vec![f64::NAN; n];
-    for i in (timeperiod - 1)..n {
-        let window = &prices[(i + 1 - timeperiod)..=i];
-        let (slope, intercept) = linreg(window);
-        // Forecast one period ahead of the last point in the window
-        result[i] = intercept + slope * timeperiod as f64;
-    }
+    let forecast_x = timeperiod as f64;
+    let result = rolling_linreg_apply(prices, timeperiod, |slope, intercept| {
+        intercept + slope * forecast_x
+    });
     Ok(result.into_pyarray(py))
 }
diff --git a/src/statistic/mod.rs b/src/statistic/mod.rs
index d4a0691..7fb2fb4 100644
--- a/src/statistic/mod.rs
+++ b/src/statistic/mod.rs
@@ -2,7 +2,7 @@
 //! Each function (or closely related group) lives in its own file.
 
 mod beta;
-mod common;
+pub(crate) mod common;
 mod correl;
 mod linearreg;
 mod stddev;
diff --git a/tests/unit/indicators/test_statistic.py b/tests/unit/indicators/test_statistic.py
index 0d537b3..96e2953 100644
--- a/tests/unit/indicators/test_statistic.py
+++ b/tests/unit/indicators/test_statistic.py
@@ -27,6 +27,68 @@ LINDATA = np.arange(1.0, 6.0)  # [1,2,3,4,5]
 CONSTDATA = np.ones(10)  # all 1.0
 
 
+def _naive_linreg_window(window: np.ndarray) -> tuple[float, float]:
+    x = np.arange(len(window), dtype=np.float64)
+    sum_x = float(np.sum(x))
+    sum_y = float(np.sum(window))
+    sum_xy = float(np.sum(x * window))
+    sum_x2 = float(np.sum(x * x))
+    n = float(len(window))
+    denom = n * sum_x2 - sum_x * sum_x
+    slope = (n * sum_xy - sum_x * sum_y) / denom if denom != 0.0 else 0.0
+    intercept = (sum_y - slope * sum_x) / n
+    return slope, intercept
+
+
+def _naive_linearreg(series: np.ndarray, timeperiod: int, x_value: float) -> np.ndarray:
+    out = np.full(len(series), np.nan, dtype=np.float64)
+    for end in range(timeperiod - 1, len(series)):
+        slope, intercept = _naive_linreg_window(series[end + 1 - timeperiod : end + 1])
+        out[end] = intercept + slope * x_value
+    return out
+
+
+def _naive_correl(x: np.ndarray, y: np.ndarray, timeperiod: int) -> np.ndarray:
+    out = np.full(len(x), np.nan, dtype=np.float64)
+    for end in range(timeperiod - 1, len(x)):
+        x_window = x[end + 1 - timeperiod : end + 1]
+        y_window = y[end + 1 - timeperiod : end + 1]
+        mean_x = float(np.sum(x_window)) / timeperiod
+        mean_y = float(np.sum(y_window)) / timeperiod
+        cov = float(np.sum((x_window - mean_x) * (y_window - mean_y)))
+        std_x = float(np.sqrt(np.sum((x_window - mean_x) ** 2)))
+        std_y = float(np.sqrt(np.sum((y_window - mean_y) ** 2)))
+        denom = std_x * std_y
+        out[end] = cov / denom if denom != 0.0 else np.nan
+    return out
+
+
+def _naive_beta(x: np.ndarray, y: np.ndarray, timeperiod: int) -> np.ndarray:
+    out = np.full(len(x), np.nan, dtype=np.float64)
+    for end in range(timeperiod, len(x)):
+        start = end - timeperiod
+        rx = np.array(
+            [
+                x[idx + 1] / x[idx] - 1.0 if x[idx] != 0.0 else np.nan
+                for idx in range(start, end)
+            ],
+            dtype=np.float64,
+        )
+        ry = np.array(
+            [
+                y[idx + 1] / y[idx] - 1.0 if y[idx] != 0.0 else np.nan
+                for idx in range(start, end)
+            ],
+            dtype=np.float64,
+        )
+        mean_x = float(np.sum(rx)) / timeperiod
+        mean_y = float(np.sum(ry)) / timeperiod
+        cov = float(np.sum((rx - mean_x) * (ry - mean_y))) / timeperiod
+        var_x = float(np.sum((rx - mean_x) ** 2)) / timeperiod
+        out[end] = cov / var_x if var_x != 0.0 else np.nan
+    return out
+
+
 # ---------------------------------------------------------------------------
 # STDDEV
 # ---------------------------------------------------------------------------
@@ -100,6 +162,11 @@ class TestLINEARREG:
     def test_length(self):
         assert len(LINEARREG(_A, 14)) == N
 
+    def test_matches_naive_regression(self):
+        expected = _naive_linearreg(_A, timeperiod=14, x_value=13.0)
+        result = LINEARREG(_A, timeperiod=14)
+        np.testing.assert_allclose(result, expected, equal_nan=True)
+
 
 # ---------------------------------------------------------------------------
 # LINEARREG_SLOPE
@@ -179,6 +246,11 @@ class TestBETA:
         valid = result[~np.isnan(result)]
         assert np.all(np.isfinite(valid))
 
+    def test_matches_naive_beta(self):
+        expected = _naive_beta(_A, _B, timeperiod=5)
+        result = BETA(_A, _B, timeperiod=5)
+        np.testing.assert_allclose(result, expected, equal_nan=True)
+
 
 # ---------------------------------------------------------------------------
 # CORREL
@@ -205,6 +277,11 @@ class TestCOREL:
     def test_length(self):
         assert len(CORREL(_A, _B, 10)) == N
 
+    def test_matches_naive_correlation(self):
+        expected = _naive_correl(_A, _B, timeperiod=10)
+        result = CORREL(_A, _B, timeperiod=10)
+        np.testing.assert_allclose(result, expected, equal_nan=True)
+
 
 # ---------------------------------------------------------------------------
 # TSF
@@ -226,3 +303,8 @@ class TestTSF:
 
     def test_length(self):
         assert len(TSF(_A, 14)) == N
+
+    def test_matches_naive_tsf(self):
+        expected = _naive_linearreg(_A, timeperiod=14, x_value=14.0)
+        result = TSF(_A, timeperiod=14)
+        np.testing.assert_allclose(result, expected, equal_nan=True)
diff --git a/tests/unit/test_data_pipeline.py b/tests/unit/test_data_pipeline.py
index ca9cfc0..ec2876f 100644
--- a/tests/unit/test_data_pipeline.py
+++ b/tests/unit/test_data_pipeline.py
@@ -323,6 +323,21 @@ class TestSignalComposition:
         score = compose(sigs, method="rank")
         assert score.shape == (30,)
 
+    def test_compose_rank_matches_manual_column_ranks(self):
+        from ferro_ta.analysis.signals import compose
+
+        sigs = np.array(
+            [
+                [3.0, 1.0],
+                [1.0, 2.0],
+                [2.0, 2.0],
+            ],
+            dtype=np.float64,
+        )
+        score = compose(sigs, method="rank")
+        expected = np.array([4.0, 3.5, 4.5], dtype=np.float64)
+        np.testing.assert_allclose(score, expected)
+
     def test_compose_equal_weights_default(self):
         from ferro_ta.analysis.signals import compose
 
@@ -574,6 +589,75 @@ class TestFeatureMatrix:
         fm = feature_matrix(ohlcv, ["SMA"])
         assert "SMA" in fm
 
+    def test_feature_matrix_mixed_fastpath_and_multi_output(self):
+        from ferro_ta.analysis.features import feature_matrix
+
+        o, h, l, c, v = _make_ohlcv(80)
+        ohlcv = {"close": c, "high": h, "low": l, "open": o, "volume": v}
+        fm = feature_matrix(
+            ohlcv,
+            [
+                ("SMA", {"timeperiod": 10}),
+                ("ATR", {"timeperiod": 14}),
+                ("BBANDS", {"timeperiod": 10}, 1),
+            ],
+        )
+        assert "SMA" in fm
+        assert "ATR" in fm
+        assert "BBANDS_1" in fm
+
+
+class TestComputeMany:
+    def test_close_indicators_match_public_api(self):
+        from ferro_ta import EMA, RSI, SMA
+        from ferro_ta.data.batch import compute_many
+
+        _, _, _, close, _ = _make_ohlcv(80)
+        results = compute_many(
+            [
+                ("SMA", {"timeperiod": 10}),
+                ("EMA", {"timeperiod": 12}),
+                ("RSI", {"timeperiod": 14}),
+            ],
+            close=close,
+        )
+
+        np.testing.assert_allclose(results[0], SMA(close, timeperiod=10), equal_nan=True)
+        np.testing.assert_allclose(results[1], EMA(close, timeperiod=12), equal_nan=True)
+        np.testing.assert_allclose(results[2], RSI(close, timeperiod=14), equal_nan=True)
+
+    def test_hlc_indicators_match_public_api(self):
+        from ferro_ta import ADX, ATR
+        from ferro_ta.data.batch import compute_many
+
+        _, high, low, close, _ = _make_ohlcv(80)
+        results = compute_many(
+            [
+                ("ATR", {"timeperiod": 14}),
+                ("ADX", {"timeperiod": 14}),
+            ],
+            close=close,
+            high=high,
+            low=low,
+        )
+
+        np.testing.assert_allclose(
+            results[0], ATR(high, low, close, timeperiod=14), equal_nan=True
+        )
+        np.testing.assert_allclose(
+            results[1], ADX(high, low, close, timeperiod=14), equal_nan=True
+        )
+
+    def test_unsupported_kwargs_fall_back_cleanly(self):
+        from ferro_ta import STDDEV
+        from ferro_ta.data.batch import compute_many
+
+        _, _, _, close, _ = _make_ohlcv(80)
+        result = compute_many([("STDDEV", {"timeperiod": 10, "nbdev": 2.0})], close=close)
+        np.testing.assert_allclose(
+            result[0], STDDEV(close, timeperiod=10, nbdev=2.0), equal_nan=True
+        )
+
 
 # ---------------------------------------------------------------------------
 # Viz (smoke tests)
diff --git a/tests/unit/test_infrastructure.py b/tests/unit/test_infrastructure.py
index f45c6f9..a68ed65 100644
--- a/tests/unit/test_infrastructure.py
+++ b/tests/unit/test_infrastructure.py
@@ -297,6 +297,40 @@ class TestBacktest:
             result_no_slip.n_trades == 0
         )
 
+    def test_commission_matches_reference_loop(self):
+        close = np.array([100.0, 102.0, 101.0, 104.0, 103.0, 105.0], dtype=np.float64)
+        raw_signals = np.array([0.0, 1.0, 1.0, -1.0, -1.0, 0.0], dtype=np.float64)
+
+        def strategy(_, **__):
+            return raw_signals
+
+        commission = 0.02
+        result = backtest(close, strategy=strategy, commission_per_trade=commission)
+
+        expected_positions = np.array(
+            [0.0, 0.0, 1.0, 1.0, -1.0, -1.0], dtype=np.float64
+        )
+        expected_returns = np.empty_like(close)
+        expected_returns[0] = 0.0
+        expected_returns[1:] = np.diff(close) / close[:-1]
+        expected_strategy_returns = expected_positions * expected_returns
+        position_changed = np.concatenate(
+            [[False], expected_positions[1:] != expected_positions[:-1]]
+        )
+
+        expected_equity = np.empty_like(close)
+        expected_equity[0] = 1.0
+        for i in range(1, len(close)):
+            expected_equity[i] = expected_equity[i - 1] * (
+                1.0 + expected_strategy_returns[i]
+            )
+            if position_changed[i]:
+                expected_equity[i] -= commission
+
+        np.testing.assert_allclose(result.positions, expected_positions)
+        np.testing.assert_allclose(result.strategy_returns, expected_strategy_returns)
+        np.testing.assert_allclose(result.equity, expected_equity)
+
 
 # ---------------------------------------------------------------------------
 # Plugin / Registry
@@ -509,7 +543,13 @@ class TestChoppinessIndex:
 # ---------------------------------------------------------------------------
 
 from ferro_ta import EMA, RSI, SMA
-from ferro_ta.data.batch import batch_apply, batch_ema, batch_rsi, batch_sma
+from ferro_ta.data.batch import (
+    batch_apply,
+    batch_atr,
+    batch_ema,
+    batch_rsi,
+    batch_sma,
+)
 
 
 class TestBatchSMA:
@@ -587,6 +627,15 @@ class TestBatchApply:
             batch_apply(np.zeros((5, 5, 5)), SMA, timeperiod=3)
 
 
+class TestBatchShapeValidation:
+    def test_batch_atr_shape_mismatch_raises(self):
+        high = np.ones((5, 2), dtype=np.float64)
+        low = np.ones((4, 2), dtype=np.float64)
+        close = np.ones((5, 2), dtype=np.float64)
+        with pytest.raises(ValueError, match="shape"):
+            batch_atr(high, low, close, timeperiod=3)
+
+
 # ---------------------------------------------------------------------------
 # Release playbook and version consistency
 # ---------------------------------------------------------------------------
diff --git a/uv.lock b/uv.lock
index 9320caa..71931aa 100644
--- a/uv.lock
+++ b/uv.lock
@@ -609,7 +609,7 @@ wheels = [
 
 [[package]]
 name = "ferro-ta"
-version = "1.0.0"
+version = "1.0.2"
 source = { editable = "." }
 dependencies = [
     { name = "numpy" },
diff --git a/wasm/Cargo.toml b/wasm/Cargo.toml
index b7472e8..dc524aa 100644
--- a/wasm/Cargo.toml
+++ b/wasm/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "ferro_ta_wasm"
-version = "1.0.0"
+version = "1.0.2"
 edition = "2021"
 description = "WebAssembly bindings for ferro-ta technical analysis indicators"
 license = "MIT"
diff --git a/wasm/bench.js b/wasm/bench.js
new file mode 100644
index 0000000..1f52b48
--- /dev/null
+++ b/wasm/bench.js
@@ -0,0 +1,113 @@
+const fs = require("node:fs");
+const path = require("node:path");
+const { performance } = require("node:perf_hooks");
+
+const wasm = require("./pkg/ferro_ta_wasm.js");
+
+function parseArgs(argv) {
+  const args = { bars: 100000, json: null };
+  for (let idx = 0; idx < argv.length; idx += 1) {
+    const token = argv[idx];
+    if (token === "--bars") {
+      args.bars = Number(argv[idx + 1]);
+      idx += 1;
+    } else if (token === "--json") {
+      args.json = argv[idx + 1];
+      idx += 1;
+    }
+  }
+  return args;
+}
+
+function makeSeries(length) {
+  const close = new Float64Array(length);
+  const high = new Float64Array(length);
+  const low = new Float64Array(length);
+  let value = 100.0;
+  for (let idx = 0; idx < length; idx += 1) {
+    value += Math.sin(idx / 13.0) * 0.35 + Math.cos(idx / 29.0) * 0.18;
+    close[idx] = value;
+    high[idx] = value + 1.25;
+    low[idx] = value - 1.10;
+  }
+  return { close, high, low };
+}
+
+function timeMin(fn, rounds = 7) {
+  fn();
+  let best = Number.POSITIVE_INFINITY;
+  for (let round = 0; round < rounds; round += 1) {
+    const started = performance.now();
+    fn();
+    best = Math.min(best, performance.now() - started);
+  }
+  return best;
+}
+
+function runBenchmark({ bars }) {
+  const { close, high, low } = makeSeries(bars);
+  const cases = [
+    ["SMA", () => wasm.sma(close, 20)],
+    ["EMA", () => wasm.ema(close, 20)],
+    ["RSI", () => wasm.rsi(close, 14)],
+    ["ATR", () => wasm.atr(high, low, close, 14)],
+    ["BBANDS", () => wasm.bbands(close, 20, 2.0, 2.0)],
+  ];
+
+  const results = cases.map(([name, fn]) => {
+    const elapsedMs = timeMin(fn);
+    return {
+      indicator: name,
+      elapsed_ms: Number(elapsedMs.toFixed(4)),
+      ns_per_bar: Number(((elapsedMs * 1e6) / bars).toFixed(2)),
+      million_bars_per_second: Number((((bars / 1e6) / (elapsedMs / 1000))).toFixed(2)),
+    };
+  });
+
+  return {
+    metadata: {
+      suite: "wasm",
+      runtime: {
+        generated_at_utc: new Date().toISOString(),
+        node_version: process.version,
+        platform: process.platform,
+        arch: process.arch,
+      },
+      dataset: {
+        bars,
+      },
+    },
+    results,
+  };
+}
+
+function printResults(payload) {
+  const bars = payload.metadata.dataset.bars;
+  console.log(`WASM Benchmark: ${bars} bars`);
+  console.log("----------------------------------------------------------------");
+  console.log(
+    `${"Indicator".padEnd(12)}${"Elapsed (ms)".padStart(14)}${"ns/bar".padStart(12)}${"M bars/s".padStart(12)}`
+  );
+  console.log("----------------------------------------------------------------");
+  for (const row of payload.results) {
+    console.log(
+      `${row.indicator.padEnd(12)}${row.elapsed_ms.toFixed(2).padStart(14)}${row.ns_per_bar
+        .toFixed(2)
+        .padStart(12)}${row.million_bars_per_second.toFixed(2).padStart(12)}`
+    );
+  }
+}
+
+function main() {
+  const args = parseArgs(process.argv.slice(2));
+  const payload = runBenchmark(args);
+  printResults(payload);
+
+  if (args.json) {
+    const outputPath = path.resolve(args.json);
+    fs.writeFileSync(outputPath, `${JSON.stringify(payload, null, 2)}\n`, "utf8");
+    console.log(`\nWrote JSON results to ${outputPath}`);
+  }
+}
+
+main();
diff --git a/wasm/package.json b/wasm/package.json
index 6c4f0ce..f6a78ca 100644
--- a/wasm/package.json
+++ b/wasm/package.json
@@ -1,12 +1,13 @@
 {
   "name": "ferro-ta-wasm",
-  "version": "1.0.1",
+  "version": "1.0.2",
   "description": "WebAssembly bindings for ferro-ta technical analysis indicators",
   "main": "pkg/ferro_ta_wasm.js",
   "types": "pkg/ferro_ta_wasm.d.ts",
   "files": ["pkg"],
   "scripts": {
     "build": "wasm-pack build --target nodejs --out-dir pkg",
+    "bench": "node bench.js",
     "prepack": "npm run build && node -e \"require('fs').rmSync('pkg/.gitignore', { force: true })\"",
     "test": "wasm-pack test --node"
   },