Installation

npm install volume-anomaly

Quickstart (no math required)

Feed 15–30 minutes of trades (oldest first) into one call and read the answer in plain language — scan() splits baseline/recent by time itself, so there is nothing to slice and nothing to get wrong:

import { scan, explain } from 'volume-anomaly';

const result = scan(trades);       // trades: last 15–30+ min, oldest first

if (result.anomaly) {
  console.log(explain(result));
  // Volume anomaly detected (severity: extreme) — confidence 0.97 vs alert threshold 0.75.
  // Strongest signal: volume ran ~81 robust sigma above the recent typical level at the 5s scale.
  // Peak at 2025-03-01T07:06:10.323Z.
  // Order flow at the peak: 74% sell-side.
  // Follow-through ranking (moveScore): 0.92 — top-tier; historically precedes real price movement...
}

switch (result.severity) {        // 'none' | 'notable' | 'strong' | 'extreme'
  case 'extreme': /* page someone */ break;
  case 'strong':  /* alert */        break;
}

Runnable on real data committed to this repo: node examples/quickstart.mjs.

Field glossary — what the numbers mean in plain words:

| Field | Plain meaning | |-------|---------------| | confidence | How unusual this moment is, 0–1. Alert when ≥ your threshold (default 0.75) | | severity | The same as an enum: none / notable / strong / extreme | | moveScore | How strongly moments like this have historically preceded actual price movement. Use to rank/prioritize alerts, not as the alert itself | | direction | Which side drove the burst (long = buyers). Describes the event; measured on real data it does not predict the next move's direction | | burstImbalance | Buy/sell balance inside the burst itself, −1 (all sells) … +1 (all buys) | | peakTs | When the anomaly peaked (Unix ms) | | stats.zVol / zRate | How many "sigmas" volume / trade-arrival rate ran above what was typical recently, per time scale | | hawkesLambda | Instantaneous trade-arrival intensity under the fitted model | | runLength / cusumStat | How long since the order-flow regime last changed / accumulated drift evidence |

detector.calibrationReport tells you (in words) whether the detector could adapt to your data or is running on universal defaults, and what to feed it to fix that.

Overview

The library detects abnormal surges in trade flow — sudden acceleration of arrivals and volume waves — from a raw stream of aggregated trades. The direction of the trade must come from your own analysis (fundamental, technical). This library answers a narrower question: is right now a statistically unusual moment in market microstructure?

The primary statistic is a self-calibrated robust z-score of trade rate and volume rate: the detector measures "trades per 5 s / per 30 s" and "quantity per 5 s / per 30 s" in the detection window and scores the peaks against the median/MAD of the same statistic over the training window — how many robust σ above the recent typical level. Thresholds are calibrated on a full day of real BTCUSDT aggTrades (1.49 M trades, see test/eval.test.ts): at the default confidence = 0.75 the detector catches ≈ 92 % of locally-strong volume anomalies (robust z ≥ 8 vs the trailing hour) and ≈ 95 % of anomaly events, with ≈ 2.5 % false alarms on normal 30-second windows.

Two secondary detectors (CUSUM and BOCPD on order-flow imbalance) are computed and reported in scores/signals, but by default carry zero weight in the combined confidence: on real data they track flow-regime shifts — a different phenomenon — and any additive weight on them measurably dilutes recall. Re-weight via scoreWeights if your use case targets flow shifts specifically.

API

detect(historical, recent, confidence?)

One-shot convenience function. Trains on historical data, then evaluates the recent window. Returns a DetectionResult.

import { detect } from 'volume-anomaly';
import type { IAggregatedTradeData } from 'volume-anomaly';

// historical: 15–30 MINUTES of trades (time span matters more than count —
// the baseline must cover enough market time to define "typical")
const historical: IAggregatedTradeData[] = await getAggregatedTrades('BTCUSDT', 10_000);
const recent:     IAggregatedTradeData[] = await getAggregatedTrades('BTCUSDT', 300);

const result = detect(historical, recent, 0.75);
// {
//   anomaly:      true,
//   confidence:   0.98,          // combined score (volume channel by default)
//   scores:       { hawkes: 0.98, cusum: 0.61, bocpd: 0.00 },  // raw sub-scores
//   stats:        { zRate: 22.3, zVol: 104.1, zRateSlow: 25.3, zVolSlow: 87.2, lambdaRatio: 3.1 },
//   imbalance:    0.72,          // buy-side dominance
//   hawkesLambda: 61.8,          // peak intensity (trades/sec)
//   cusumStat:    3.1,           // CUSUM accumulator
//   runLength:    2,             // periods since last changepoint
//   signals: [
//     { kind: 'volume_spike',    score: 0.98, meta: { lambda: 61.8, zRate: 22.3, zVol: 104.1, ... } },
//     { kind: 'imbalance_shift', score: 0.72, meta: { imbalance: 0.72, absImbalance: 0.72 } },
//   ]
// }

Parameters:

| Parameter | Type | Default | Description | |-----------|------|---------|-------------| | historical | IAggregatedTradeData[] | required | Baseline window for training. Must span 15–30 minutes of market time (code minimum: 50 trades). A count-based baseline ("last 500 trades") adapts its duration to market pace and masks the very burst being detected | | recent | IAggregatedTradeData[] | required | Window to evaluate. Typically 100–500 trades; must span at least rateHorizonSec (5 s) for the fast channel to engage | | confidence | number | 0.75 | Threshold in (0, 1). result.anomaly = result.confidence >= confidence |

Returns: DetectionResult

interface DetectionResult {
  anomaly:      boolean;       // confidence >= threshold
  confidence:   number;        // composite score [0,1]
  scores:       { hawkes: number; cusum: number; bocpd: number }; // raw sub-scores
  stats:        { zRate: number; zVol: number;                    // robust z of peak
                  zRateSlow: number; zVolSlow: number;            // rolling rates
                  lambdaRatio: number };                          // peak λ vs training
  signals:      AnomalySignal[]; // which sub-detectors fired
  imbalance:    number;        // buy/sell balance [-1, +1]
  hawkesLambda: number;        // conditional intensity at last trade (trades/sec)
  cusumStat:    number;        // max(S⁺, S⁻) — CUSUM accumulator
  runLength:    number;        // MAP run length — periods since last changepoint
}

predict(historical, recent, confidence?, imbalanceThreshold?)

One-shot convenience function. Wraps detect() and adds a directional signal derived from the burst-local order flow.

import { predict } from 'volume-anomaly';

const result = predict(historical, recent, 0.75, 0.3);
// {
//   anomaly:        true,
//   confidence:     0.81,
//   direction:      'long',    // 'long' | 'short' | 'neutral'
//   imbalance:      0.42,      // full-window flow balance
//   burstImbalance: 0.72,      // flow balance INSIDE the peak burst window
//   moveScore:      0.68,      // predictive ranking score for forward price response
// }

Parameters:

| Parameter | Type | Default | Description | |-----------|------|---------|-------------| | historical | IAggregatedTradeData[] | required | Baseline window for training (≥ 50 trades) | | recent | IAggregatedTradeData[] | required | Window to evaluate | | confidence | number | 0.75 | Anomaly threshold [0,1] | | imbalanceThreshold | number | (trained) | Override the trained directional threshold. Omit to use the value derived automatically from training data (p75 of the rolling signed imbalance series) |

Direction logic:

thr = max(0, imbalanceThreshold)  (if provided explicitly)
    = max(0, detector.trainedModels.imbalanceThreshold)  (otherwise — p75 from training)

direction = 'long'    if anomaly && burstImbalance >  +thr
direction = 'short'   if anomaly && burstImbalance < −thr
direction = 'neutral' otherwise (no anomaly, or balanced flow at the burst)

Direction reads the burst-local imbalance — the order flow inside the rolling window that actually produced the anomaly score — not the full-window average, which dilutes a burst's onset with surrounding two-way flow (measured on real data: a pure sell burst reads −0.9 at the burst but only −0.42 over the full bucket). The burst imbalance is additionally shrunk toward neutral by effective sample size (Kish n_eff over qty), so a three-trade or single-whale "burst" cannot fake directional conviction.

On a neutral/balanced market thr will be near 0 (most windows have close-to-zero imbalance, p75 ≈ 0.1–0.2). On a trending market the p75 shifts upward with the trend, so the bar for direction=long rises accordingly — preventing chronic false long signals during a bull run where sustained buy imbalance is normal, not anomalous. The threshold is clamped at zero, so 'long' always implies buy-side burst flow and 'short' sell-side.

⚠️ direction is descriptive, not predictive. It identifies which side drove the burst — essential for interpreting the event — but measured on the full-day benchmark, no flow-direction statistic (hard-window imbalance, exponentially-weighted flow, or a full bivariate buy/sell Hawkes excitation model — all evaluated) predicted the sign of the forward 1–5 minute return better than chance (~50% hit rate): the price impact of aggressive flow is contemporaneous, and by detection time it is already in the price. Use direction to understand the event and moveScore to rank expected follow-through magnitude; do not trade the direction sign blindly. The bivariate Hawkes model remains available for research via the math exports (biHawkesFit, biExcitationShare).

Returns: PredictionResult

interface PredictionResult {
  anomaly:        boolean;    // confidence >= threshold
  confidence:     number;     // composite score [0,1]
  direction:      Direction;  // 'long' | 'short' | 'neutral'
  imbalance:      number;     // full-window buy/sell balance [-1, +1]
  burstImbalance: number;     // buy/sell balance inside the peak burst window
  moveScore:      number;     // predictive ranking score [0,1) — see below
}

moveScore — ranking alerts by expected follow-through. confidence answers "is this an anomaly?" and is optimized for the detection frontier; for ranking which alerts precede actual price movement it is handicapped by design (per-baseline adaptive levels re-zero the scale each training window, and fast channels are ranking noise). moveScore is the peak long-scale volume z through a fixed universal mapping — no baseline adaptation, comparable across windows and across time. Measured on the full-day benchmark: ranks the forward 1-minute price range at AUC ≈ 0.64 vs ≈ 0.60 for confidence (naive momentum ≈ 0.57). Use it to prioritize or size — never as the detection threshold.

Practical usage with getAggregatedTrades from backtest-kit:

⚠️ Never pass the same trades in both historical and recent. Training calibrates the baseline. If recent overlaps with historical, any anomaly in that period is absorbed into the baseline — the detector learns to treat it as normal and misses it. Always slice so that recent starts where historical ends.

import { predict } from 'volume-anomaly';
import type { IAggregatedTradeData } from 'volume-anomaly';

// Your data-fetching function — returns the last `limit` trades, oldest first:
declare function getAggregatedTrades(
  symbol: string,
  limit:  number,
): Promise<IAggregatedTradeData[]>;

// ── One-shot (single API call, zero overlap) ──────────────────────────────────

const N_train  = 1200;  // calibration window
const N_detect = 200;   // window to evaluate

const all        = await getAggregatedTrades('BTCUSDT', N_train + N_detect);
const historical = all.slice(0, N_train);  // older 1200 trades — baseline
const recent     = all.slice(N_train);     // newest 200 trades — no overlap

const result = predict(historical, recent, 0.75);
// {
//   anomaly:    true,
//   confidence: 0.83,
//   direction:  'long',   // 'long' | 'short' | 'neutral'
//   imbalance:  0.61,
// }

if (result.anomaly) {
  console.log(`direction=${result.direction}  confidence=${result.confidence.toFixed(2)}`);
}

predict() trains a fresh detector on every call. For continuous monitoring (many detect() calls from one trained model) use VolumeAnomalyDetector directly — see the class API below.

scan(trades, options?) / explain(result, threshold?)

The no-slicing convenience layer (see Quickstart). scan(trades, { recentSec?, confidence?, ...DetectorConfig }) sorts the stream, evaluates the last recentSec seconds (default 30) against a baseline trained on everything before — overlap is impossible by construction; throws a clear error when fewer than 50 baseline trades remain. Returns DetectionResult & { direction }. explain() renders any result (detect/scan/predict) as plain-language lines: verdict + severity, the strongest channel and scale, peak time, order flow at the peak, and how to read moveScore/direction.

new VolumeAnomalyDetector(config?)

Stateful class. Use when you need to re-use fitted models across multiple detect() calls without re-training, or when you want to tune individual model parameters.

import { VolumeAnomalyDetector } from 'volume-anomaly';

const detector = new VolumeAnomalyDetector({
  windowSize:     50,        // trades per imbalance window (CUSUM/BOCPD)
  rateHorizonSec: 5,         // fast rate/volume horizon (seconds)
  slowHorizonSec: 30,        // slow rate/volume horizon (seconds)
  hazardLambda:   200,       // expected periods between changepoints
  cusumKSigmas:   0.5,       // CUSUM slack k = 0.5 · σ
  cusumHSigmas:   5,         // CUSUM alarm h floor = 5 · σ
  scoreWeights:   [1, 0, 0], // [volume channel, CUSUM, BOCPD]
});

detector.train(historicalTrades);
const result = detector.detect(recentTrades, 0.75);

// Inspect fitted parameters for debugging:
const { hawkesParams, cusumParams, bocpdPrior } = detector.trainedModels!;

Config parameters:

| Parameter | Type | Default | Description | |-----------|------|---------|-------------| | windowSize | number | 50 | Number of trades per rolling imbalance window (CUSUM/BOCPD). Smaller = more reactive to local shifts, larger = smoother signal | | rateHorizonSec | number | auto (≥ 5) | Fast time horizon for the volume channel: "trades per horizon" / "qty per horizon". Catches brief intense bursts. When omitted, chosen on the fly: max(5 s, 25 × median inter-trade gap), capped by the training span — sparse instruments automatically get a horizon that carries a meaningful rate. Pass a value to pin it | | slowHorizonSec | number | auto (≥ 30) | Slow horizon for the same statistic. Catches sustained volume waves spread too evenly to concentrate at the fast horizon. Auto: 6 × fast horizon, floored at 30 s, capped by the training span | | hazardLambda | number | 200 | Expected number of windows between changepoints (BOCPD hazard rate H = 1/λ). Set lower for more frequent regime changes | | cusumKSigmas | number | 0.5 | CUSUM allowable slack k in σ units. Controls sensitivity: lower k = faster response but more false positives | | cusumHSigmas | number | 5 | CUSUM alarm threshold floor in σ units; the effective h is self-calibrated as max(hSigmas·σ, 2 × max training excursion) because the rolling series is heavily autocorrelated | | scoreWeights | [n, n, n] | [1, 0, 0] | Weights for [volume channel, CUSUM, BOCPD]. Must sum to 1. The default gives the flow-shift detectors zero weight — on real data they dilute volume-anomaly recall (see test/eval.test.ts); their scores remain available in result.scores | | imbalancePercentile | number | 75 | Percentile of the training rolling signed imbalance used as the directional threshold. p75 = direction fires only when imbalance exceeds the 75th percentile of the training distribution. Clamped at zero before use, so 'long' always implies imbalance > 0 and 'short' implies imbalance < 0 |

Every numeric option is validated at construction — a bad value (e.g. hazardLambda: 0.5, which would silently collapse the BOCPD state, or a non-integer windowSize) throws immediately with the option name and received value instead of producing a miscalibrated detector. Two input mistakes that never crash but silently rescale every time horizon 1000× are also rejected: timestamps that look like Unix seconds or microseconds (they must be milliseconds), and a percent-style confidence (75 instead of 0.75).

toJSON() / VolumeAnomalyDetector.fromJSON(snapshot)

Trained models are plain finite numbers, so a detector serializes losslessly: train once (e.g. in a worker on a schedule), ship the snapshot anywhere, and detect without re-training.

const detector = new VolumeAnomalyDetector();
detector.train(historicalTrades);

const snapshot = JSON.stringify(detector);          // uses toJSON() automatically

// … elsewhere (another process, a file, Redis, …):
const restored = VolumeAnomalyDetector.fromJSON(snapshot); // string or parsed object
const result   = restored.detect(recentTrades, 0.75);      // identical to the original

fromJSON validates the snapshot structurally — config re-passes constructor validation, and every numeric leaf of the models must be finite (JSON silently turns NaN/Infinity into null, so a corrupted snapshot fails loudly here, with the offending path in the message, rather than as NaN confidence inside detect()). Auto-chosen horizons stay auto after restore: a restored detector re-derives them if you call train() again.

confidence — how it works (critical difference from garch)

This is fundamentally different from garch's confidence parameter.

In garch, confidence is a two-sided probability passed through probit to get a z-score, which then scales a log-normal price corridor: confidence → z = probit((1+confidence)/2) → P·exp(±z·σ). The confidence there controls band width, not a classification threshold.

Here, confidence is a hard threshold on a composite score. The formula is:

score_final = w_H · score_volume + w_C · score_cusum + w_B · score_bocpd
            = score_volume                                    (default weights [1,0,0])
anomaly     = score_final >= confidence

The volume-channel score is built in three steps:

1. Multi-scale scan. The peak rolling arrival rate and volume rate are measured at every horizon in the trained family (fast ≈ 5 s, a geometric-mean mid scale, slow ≈ 30 s, and an extended scale up to 3× slow when the training span allows) and standardized per channel: t = (z − c − u·spanShift) / u, where z is the robust z of the peak against the training median/MAD.
2. Cross-type corroboration (Stouffer). The best arrival-rate t and the best volume t combine as max(t_lead, (t_lead + max(t_other, 0)) / √2): a moderate rate excess AND a moderate volume excess together outrank either alone. Scales of the same type never corroborate each other — they see the same wiggle at different bandwidths.
3. Rational sigmoid. score = 0.5 + 0.5·t/(1+|t|): 0.5 at the calibrated level, 0.75 one tail unit beyond it, approaching 1 harmonically — the top of the scale keeps ranking extreme events instead of saturating into ties, which matters when the score is used predictively.

The level c is chosen on the fly at train() time: the P85 of the peak-z values the baseline itself produced in alert-sized stretches, made robust to events inside the training window by a Gumbel fit to the LEFT quantiles of those maxima (contamination inflates only the upper quantiles; the level is min(empirical P85, Gumbel-extrapolated P85)), and floored by per-channel universal levels found by dense brute-force search on a full day of real BTCUSDT trades (test/eval.test.ts): 14 for arrival-rate channels, 6.5 for volume channels (volume z separates anomalies much more cleanly, so it earns a lower bar). The tail unit u = 5.5 is universal in z-space. On a typical baseline confidence = 0.75 fires at volume z ≥ 12 or rate z ≥ 19.5; a hotter/noisier baseline automatically raises the bar.

Practical guidance (measured on the real-data benchmark, 30 s buckets):

| confidence | Fires at (t units beyond level) | Recall (strong anomalies) | False alarms (normal buckets) | |-------------|----------------------------------|---------------------------|-------------------------------| | 0.5 | t ≥ 0 (vol z ≥ 6.5 / rate z ≥ 14) | higher | higher | | 0.75 (default) | t ≥ 1 (vol z ≥ 12 / rate z ≥ 19.5) | ≈ 96 % (events ≈ 97 %) | ≈ 3.2 % | | 0.9 | t ≥ 4 (vol z ≥ 28.5 / rate z ≥ 36) | lower | lower |

A result with confidence = 0.74 at a threshold of 0.75 means the moment is borderline — the combined statistic sits just under one tail unit beyond the calibrated level.

Input data

interface IAggregatedTradeData {
  id:           string;   // Binance aggTradeId
  price:        number;   // Execution price
  qty:          number;   // Trade size (base asset)
  timestamp:    number;   // Unix milliseconds
  isBuyerMaker: boolean;  // true  → sell aggressor (taker sold into bid)
                          // false → buy  aggressor (taker bought ask)
}

isBuyerMaker semantics — this field trips people up. In a limit order book, the maker posts the resting order. When isBuyerMaker = true, the buyer is the maker (passive bid), meaning the seller was the aggressive taker. From an order flow perspective: isBuyerMaker = true → sell aggression, isBuyerMaker = false → buy aggression.

Math

Volume Imbalance

The first quantity derived from raw trades — used as the input series for both CUSUM and BOCPD.

buyVol  = Σ qty_i   for all i where isBuyerMaker = false
sellVol = Σ qty_i   for all i where isBuyerMaker = true
imbalance = (buyVol - sellVol) / (buyVol + sellVol)

Result is in [-1, +1]. +1 = all volume is buy-aggressor. -1 = all volume is sell-aggressor. 0 = balanced. Empty input returns 0.

The key design choice: weighted by qty, not trade count. A single 50 BTC block trade counts 50× more than a 1 BTC retail fill. This makes the imbalance measure resistant to spoofing via many tiny orders.

Each call to detect() computes the imbalance for the full recent window (directional signal, returned as result.imbalance), and also a rolling imbalance series with windowSize-trade sliding windows (used as input to CUSUM and BOCPD):

rolling[i] = imbalance(trades[i - windowSize : i])   for i = windowSize, ..., n

The rolling series converts a raw trade stream into a time series of local imbalances, making it suitable for the sequential change detectors below.

Hawkes Process

Model: Univariate Hawkes process with exponential kernel (Hawkes, 1971).

λ(t) = μ + Σ_{tᵢ < t} α · exp(−β · (t − tᵢ))

• μ > 0 — background intensity (trades/sec in quiet market)
• α > 0 — excitation magnitude (how much each trade boosts future intensity)
• β > 0 — decay rate (how fast the excitement fades)
• α < β — stationarity constraint (subcritical process)

The model says: trades arrive at a baseline rate μ, but each arriving trade triggers a burst of additional arrivals that decays exponentially. This captures the empirical clustering of order flow — a large trade tends to be followed by a flurry of reactive orders.

Unconditional mean:

E[λ] = μ / (1 − α/β)

The ratio α/β is the branching ratio — the expected number of secondary events triggered by one primary event. At α/β = 0.6, each trade triggers on average 0.6 follow-on trades. At α/β → 1 the process becomes supercritical (explosive).

Log-likelihood — O(n) recursive computation:

The naive LL is O(n²) since each event sees all previous events. Ogata (1988) reduced this to O(n) with the recursive compensator trick:

A(0) = 0
A(i) = exp(−β · (tᵢ − tᵢ₋₁)) · (1 + A(i−1))

λ(tᵢ) = μ + α · A(i)

ln L = −μ·T − (α/β)·Σᵢ(1 − exp(−β·(T−tᵢ))) + Σᵢ ln λ(tᵢ)

The second term is the compensator — the expected number of events the model predicts over [0,T]. The third term is the data likelihood under the fitted intensity. Maximum likelihood balances these two.

Fitting via Nelder-Mead MLE:

Parameters [μ, α, β] are estimated by minimising the negative log-likelihood. Starting point:

T    = t_n − t_0          (observation window length)
μ₀   = 0.5 · n/T          (half the empirical rate)
α₀   = 0.4 · n/T          (40% excitation share)
β₀   = n/T                (rate = decay)

Constraints enforced inside the objective: if μ ≤ 0 or α ≤ 0 or β ≤ 0 or α ≥ β, return 1e10 (hard wall). This keeps the optimizer in the subcritical stationary region.

Peak intensity over the detection window:

Instead of evaluating λ at the last event only, the detector takes the maximum λ(tᵢ) seen at any event in the window using the same O(n) recursive trick:

A(0) = 0
A(i) = exp(−β · (tᵢ − tᵢ₋₁)) · (1 + A(i−1))
λ(tᵢ) = μ + α · A(i)

peakLambda = max over i of λ(tᵢ)

This ensures that a burst occurring in the middle of the window is detected even after the kernel has decayed by the last event.

How the detector actually scores (self-calibrated z-channels):

The theoretical baselines above (E[λ], fitted μ) are exported in the math API (hawkesAnomalyScore) but are not what VolumeAnomalyDetector uses — on real trade streams rates fluctuate several-fold between adjacent windows of a perfectly normal market, so "2× the fitted μ" is routine noise (measured: 28 % false alarms on real BTCUSDT). Instead the detector scores robust z of rolling rate and volume-rate statistics against baselines measured on the training window itself:

for horizon h ∈ { rateHorizonSec (5 s), slowHorizonSec (30 s) }:
  rate_i    = trades in (tᵢ − h, tᵢ] / h        (one sample per trade)
  volRate_i = qty    in (tᵢ − h, tᵢ] / h

  z = (max over detection window − median over training) /
      (1.4826 · max(MAD_training, 0.1 · median))

score_volume = max over 4 channels of  σ((z − 12) · 0.4)

The fitted Hawkes model still provides hawkesLambda / lambdaRatio diagnostics (peak conditional intensity vs the training peak), reported in result.stats and signal meta for transparency.

What the volume channel captures: arrival-rate acceleration (many trades per second) and volume waves (few huge trades) at both brief and sustained horizons. It is blind to the direction of trades — direction comes from imbalance.

CUSUM — Sequential Change Detection

Model: Cumulative Sum Control Chart (Page, 1954). Applied to the rolling imbalance series.

The input series is xₜ = |imbalance(window_t)| — absolute imbalance magnitude. The two-sided CUSUM tracks:

S⁺ₜ = max(0,  S⁺_{t-1} + xₜ − μ₀ − k)
S⁻ₜ = max(0,  S⁻_{t-1} − xₜ + μ₀ − k)

• μ₀ — in-control mean of |imbalance| (from training data)
• k — allowable slack = cusumKSigmas · σ₀ (filters out noise below k)
• h — alarm threshold = cusumHSigmas · σ₀ (fires when S ≥ h)

S⁺ accumulates evidence that the series has shifted above its historical mean. S⁻ accumulates evidence of a downward shift. Both are reset to zero after an alarm.

Why absolute imbalance: using |imbalance| instead of signed imbalance means both extreme buy pressure and extreme sell pressure register as anomalies. The direction comes from result.imbalance (signed), not from CUSUM.

Training — parameter estimation:

μ₀   = mean(|imbalance|)       over the training window
σ₀²  = var(|imbalance|)        sample variance
k    = cusumKSigmas · σ₀       (default 0.5σ)
h    = cusumHSigmas · σ₀       (default 5σ)

Average run length under H₀ (ARL₀): the expected number of observations before a false alarm. For Gaussian series, the approximate relationship between h, k and ARL₀ is:

ARL₀ ≈ exp(2·k·h / σ₀²)

At the defaults k = 0.5σ, h = 5σ: ARL₀ ≈ exp(5) ≈ 148. Raising cusumHSigmas to 6 gives ARL₀ ≈ e⁶ ≈ 403. Lowering to 4 gives ARL₀ ≈ exp(4) ≈ 55 — fires quickly but with more false positives.

CUSUM anomaly score:

score_cusum = min(max(S⁺, S⁻) / h,  1)

Linear: 0 at no accumulation, 1 at the alarm boundary. The score reaches 1 exactly when CUSUM would fire, then resets. Between resets it linearly grows with accumulated evidence.

Important nuance — auto-reset on alarm: when the alarm fires, both S⁺ and S⁻ are reset to zero and observation counter n resets. The score thus drops to 0 right after a confirmed alarm. This means score_cusum = 1 is momentary: the next observation after a fire starts fresh. If you see cusumStat close to h but not quite there, the moment is building.

BOCPD — Bayesian Online Changepoint Detection

Model: Adams & MacKay (2007). Computes the posterior distribution over run lengths — the number of observations since the last changepoint — updated online with each new observation.

The fundamental recursion (from the paper):

P(rₜ | x₁:ₜ) ∝ Σ_{rₜ₋₁} P(xₜ | rₜ₋₁, x_{t-r:t}) · P(rₜ | rₜ₋₁) · P(rₜ₋₁ | x₁:ₜ₋₁)

There are two possible transitions at each step:

• Growth P(rₜ = r+1 | rₜ₋₁ = r) = 1 − H — run continues, length grows
• Changepoint P(rₜ = 0 | rₜ₋₁ = r) = H — run resets to zero

The hazard function H = 1 / hazardLambda is constant (geometric / memoryless gaps between changepoints). hazardLambda = 200 means the model expects a changepoint every 200 windows on average.

Underlying observation model — Normal-Gamma conjugate:

Each run-length hypothesis r maintains a separate Normal-Gamma posterior over the mean and precision of {xₜ₋ᵣ, ..., xₜ}. The predictive probability of a new observation given run length r is a Student-t:

p(xₜ | rₜ₋₁ = r, x_{t-r:t}) = Student-t(2αN, μN, βN(κN+1)/(αN·κN))

Where the posterior hyperparameters after n = r observations are updated from prior (μ₀, κ₀, α₀, β₀):

κN = κ₀ + n
αN = α₀ + n/2
μN = (κ₀·μ₀ + n·x̄) / κN
βN = β₀ + 0.5·M₂ + κ₀·n·(x̄ − μ₀)² / (2·κN)

M₂ = Σ(xᵢ − x̄)² is maintained via Welford's online algorithm (numerically stable, O(1) per update).

Prior hyperparameters — derived from training:

μ₀    = mean(|imbalance|)       training window mean
κ₀    = 1                       weak prior (1 pseudo-observation)
α₀    = 1                       weak prior on precision
β₀    = var(|imbalance|)        training window variance

A prior with κ₀ = 1 means the prior contributes the equivalent of one observation. After 10 real observations the likelihood dominates; the prior is only important for brand-new run-length hypotheses (segments just started).

Log-domain computation for numerical stability:

All probabilities are maintained as log-probabilities. The changepoint mass is accumulated via log-sum-exp:

log P(rₜ = 0) = logSumExp over r of [log P(rₜ₋₁ = r) + log p(xₜ | r) + log H]
log P(rₜ = r+1) = log P(rₜ₋₁ = r) + log p(xₜ | r) + log(1-H)

After each update, all log-probs are normalised by subtracting logSumExp(all). This keeps the distribution proper and prevents underflow.

Pruning: hypotheses with log P(rₜ = r) < −30 (probability < 1e-13) are discarded. This bounds memory and computation: in practice the active set stays O(1) to O(few hundred) even after thousands of observations.

Diagnostics returned:

• cpProbability = P(rₜ = 0 | x₁:ₜ) — probability that a changepoint occurred exactly at observation t
• mapRunLength — the run length with highest posterior probability (MAP estimator)

BOCPD anomaly score — relative run-length drop:

cpProbability is approximately equal to the constant prior hazard H = 1/hazardLambda and does not spike at genuine changepoints — it is dominated by the prior, not the data. The real signal is mapRunLength: in a stable process it grows monotonically; a changepoint resets it to near zero.

The score measures the relative drop from the previous step:

drop = clamp((prevRunLength − mapRunLength) / prevRunLength,  0, 1)

score_bocpd = 1 / (1 + exp(−(drop − 0.5) · 8))

Typical values:

• drop = 0 (run length grew — stable) → score ≈ 0.018
• drop = 0.5 (run length halved) → score = 0.50
• drop ≥ 0.9 (e.g. 90 → 1 after reset) → score ≈ 0.98

The sigmoid is centred at drop = 0.5 with steepness 8. The score is taken as the peak over the entire detection window, so changepoints that occurred mid-window are still captured.

What BOCPD captures: regime shifts — moments where the distribution of imbalance itself changes, not just its current level. A market transitioning from choppy balanced flow to sustained directional flow will register here, often before the imbalance crosses an absolute threshold.

Composite Score and Signal Thresholds

The volume channel is the score that matters by default:

score_volume = max over channels of  σ((z − c_ch − u·spanShift) · ln3 / u)
  channels:  z of peak "trades per fast horizon"   (fast arrival bursts)
             z of peak "qty per fast horizon"      (fast volume bursts, block trades)
             z of peak "trades per slow horizon"   (sustained arrival waves)
             z of peak "qty per slow horizon"      (sustained volume waves)

  z    = (peak_detect − median_train) / (1.4826 · max(MAD_train, 0.1·median_train))
  c_ch = max(P85 of the peak-z the baseline itself produced in slow-horizon
             stretches, universal per-channel floor:
             14 for rate channels, 6.5 for volume channels)  → score 0.5
  u    = 5.5 (universal tail unit); c_ch + u                  → score 0.75
  spanShift = log10(max(1, detectionSpan / slowHorizon))      (longer window ⇒
              more independent looks at the max ⇒ higher bar)

confidence_score = w_H · score_volume + w_C · score_cusum + w_B · score_bocpd
                 = score_volume                          (default weights [1, 0, 0])
anomaly          = confidence_score >= confidence_threshold

Everything instrument-specific is chosen on the fly at train() time:

• Horizons — scaled up from the median inter-trade gap for sparse streams (config values are floors; pass explicit values to pin them).
• Yardstick (median/MAD) — measured per channel on the training window.
• Score level c_ch — the baseline's own null distribution of window maxima, floored by the universal mapping validated on a full day of real BTCUSDT trades. A baseline that itself contains recurring bursts absorbs them into its null quantiles, so a repeat of a known pattern scores ≈ 0.5 rather than 1.0.

There are no theoretical thresholds ("2× the fitted μ", "5σ under i.i.d.") — those were tested on real data and misfire badly (28 % false alarms), because real trade streams are heavily autocorrelated and heavy-tailed. The remaining fixed anchors (per-channel floor levels and the tail unit) live in standardized z-space, were found by a dense brute-force search over the full mapping space on the benchmark (44 100 configurations; the chosen point sits on a stable Pareto plateau), and act only as a sensitivity limit — the data can make the detector stricter, never more trigger-happy.

Signals are individual detector firings appended to result.signals when:

| Signal kind | Fires when | Score attached | |-------------|-----------|----------------| | volume_spike | score_volume > 0.5 | max of the four z-channel sigmoids | | imbalance_shift | \|imbalance\| > 0.4 | Raw absolute imbalance | | cusum_alarm | score_cusum > 0.7 | Linear ratio max(S⁺, S⁻) / h (h self-calibrated) | | bocpd_changepoint | score_bocpd > 0.3 | Run-length-drop score above the training noise floor |

A signal in result.signals does not require result.anomaly = true. Signals (and the raw result.scores / result.stats) let you understand why the composite score is what it is.

Design property: at default weights the anomaly decision is purely "is the trade rate or volume rate many robust σ above what this instrument did in the last 15–30 minutes". CUSUM/BOCPD track order-flow imbalance regime shifts — a different phenomenon that correlates poorly with volume anomalies on real data (measured: any additive weight on them reduces recall and adds false alarms). Give them weight only if flow shifts are what you actually want to detect.

Math internals — exported for testing

All internal functions are accessible via the math export path for unit testing and experimentation:

import {
  // Hawkes
  volumeImbalance,
  hawkesLogLikelihood,
  hawkesFit,
  hawkesLambda,
  hawkesPeakLambda,  // max λ(tᵢ) over window — used by the detector
  hawkesAnomalyScore,

  // CUSUM
  cusumFit,
  cusumUpdate,       // returns { state, alarm, preResetState }
  cusumInitState,
  cusumAnomalyScore,
  cusumBatch,

  // BOCPD
  bocpdUpdate,
  bocpdInitState,
  bocpdAnomalyScore,
  bocpdBatch,
  defaultPrior,

  // Optimizer
  nelderMead,
} from 'volume-anomaly/math';

hawkesLogLikelihood(timestamps, params)

Raw Ogata log-likelihood. timestamps must be sorted ascending in seconds. Returns -Infinity if params are invalid (μ ≤ 0, α ≤ 0, β ≤ 0). Does not enforce α < β — that constraint is applied only inside hawkesFit.

hawkesFit(timestamps)

Returns { params, logLik, stationarity, converged }. stationarity = α/β. If timestamps.length < 10, returns a flat Poisson fallback with converged: false.

hawkesLambda(t, timestamps, params)

Evaluates λ(t) at a specific time given a history of prior events. All timestamps must be < t.

hawkesPeakLambda(timestamps, params)

Returns the maximum λ(tᵢ) over all events in timestamps using the O(n) recursive A(i) trick. This is what the detector uses internally instead of hawkesLambda — a burst that decayed by the last event is still captured. hawkesLambda evaluates at a single point; hawkesPeakLambda scans the full window.

cusumUpdate(state, x, params)

Pure function. Returns { state: CusumState, alarm: boolean, preResetState: CusumState }. Does not mutate the input state. preResetState holds the accumulator values before the alarm reset — use it for scoring, since state.sPos/sNeg are zeroed when alarm = true.

bocpdUpdate(state, x, prior, hazardLambda?)

Returns { state, mapRunLength, cpProbability }. The returned state contains pruned log-probability arrays. Pass hazardLambda in windows (same units as your observation index).

nelderMead(f, x0, options?)

Nelder-Mead simplex optimizer. Used internally by Hawkes fitting. Returns { x, fx, iters, converged }.

Optimization

Hawkes parameter estimation uses single-start Nelder-Mead (3 parameters: μ, α, β). The starting point is empirically derived from the data (empirical rate as μ₀) so a single start is typically sufficient. CUSUM and BOCPD use closed-form estimation (sample mean/variance for CUSUM, Welford update for BOCPD) — no optimizer needed.

| Component | Parameters | Optimization | Complexity | |-----------|-----------|--------------|------------| | Hawkes | 3 (μ, α, β) | Nelder-Mead, 1 start, 1000 iter | O(n) per LL eval | | CUSUM | 4 (μ₀, σ₀, k, h) | Closed-form (sample stats) | O(n) | | BOCPD | 4 (μ₀, κ₀, α₀, β₀) | Closed-form (sample stats) | O(n) per update, pruned |

BOCPD update is technically O(r_max) where r_max is the number of surviving run-length hypotheses. The pruning threshold log P < −30 keeps this bounded in practice (typically < 300 hypotheses even after 10,000 observations).

Training and detection window sizes

train() — historical window

The time span matters more than the trade count. The volume-channel baselines (median/MAD of rolling rates) are meaningful only if the training window covers enough market time to define "typical":

| historical spans | Volume-channel baselines | Notes | |--------------------|--------------------------|-------| | < 1 minute | Unreliable | The slow horizon (30 s) has almost no independent samples; a hot pre-burst minute masks the burst itself | | 1–5 minutes | Weak | Usable for the fast horizon only | | 15–30 minutes | Good | Recommended. This is what the calibration in test/eval.test.ts uses | | Multiple hours | Robust but stale | Beware regime staleness — baselines may span very different market conditions |

Never build the baseline as "the last N trades": a count-based window adapts its duration to market pace — right before a burst it covers seconds of already-hot market and the inflated baseline masks the very anomaly being detected, while in a calm period it covers narrow minutes and routine wiggles look anomalous. Feed a fixed time span of recent history instead (internally the expensive fits are capped: Hawkes MLE uses the most recent 2 000 trades, the imbalance series the most recent 1 000; the O(n) rate baselines use everything).

train() throws below 50 trades. The training window should represent normal market conditions: an anomaly inside the baseline inflates the yardstick and reduces sensitivity to similar events.

detect() — recent window

The volume channel needs the recent window to span at least rateHorizonSec (5 s default) — shorter windows fall back to a single whole-window rate sample. CUSUM/BOCPD additionally need trades ≥ windowSize for their rolling imbalance series; when the window is shorter, score weights are renormalized over the detectors that actually ran, so the volume channel can still reach confidence 1 alone.

| Trades in recent | Notes | |--------------------|-------| | < 5 s of market time | Fast channel falls back to whole-window rate; slow channel off | | 100–300 trades | Recommended — what the real-data calibration uses | | 300–1000 trades | Fine; peak-based statistics don't dilute with window growth |

Rule of thumb: call detect() every few seconds with the last ~300 trades; re-train() every few minutes with the trailing 15–30 minutes.

windowSize guidance

| windowSize | Sensitivity | Lag | Minimum train() | Minimum detect() for full signal | |-------------|-------------|-----|-------------------|-------------------------------------| | 20 | Very high | Low | 50 trades (code minimum) | 40 trades | | 50 (default) | Balanced | Moderate | 100 trades (recommended) | 100 trades | | 100 | Lower | Higher | 200 trades | 200 trades | | 200 | Low | High | 400 trades | 400 trades |

On high-volume pairs (BTC/USDT perpetual), 50 trades may span only 1–2 seconds. On low-volume pairs, 50 trades may span minutes. Calibrate to the effective time scale that matters for your entry.

Integration with garch

Typical workflow combining both libraries:

import { predict }                  from 'garch';
import { detect, VolumeAnomalyDetector } from 'volume-anomaly';

// 1. Train volume detector once per session
const detector = new VolumeAnomalyDetector({ windowSize: 50 });
detector.train(await getAggregatedTrades('BTCUSDT', 2000));

// 2. On each new candle close:
async function onCandle(candles: Candle[], recentTrades: IAggregatedTradeData[]) {
  // Volume anomaly — entry timing
  const vol = detector.detect(recentTrades, 0.75);

  if (!vol.anomaly) return; // not an anomalous moment

  // Directional filter from your fundamental analysis
  const isBuySignal = vol.imbalance > 0.3 && myFundamentalBullish();
  if (!isBuySignal) return;

  // garch — TP/SL sizing
  const { upperPrice, lowerPrice, sigma, reliable } = predict(candles, '15m');
  if (!reliable) return;

  const entry = currentPrice;
  const tp    = upperPrice;         // +1σ target (default confidence 0.6827)
  const sl    = lowerPrice;         // -1σ stop

  // Or use 95% VaR for wider stop:
  const { lowerPrice: sl95 } = predict(candles, '15m', undefined, 0.95);

  placeOrder({ entry, tp, sl: sl95 });
}

garch.predict answers: how big is the next normal move? volume-anomaly.detect answers: is this moment abnormal enough to act on? They are complementary and independent.

Tests

735 tests across 18 test files. All passing. 100% statement/function/line coverage, 98.72% branch (two unreachable ?? guards).

| File | Tests | What is covered | |------|-------|-----------------| | hawkes.test.ts | 20 | Imbalance formula, LL computation, MLE fitting, λ evaluation and decay, anomaly score monotonicity and supercritical clamp | | cusum.test.ts | 15 | Parameter estimation, state update (pure function), accumulation, alarm + reset, score range, batch detection | | bocpd.test.ts | 13 | Init state, update, probability normalisation, run length growth in stable regime, CP spike on distribution shift, immutability, batch | | detector.test.ts | 20 | Pre-train guard, isTrained flag, minimum training size, DetectionResult fields, confidence range, empty window, signal score range | | detect.test.ts | 36 | End-to-end anomaly detection, confidence thresholds, signal composition, edge inputs | | seeded.test.ts | 67 | Deterministic seeded scenarios covering long/short/neutral bursts across parameter space | | predict.test.ts | 24 | Direction assignment, trained imbalanceThreshold, imbalancePercentile config, trending vs balanced threshold | | invariants.test.ts | 29 | Monotonicity, score bounds, immutability, score weight validation | | adversarial.test.ts | 58 | Adversarial inputs: NaN propagation, extreme values, Inf timestamps, zero-qty trades | | falsepositive.test.ts | 18 | Scenarios that must NOT trigger: gradual drift, HFT clusters, trending market, whale trades, overnight gaps | | edgecases.test.ts | 80 | Boundary conditions, signal threshold exact values (strict >), detect < windowSize bypass, train twice, cusumBatch multiple alarms | | realdata.test.ts | 23 | Real BTCUSDT-2025-03-01 data: 4 spike windows + 1 calm baseline | | robustness.test.ts | 66 | Mathematical invariants: range/symmetry/monotonicity for all functions, BOCPD normalisation Σexp(lp) ≤ 1, 100-case property-based detector test | | extreme.test.ts | 52 | Stuck-at-extremum: hazardLambda edge cases, μ = 0, extreme β, degenerate Nelder-Mead, Welford drift, β₀ = 0 | | newextreme.test.ts | 58 | NaN propagation in CUSUM/BOCPD, hawkesAnomalyScore extremes, cusumAnomalyScore h = NaN, prevRL = Inf/NaN, kappa0 = 0 | | thirdextreme.test.ts | 74 | hazardLambda = 0 collapse, β ≤ 0, hawkesFit n = 0/T = 0, hawkesPeakLambda n = 1/β = 0, volumeImbalance NaN qty, cusumFit NaN filter | | fourthextreme.test.ts | 63 | hawkesAnomalyScore NaN peak + valid params, cusumUpdate NaN params, cusumAnomalyScore NaN state, bocpdUpdate beta0 = Inf, Infinity qty | | perf.test.ts | 19 | Latency P95 bounds, throughput (detect(200) ≥ 800/s), scaling ratios, stability over 500 sequential calls |

npm test

License

MIT