Skip to content

API Reference

This page combines a short API map with generated reference documentation from the package docstrings.

Which API Should I Use?

Workflow API
Standard one-shot detection cfad.detect
Custom rolling detector configuration cfad.detection.RollingDetector
Online updates cfad.detection.StreamDetector
Train/test temporal evaluation cfad.backtest.WalkForwardBacktest
ECF model comparison cfad.compare_models
ECF goodness-of-fit tests and distances cfad.gof
Parameter sweeps cfad.sensitivity
Plotting diagnostics cfad.viz

Return Objects

detect() and RollingDetector.fit_transform() return an AnomalyReport. The most important fields are scores, cusum_pos, cusum_neg, alarm_indices, window_end_indices, mu0, sigma0, and threshold.

WalkForwardBacktest.run() returns a BacktestResult, which can be summarized with summary() or converted to a fold-concatenated DataFrame with to_dataframe().

High-Level API

Public API entry points for common CFAD workflows.

AnomalyReport dataclass 

AnomalyReport(scores: NDArray[float64], cusum_pos: NDArray[float64], cusum_neg: NDArray[float64], alarm_indices: NDArray[int64], window_end_indices: NDArray[int64], dates: Optional[DatetimeIndex] = None, mu0: float = 0.0, sigma0: float = 1.0, threshold: float = 5.0)

Container for detector output.

alarm_dates property 

alarm_dates: Optional[DatetimeIndex]

Return dates corresponding to alarm-window endpoints.

window_end_indices use Python's half-open convention and therefore point one position beyond the final observation in each rolling window.

summary 

summary() -> str

Return a compact human-readable report summary.

Source code in cfad/detection.py 
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def summary(self) -> str:
    """Return a compact human-readable report summary."""
    n_alarms = len(self.alarm_indices)
    lines = [
        "CFAD Anomaly Report",
        f"  Windows evaluated : {len(self.scores)}",
        f"  Alarms fired      : {n_alarms}",
        f"  In-control mean   : {self.mu0:.4f}",
        f"  In-control std    : {self.sigma0:.4f}",
        f"  CUSUM threshold   : {self.threshold}",
    ]
    alarm_dates = self.alarm_dates
    if alarm_dates is not None and len(alarm_dates) > 0:
        lines.append(f"  First alarm       : {alarm_dates[0].date()}")
    return "\n".join(lines)

RollingDetector 

RollingDetector(window: int = 60, xi_min: float = -10.0, xi_max: float = 10.0, n_xi: int = 128, step: int = 1, calibration_frac: float = 0.3, k: float = 0.5, h: float = 5.0)

Rolling ECF shape detector followed by a two-sided Page-CUSUM.

Source code in cfad/detection.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def __init__(
    self,
    window: int = 60,
    xi_min: float = -10.0,
    xi_max: float = 10.0,
    n_xi: int = 128,
    step: int = 1,
    calibration_frac: float = 0.3,
    k: float = 0.5,
    h: float = 5.0,
) -> None:
    if window <= 1:
        raise ValueError("window must be greater than 1")
    if xi_max <= xi_min:
        raise ValueError("xi_max must be greater than xi_min")
    if n_xi < 4:
        raise ValueError("n_xi must be at least 4")
    if step <= 0:
        raise ValueError("step must be positive")
    if not (0.0 < calibration_frac < 1.0):
        raise ValueError("calibration_frac must lie in (0, 1)")
    if k < 0.0:
        raise ValueError("k must be non-negative")
    if h <= 0.0:
        raise ValueError("h must be positive")

    self.window = int(window)
    self.xi_grid = np.linspace(xi_min, xi_max, int(n_xi), dtype=np.float64)
    self.step = int(step)
    self.calibration_frac = float(calibration_frac)
    self.k = float(k)
    self.h = float(h)
    self._mu0: Optional[float] = None
    self._sigma0: Optional[float] = None

score_windows 

score_windows(returns: NDArray[float64]) -> tuple[NDArray[np.float64], NDArray[np.int64]]

Compute rolling ECF-shape scores without fitting CUSUM calibration.

This separation is useful for walk-forward evaluation: a test fold can be scored without estimating any in-control parameter from the test data.

Source code in cfad/detection.py 
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
def score_windows(
    self,
    returns: NDArray[np.float64],
) -> tuple[NDArray[np.float64], NDArray[np.int64]]:
    """Compute rolling ECF-shape scores without fitting CUSUM calibration.

    This separation is useful for walk-forward evaluation: a test fold can
    be scored without estimating any in-control parameter from the test data.
    """
    values = np.asarray(returns, dtype=np.float64)
    if values.ndim != 1:
        raise ValueError("returns must be one-dimensional")
    if values.size < self.window:
        raise ValueError("returns must contain at least one full window")
    if not np.all(np.isfinite(values)):
        raise ValueError("returns must contain only finite values")

    ecf_mat, end_idx = rolling_ecf(
        values,
        self.xi_grid,
        self.window,
        self.step,
    )
    means, stds = _rolling_moments(values, self.window, self.step)
    scores = gaussian_ecf_distance_scores(
        ecf_mat,
        self.xi_grid,
        means,
        stds,
    )
    return (
        np.asarray(scores, dtype=np.float64),
        np.asarray(end_idx, dtype=np.int64),
    )

apply_calibration 

apply_calibration(scores: NDArray[float64], mu0: float, sigma0: float) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.int64]]

Apply this detector's CUSUM settings to externally calibrated scores.

Source code in cfad/detection.py 
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
def apply_calibration(
    self,
    scores: NDArray[np.float64],
    mu0: float,
    sigma0: float,
) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.int64]]:
    """Apply this detector's CUSUM settings to externally calibrated scores."""
    cusum_fn = _cusum_c if _HAS_C_EXT else _cusum_python
    positive, negative, alarms = cusum_fn(
        np.asarray(scores, dtype=np.float64),
        float(mu0),
        float(sigma0),
        self.k,
        self.h,
    )
    return (
        np.asarray(positive, dtype=np.float64),
        np.asarray(negative, dtype=np.float64),
        np.asarray(alarms, dtype=np.int64),
    )

fit_transform 

fit_transform(returns: NDArray[float64], dates: Optional[DatetimeIndex] = None) -> AnomalyReport

Score windows, estimate in-control moments, and apply CUSUM.

Source code in cfad/detection.py 
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
def fit_transform(
    self,
    returns: NDArray[np.float64],
    dates: Optional[pd.DatetimeIndex] = None,
) -> AnomalyReport:
    """Score windows, estimate in-control moments, and apply CUSUM."""
    scores, end_idx = self.score_windows(returns)

    n_cal = max(10, int(self.calibration_frac * len(scores)))
    n_cal = min(n_cal, len(scores))
    self._mu0 = float(np.mean(scores[:n_cal]))
    sigma0 = float(np.std(scores[:n_cal], ddof=1)) if n_cal > 1 else 0.0
    self._sigma0 = max(sigma0, 1e-12)

    positive, negative, alarms = self.apply_calibration(
        scores,
        self._mu0,
        self._sigma0,
    )
    return AnomalyReport(
        scores=scores,
        cusum_pos=positive,
        cusum_neg=negative,
        alarm_indices=alarms,
        window_end_indices=end_idx,
        dates=dates,
        mu0=self._mu0,
        sigma0=self._sigma0,
        threshold=self.h,
    )

detect 

detect(returns: Union[NDArray, Series], window: int = 60, xi_range: tuple[float, float] = (-10.0, 10.0), n_xi: int = 128, step: int = 1, calibration_frac: float = 0.3, k: float = 0.5, h: float = 5.0, *, height: float | None = None) -> AnomalyReport

Detect distributional-shape changes in a financial return series.

Each rolling empirical characteristic function is compared with the Gaussian characteristic function fitted to the same window. The resulting real-frequency L2 distance is monitored with a two-sided Page-CUSUM.

Parameters:

Name Type Description Default
returns array - like or Series

One-dimensional return series.

 required
window int

Rolling window size for ECF estimation.

 60
xi_range tuple[float, float]

Real-frequency grid bounds.

 (-10, 10)
n_xi int

Number of frequency grid points.

 128
step int

Rolling step.

 1
calibration_frac float

Fraction of score windows used to estimate the in-control score mean and standard deviation.

 0.3
k float

Dimensionless Page-CUSUM reference value on standardized scores.

 0.5
h float

CUSUM decision threshold.

 5.0
height (float or None, keyword - only)

Deprecated compatibility argument from the former empirical-contour implementation. It is ignored by the corrected real-frequency score and will be removed in a future breaking release.

 None

Returns:

Type Description
AnomalyReport

Scores, CUSUM statistics, alarm indices, and optional dates.
Source code in cfad/api.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
def detect(
    returns: Union[NDArray, pd.Series],
    window: int = 60,
    xi_range: tuple[float, float] = (-10.0, 10.0),
    n_xi: int = 128,
    step: int = 1,
    calibration_frac: float = 0.3,
    k: float = 0.5,
    h: float = 5.0,
    *,
    height: float | None = None,
) -> AnomalyReport:
    """Detect distributional-shape changes in a financial return series.

    Each rolling empirical characteristic function is compared with the
    Gaussian characteristic function fitted to the same window. The resulting
    real-frequency L2 distance is monitored with a two-sided Page-CUSUM.

    Parameters
    ----------
    returns : array-like or pandas.Series
        One-dimensional return series.
    window : int, default=60
        Rolling window size for ECF estimation.
    xi_range : tuple[float, float], default=(-10, 10)
        Real-frequency grid bounds.
    n_xi : int, default=128
        Number of frequency grid points.
    step : int, default=1
        Rolling step.
    calibration_frac : float, default=0.3
        Fraction of score windows used to estimate the in-control score mean and
        standard deviation.
    k : float, default=0.5
        Dimensionless Page-CUSUM reference value on standardized scores.
    h : float, default=5.0
        CUSUM decision threshold.
    height : float or None, keyword-only
        Deprecated compatibility argument from the former empirical-contour
        implementation. It is ignored by the corrected real-frequency score and
        will be removed in a future breaking release.

    Returns
    -------
    AnomalyReport
        Scores, CUSUM statistics, alarm indices, and optional dates.
    """
    if height is not None:
        warnings.warn(
            "height is deprecated and ignored; CFAD now uses a real-frequency "
            "ECF shape score. Tune xi_range instead.",
            DeprecationWarning,
            stacklevel=2,
        )

    dates = None
    if isinstance(returns, pd.Series):
        dates = returns.index if hasattr(returns.index, "to_pydatetime") else None
        returns_arr = returns.to_numpy(dtype=np.float64)
    else:
        returns_arr = np.asarray(returns, dtype=np.float64)

    detector = RollingDetector(
        window=window,
        xi_min=xi_range[0],
        xi_max=xi_range[1],
        n_xi=n_xi,
        step=step,
        calibration_frac=calibration_frac,
        k=k,
        h=h,
    )
    return detector.fit_transform(returns_arr, dates=dates)

compare_models 

compare_models(returns: NDArray[float64], xi: Optional[NDArray[float64]] = None) -> dict[str, object]

Fit Gaussian and NIG models and compare real-frequency ECF distance.

This model comparison is descriptive evidence about distributional fit. It must not be interpreted as a test for branch cuts or population-CF singularities from a finite-sample empirical characteristic function.

Source code in cfad/api.py 
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
def compare_models(
    returns: NDArray[np.float64],
    xi: Optional[NDArray[np.float64]] = None,
) -> dict[str, object]:
    """Fit Gaussian and NIG models and compare real-frequency ECF distance.

    This model comparison is descriptive evidence about distributional fit. It
    must not be interpreted as a test for branch cuts or population-CF
    singularities from a finite-sample empirical characteristic function.
    """
    from cfad.empirical_cf import ecf_at
    from cfad.models.gaussian import GaussianCF
    from cfad.models.nig import NIGCF

    values = np.asarray(returns, dtype=np.float64)
    if values.ndim != 1 or values.size < 2:
        raise ValueError("returns must be one-dimensional with at least 2 values")
    if not np.all(np.isfinite(values)):
        raise ValueError("returns must contain only finite values")

    grid = (
        np.linspace(-15.0, 15.0, 256, dtype=np.float64)
        if xi is None
        else np.asarray(xi, dtype=np.float64)
    )
    if grid.ndim != 1 or grid.size < 4:
        raise ValueError("xi must be one-dimensional with at least 4 values")

    empirical = ecf_at(values, grid)
    gaussian = GaussianCF().fit(values)
    nig = NIGCF().fit(values)

    gaussian_distance = float(np.mean(np.abs(empirical - gaussian.cf(grid)) ** 2))
    nig_distance = float(np.mean(np.abs(empirical - nig.cf(grid)) ** 2))

    return {
        "gaussian": {
            "model": gaussian,
            "ecf_l2": gaussian_distance,
            "aic": gaussian.aic(values),
        },
        "nig": {
            "model": nig,
            "ecf_l2": nig_distance,
            "aic": nig.aic(values),
        },
        "winner": "nig" if nig_distance < gaussian_distance else "gaussian",
    }

Detection Objects

Rolling ECF shape detector followed by a two-sided Page-CUSUM.

Source code in cfad/detection.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def __init__(
    self,
    window: int = 60,
    xi_min: float = -10.0,
    xi_max: float = 10.0,
    n_xi: int = 128,
    step: int = 1,
    calibration_frac: float = 0.3,
    k: float = 0.5,
    h: float = 5.0,
) -> None:
    if window <= 1:
        raise ValueError("window must be greater than 1")
    if xi_max <= xi_min:
        raise ValueError("xi_max must be greater than xi_min")
    if n_xi < 4:
        raise ValueError("n_xi must be at least 4")
    if step <= 0:
        raise ValueError("step must be positive")
    if not (0.0 < calibration_frac < 1.0):
        raise ValueError("calibration_frac must lie in (0, 1)")
    if k < 0.0:
        raise ValueError("k must be non-negative")
    if h <= 0.0:
        raise ValueError("h must be positive")

    self.window = int(window)
    self.xi_grid = np.linspace(xi_min, xi_max, int(n_xi), dtype=np.float64)
    self.step = int(step)
    self.calibration_frac = float(calibration_frac)
    self.k = float(k)
    self.h = float(h)
    self._mu0: Optional[float] = None
    self._sigma0: Optional[float] = None

window instance-attribute 

window = int(window)

xi_grid instance-attribute 

xi_grid = np.linspace(xi_min, xi_max, int(n_xi), dtype=np.float64)

step instance-attribute 

step = int(step)

calibration_frac instance-attribute 

calibration_frac = float(calibration_frac)

k instance-attribute 

k = float(k)

h instance-attribute 

h = float(h)

_mu0 instance-attribute 

_mu0: Optional[float] = None

_sigma0 instance-attribute 

_sigma0: Optional[float] = None

score_windows 

score_windows(returns: NDArray[float64]) -> tuple[NDArray[np.float64], NDArray[np.int64]]

Compute rolling ECF-shape scores without fitting CUSUM calibration.

This separation is useful for walk-forward evaluation: a test fold can be scored without estimating any in-control parameter from the test data.

Source code in cfad/detection.py 
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
def score_windows(
    self,
    returns: NDArray[np.float64],
) -> tuple[NDArray[np.float64], NDArray[np.int64]]:
    """Compute rolling ECF-shape scores without fitting CUSUM calibration.

    This separation is useful for walk-forward evaluation: a test fold can
    be scored without estimating any in-control parameter from the test data.
    """
    values = np.asarray(returns, dtype=np.float64)
    if values.ndim != 1:
        raise ValueError("returns must be one-dimensional")
    if values.size < self.window:
        raise ValueError("returns must contain at least one full window")
    if not np.all(np.isfinite(values)):
        raise ValueError("returns must contain only finite values")

    ecf_mat, end_idx = rolling_ecf(
        values,
        self.xi_grid,
        self.window,
        self.step,
    )
    means, stds = _rolling_moments(values, self.window, self.step)
    scores = gaussian_ecf_distance_scores(
        ecf_mat,
        self.xi_grid,
        means,
        stds,
    )
    return (
        np.asarray(scores, dtype=np.float64),
        np.asarray(end_idx, dtype=np.int64),
    )

apply_calibration 

apply_calibration(scores: NDArray[float64], mu0: float, sigma0: float) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.int64]]

Apply this detector's CUSUM settings to externally calibrated scores.

Source code in cfad/detection.py 
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
def apply_calibration(
    self,
    scores: NDArray[np.float64],
    mu0: float,
    sigma0: float,
) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.int64]]:
    """Apply this detector's CUSUM settings to externally calibrated scores."""
    cusum_fn = _cusum_c if _HAS_C_EXT else _cusum_python
    positive, negative, alarms = cusum_fn(
        np.asarray(scores, dtype=np.float64),
        float(mu0),
        float(sigma0),
        self.k,
        self.h,
    )
    return (
        np.asarray(positive, dtype=np.float64),
        np.asarray(negative, dtype=np.float64),
        np.asarray(alarms, dtype=np.int64),
    )

fit_transform 

fit_transform(returns: NDArray[float64], dates: Optional[DatetimeIndex] = None) -> AnomalyReport

Score windows, estimate in-control moments, and apply CUSUM.

Source code in cfad/detection.py 
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
def fit_transform(
    self,
    returns: NDArray[np.float64],
    dates: Optional[pd.DatetimeIndex] = None,
) -> AnomalyReport:
    """Score windows, estimate in-control moments, and apply CUSUM."""
    scores, end_idx = self.score_windows(returns)

    n_cal = max(10, int(self.calibration_frac * len(scores)))
    n_cal = min(n_cal, len(scores))
    self._mu0 = float(np.mean(scores[:n_cal]))
    sigma0 = float(np.std(scores[:n_cal], ddof=1)) if n_cal > 1 else 0.0
    self._sigma0 = max(sigma0, 1e-12)

    positive, negative, alarms = self.apply_calibration(
        scores,
        self._mu0,
        self._sigma0,
    )
    return AnomalyReport(
        scores=scores,
        cusum_pos=positive,
        cusum_neg=negative,
        alarm_indices=alarms,
        window_end_indices=end_idx,
        dates=dates,
        mu0=self._mu0,
        sigma0=self._sigma0,
        threshold=self.h,
    )

Online detector using the same ECF-shape score as RollingDetector.

Source code in cfad/detection.py 
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
def __init__(
    self,
    window: int,
    xi_min: float,
    xi_max: float,
    n_xi: int,
    mu0: Optional[float] = None,
    sigma0: Optional[float] = None,
    warmup: Optional[int] = None,
    k: float = 0.5,
    h: float = 5.0,
) -> None:
    if window <= 1:
        raise ValueError("window must be greater than 1")
    if xi_max <= xi_min:
        raise ValueError("xi_max must be greater than xi_min")
    if n_xi < 4:
        raise ValueError("n_xi must be at least 4")
    if (mu0 is None) != (sigma0 is None):
        raise ValueError("mu0 and sigma0 must be both provided or both omitted")
    if sigma0 is not None and sigma0 <= 0.0:
        raise ValueError("sigma0 must be positive")
    if warmup is not None and warmup < 0:
        raise ValueError("warmup must be non-negative")
    if k < 0.0:
        raise ValueError("k must be non-negative")
    if h <= 0.0:
        raise ValueError("h must be positive")

    self.window = int(window)
    self.xi_grid = np.linspace(xi_min, xi_max, int(n_xi), dtype=np.float64)
    self.k = float(k)
    self.h = float(h)
    self._fixed_mu0 = None if mu0 is None else float(mu0)
    self._fixed_sigma0 = None if sigma0 is None else float(sigma0)
    self.warmup = self.window if warmup is None else int(warmup)

    self._buffer: deque[float] = deque(maxlen=self.window)
    self._warmup_scores: list[float] = []
    self.n_obs = 0
    self.cusum_pos = 0.0
    self.cusum_neg = 0.0
    self.mu0: Optional[float] = None
    self.sigma0: Optional[float] = None
    self._calibrated = False
    self.reset()

window instance-attribute 

window = int(window)

xi_grid instance-attribute 

xi_grid = np.linspace(xi_min, xi_max, int(n_xi), dtype=np.float64)

k instance-attribute 

k = float(k)

h instance-attribute 

h = float(h)

_fixed_mu0 instance-attribute 

_fixed_mu0 = None if mu0 is None else float(mu0)

_fixed_sigma0 instance-attribute 

_fixed_sigma0 = None if sigma0 is None else float(sigma0)

warmup instance-attribute 

warmup = self.window if warmup is None else int(warmup)

_buffer instance-attribute 

_buffer: deque[float] = deque(maxlen=self.window)

_warmup_scores instance-attribute 

_warmup_scores: list[float] = []

n_obs instance-attribute 

n_obs = 0

cusum_pos instance-attribute 

cusum_pos = 0.0

cusum_neg instance-attribute 

cusum_neg = 0.0

mu0 instance-attribute 

mu0: Optional[float] = None

sigma0 instance-attribute 

sigma0: Optional[float] = None

_calibrated instance-attribute 

_calibrated = False

is_calibrated property 

is_calibrated: bool

Return whether score calibration is complete.

_make_output 

_make_output(score: float, alarm: bool) -> dict[str, float | bool | int]

Create the public result object for one streamed observation.

Source code in cfad/detection.py 
327
328
329
330
331
332
333
334
335
336
337
338
339
340
def _make_output(
    self,
    score: float,
    alarm: bool,
) -> dict[str, float | bool | int]:
    """Create the public result object for one streamed observation."""
    return {
        "score": float(score),
        "cusum_pos": float(self.cusum_pos),
        "cusum_neg": float(self.cusum_neg),
        "alarm": bool(alarm),
        "n_obs": int(self.n_obs),
        "calibrated": bool(self._calibrated),
    }

_calibrate_if_ready 

_calibrate_if_ready() -> None

Estimate in-control score moments once enough score windows exist.

Source code in cfad/detection.py 
342
343
344
345
346
347
348
349
350
def _calibrate_if_ready(self) -> None:
    """Estimate in-control score moments once enough score windows exist."""
    if self._calibrated or len(self._warmup_scores) < self.warmup:
        return
    warmup_arr = np.asarray(self._warmup_scores, dtype=np.float64)
    self.mu0 = float(np.mean(warmup_arr))
    sigma0 = float(np.std(warmup_arr, ddof=1)) if warmup_arr.size > 1 else 0.0
    self.sigma0 = max(sigma0, 1e-12)
    self._calibrated = True

update 

update(r: float) -> dict[str, float | bool | int]

Ingest one return observation and update the sequential detector.

Source code in cfad/detection.py 
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
def update(self, r: float) -> dict[str, float | bool | int]:
    """Ingest one return observation and update the sequential detector."""
    if not np.isfinite(r):
        raise ValueError("streamed returns must be finite")

    self.n_obs += 1
    self._buffer.append(float(r))
    if len(self._buffer) < self.window:
        return self._make_output(np.nan, alarm=False)

    sample = np.asarray(self._buffer, dtype=np.float64)
    ecf_vec = ecf_at(sample, self.xi_grid)
    mean = np.asarray([float(np.mean(sample))], dtype=np.float64)
    std = np.asarray([float(np.std(sample, ddof=1))], dtype=np.float64)
    score = float(
        gaussian_ecf_distance_scores(
            ecf_vec[np.newaxis, :],
            self.xi_grid,
            mean,
            std,
        )[0]
    )

    if not self._calibrated:
        self._warmup_scores.append(score)
        self._calibrate_if_ready()
        return self._make_output(np.nan, alarm=False)

    if self.mu0 is None or self.sigma0 is None:
        raise RuntimeError("detector is calibrated but mu0/sigma0 is missing")

    z = (score - self.mu0) / self.sigma0
    self.cusum_pos = max(0.0, self.cusum_pos + z - self.k)
    self.cusum_neg = max(0.0, self.cusum_neg - z - self.k)
    alarm = bool(self.cusum_pos > self.h or self.cusum_neg > self.h)
    if alarm:
        self.cusum_pos = 0.0
        self.cusum_neg = 0.0
    return self._make_output(score, alarm=alarm)

update_batch 

update_batch(returns: NDArray[float64]) -> list[dict[str, float | bool | int]]

Process a one-dimensional array through :meth:update.

Source code in cfad/detection.py 
392
393
394
395
396
397
398
399
400
def update_batch(
    self,
    returns: NDArray[np.float64],
) -> list[dict[str, float | bool | int]]:
    """Process a one-dimensional array through :meth:`update`."""
    values = np.asarray(returns, dtype=np.float64)
    if values.ndim != 1:
        raise ValueError("returns must be one-dimensional")
    return [self.update(float(value)) for value in values]

reset 

reset() -> None

Reset buffer, CUSUM state, and score calibration.

Source code in cfad/detection.py 
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
def reset(self) -> None:
    """Reset buffer, CUSUM state, and score calibration."""
    self._buffer.clear()
    self._warmup_scores = []
    self.n_obs = 0
    self.cusum_pos = 0.0
    self.cusum_neg = 0.0

    if self._fixed_mu0 is not None and self._fixed_sigma0 is not None:
        self.mu0 = self._fixed_mu0
        self.sigma0 = self._fixed_sigma0
        self._calibrated = True
    else:
        self.mu0 = None
        self.sigma0 = None
        self._calibrated = False

Container for detector output.

scores instance-attribute 

scores: NDArray[float64]

cusum_pos instance-attribute 

cusum_pos: NDArray[float64]

cusum_neg instance-attribute 

cusum_neg: NDArray[float64]

alarm_indices instance-attribute 

alarm_indices: NDArray[int64]

window_end_indices instance-attribute 

window_end_indices: NDArray[int64]

dates class-attribute instance-attribute 

dates: Optional[DatetimeIndex] = None

mu0 class-attribute instance-attribute 

mu0: float = 0.0

sigma0 class-attribute instance-attribute 

sigma0: float = 1.0

threshold class-attribute instance-attribute 

threshold: float = 5.0

alarm_dates property 

alarm_dates: Optional[DatetimeIndex]

Return dates corresponding to alarm-window endpoints.

window_end_indices use Python's half-open convention and therefore point one position beyond the final observation in each rolling window.

summary 

summary() -> str

Return a compact human-readable report summary.

Source code in cfad/detection.py 
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def summary(self) -> str:
    """Return a compact human-readable report summary."""
    n_alarms = len(self.alarm_indices)
    lines = [
        "CFAD Anomaly Report",
        f"  Windows evaluated : {len(self.scores)}",
        f"  Alarms fired      : {n_alarms}",
        f"  In-control mean   : {self.mu0:.4f}",
        f"  In-control std    : {self.sigma0:.4f}",
        f"  CUSUM threshold   : {self.threshold}",
    ]
    alarm_dates = self.alarm_dates
    if alarm_dates is not None and len(alarm_dates) > 0:
        lines.append(f"  First alarm       : {alarm_dates[0].date()}")
    return "\n".join(lines)

Backtesting

Walk-forward evaluation with train-only score calibration.

Source code in cfad/backtest.py 
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def __init__(
    self,
    detector_kwargs: dict,
    n_folds: int = 5,
    train_frac: float = 0.6,
    expanding: bool = True,
) -> None:
    if n_folds <= 0:
        raise ValueError("n_folds must be positive")
    if not (0.0 < train_frac < 1.0):
        raise ValueError("train_frac must be in (0, 1)")
    if not isinstance(detector_kwargs, dict):
        raise TypeError("detector_kwargs must be a dict")

    kwargs = dict(detector_kwargs)
    kwargs.pop("height", None)
    self.detector_kwargs = kwargs
    self.n_folds = int(n_folds)
    self.train_frac = float(train_frac)
    self.expanding = bool(expanding)

detector_kwargs instance-attribute 

detector_kwargs = kwargs

n_folds instance-attribute 

n_folds = int(n_folds)

train_frac instance-attribute 

train_frac = float(train_frac)

expanding instance-attribute 

expanding = bool(expanding)

_split_folds 

_split_folds(n_obs: int, window: int) -> list[tuple[int, int, int, int]]

Create non-overlapping test folds as half-open index intervals.

Source code in cfad/backtest.py 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
def _split_folds(
    self,
    n_obs: int,
    window: int,
) -> list[tuple[int, int, int, int]]:
    """Create non-overlapping test folds as half-open index intervals."""
    initial_train = max(int(np.floor(self.train_frac * n_obs)), window)
    if initial_train >= n_obs:
        raise ValueError("Initial training fold leaves no samples for testing")

    remaining = n_obs - initial_train
    if remaining < self.n_folds:
        raise ValueError("Not enough samples in the test region for n_folds")

    base, extra = divmod(remaining, self.n_folds)
    train_size_fixed = initial_train
    folds: list[tuple[int, int, int, int]] = []
    test_start = initial_train

    for i in range(self.n_folds):
        fold_len = base + (1 if i < extra else 0)
        test_end = test_start + fold_len
        train_end = test_start
        train_start = 0 if self.expanding else max(0, train_end - train_size_fixed)

        if train_end - train_start < window:
            raise ValueError("Training fold is shorter than detector window")
        if test_end - test_start < window:
            raise ValueError("Test fold is shorter than detector window")

        folds.append((train_start, train_end, test_start, test_end))
        test_start = test_end

    return folds

_sanitize_sigma staticmethod 

_sanitize_sigma(value: float) -> float

Ensure a strictly positive finite calibration scale.

Source code in cfad/backtest.py 
159
160
161
162
163
164
@staticmethod
def _sanitize_sigma(value: float) -> float:
    """Ensure a strictly positive finite calibration scale."""
    if not np.isfinite(value) or value <= 0.0:
        return 1e-12
    return float(value)

run 

run(returns: NDArray[float64], dates: Optional[DatetimeIndex] = None) -> BacktestResult

Execute a leakage-free walk-forward evaluation.

Each fold estimates mu0 and sigma0 only from training-window scores. Test-window ECF scores are produced by score_windows and the frozen training calibration is then applied via apply_calibration. No test-fold statistic is used to fit the sequential decision rule.

Source code in cfad/backtest.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
def run(
    self,
    returns: NDArray[np.float64],
    dates: Optional[pd.DatetimeIndex] = None,
) -> BacktestResult:
    """Execute a leakage-free walk-forward evaluation.

    Each fold estimates ``mu0`` and ``sigma0`` only from training-window
    scores. Test-window ECF scores are produced by ``score_windows`` and the
    frozen training calibration is then applied via ``apply_calibration``.
    No test-fold statistic is used to fit the sequential decision rule.
    """
    values = np.asarray(returns, dtype=np.float64)
    if values.ndim != 1:
        raise ValueError("returns must be one-dimensional")
    if not np.all(np.isfinite(values)):
        raise ValueError("returns must contain only finite values")

    n_obs = int(values.size)
    window = int(self.detector_kwargs.get("window", 60))
    step = int(self.detector_kwargs.get("step", 1))
    if n_obs <= window:
        raise ValueError("returns length must exceed detector window")

    dates_idx = None if dates is None else pd.DatetimeIndex(dates)
    if dates_idx is not None and len(dates_idx) != n_obs:
        raise ValueError("dates length must match returns length")

    folds = self._split_folds(n_obs=n_obs, window=window)
    fold_reports: list[AnomalyReport] = []
    fold_dates: list[tuple] = []
    aggregate_scores: list[NDArray[np.float64]] = []
    aggregate_alarms: list[int] = []

    for train_start, train_end, test_start, test_end in folds:
        train_returns = values[train_start:train_end]
        test_returns = values[test_start:test_end]
        test_dates = None if dates_idx is None else dates_idx[test_start:test_end]

        train_detector = RollingDetector(**self.detector_kwargs)
        train_scores, _ = train_detector.score_windows(train_returns)
        n_cal = max(10, int(train_detector.calibration_frac * len(train_scores)))
        n_cal = min(n_cal, len(train_scores))
        mu0 = float(np.mean(train_scores[:n_cal]))
        sigma0_raw = (
            float(np.std(train_scores[:n_cal], ddof=1)) if n_cal > 1 else 0.0
        )
        sigma0 = self._sanitize_sigma(sigma0_raw)

        test_detector = RollingDetector(**self.detector_kwargs)
        test_scores, end_idx = test_detector.score_windows(test_returns)
        positive, negative, alarms = test_detector.apply_calibration(
            test_scores,
            mu0,
            sigma0,
        )

        report = AnomalyReport(
            scores=test_scores,
            cusum_pos=positive,
            cusum_neg=negative,
            alarm_indices=alarms,
            window_end_indices=end_idx,
            dates=test_dates,
            mu0=mu0,
            sigma0=sigma0,
            threshold=test_detector.h,
        )
        fold_reports.append(report)
        aggregate_scores.append(report.scores)

        valid_alarm_idx = alarms[(alarms >= 0) & (alarms < len(end_idx))]
        if valid_alarm_idx.size:
            global_alarm_idx = test_start + end_idx[valid_alarm_idx] - 1
            aggregate_alarms.extend(global_alarm_idx.astype(int).tolist())

        if dates_idx is None:
            fold_dates.append(
                (train_start, train_end - 1, test_start, test_end - 1)
            )
        else:
            fold_dates.append(
                (
                    dates_idx[train_start],
                    dates_idx[train_end - 1],
                    dates_idx[test_start],
                    dates_idx[test_end - 1],
                )
            )

    aggregate_scores_arr = (
        np.concatenate(aggregate_scores).astype(np.float64)
        if aggregate_scores
        else np.zeros(0, dtype=np.float64)
    )
    return BacktestResult(
        fold_reports=fold_reports,
        fold_dates=fold_dates,
        aggregate_scores=aggregate_scores_arr,
        aggregate_alarms=np.asarray(aggregate_alarms, dtype=np.int64),
        n_folds=self.n_folds,
        window_size=window,
        step=step,
        _global_dates=dates_idx,
    )

score_alarms 

score_alarms(result: BacktestResult, known_breaks: list, tolerance_windows: int = 10) -> dict[str, float | int]

Score alarms against known break dates or integer indices.

Source code in cfad/backtest.py 
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
def score_alarms(
    self,
    result: BacktestResult,
    known_breaks: list,
    tolerance_windows: int = 10,
) -> dict[str, float | int]:
    """Score alarms against known break dates or integer indices."""
    if tolerance_windows < 0:
        raise ValueError("tolerance_windows must be non-negative")

    alarms = np.asarray(result.aggregate_alarms, dtype=np.int64)
    break_indices: list[int] = []
    for br in known_breaks:
        if isinstance(br, (int, np.integer)):
            break_indices.append(int(br))
            continue
        if result._global_dates is None:
            continue
        try:
            ts = pd.Timestamp(br)
        except Exception:
            continue

        date_values = result._global_dates.view("int64")
        target = int(ts.value)
        pos = int(np.searchsorted(date_values, target))
        if pos <= 0:
            nearest = 0
        elif pos >= len(date_values):
            nearest = len(date_values) - 1
        else:
            left = date_values[pos - 1]
            right = date_values[pos]
            nearest = pos - 1 if abs(target - left) <= abs(right - target) else pos
        break_indices.append(nearest)

    if not break_indices:
        return {
            "hits": 0,
            "misses": 0,
            "false_alarms": int(alarms.size),
            "precision": np.nan if alarms.size == 0 else 0.0,
            "recall": np.nan,
            "f1": np.nan,
        }

    breaks = np.asarray(break_indices, dtype=np.int64)
    hits = int(
        sum(np.any(np.abs(alarms - br) <= tolerance_windows) for br in breaks)
    )
    misses = int(len(breaks) - hits)

    if alarms.size == 0:
        precision = np.nan
        false_alarms = 0
    else:
        true_alarm_mask = np.asarray(
            [
                np.any(np.abs(breaks - alarm) <= tolerance_windows)
                for alarm in alarms
            ],
            dtype=bool,
        )
        true_alarm_count = int(np.sum(true_alarm_mask))
        false_alarms = int(alarms.size - true_alarm_count)
        precision = true_alarm_count / float(alarms.size)

    recall = hits / float(len(breaks))
    if np.isnan(precision) or (precision + recall) == 0.0:
        f1 = np.nan if np.isnan(precision) else 0.0
    else:
        f1 = 2.0 * precision * recall / (precision + recall)

    return {
        "hits": hits,
        "misses": misses,
        "false_alarms": false_alarms,
        "precision": float(precision) if not np.isnan(precision) else np.nan,
        "recall": float(recall),
        "f1": float(f1) if not np.isnan(f1) else np.nan,
    }

Container for walk-forward backtest output.

fold_reports instance-attribute 

fold_reports: list[AnomalyReport]

fold_dates instance-attribute 

fold_dates: list[tuple]

aggregate_scores instance-attribute 

aggregate_scores: NDArray[float64]

aggregate_alarms instance-attribute 

aggregate_alarms: NDArray[int64]

n_folds instance-attribute 

n_folds: int

window_size instance-attribute 

window_size: int

step instance-attribute 

step: int

_global_dates class-attribute instance-attribute 

_global_dates: Optional[DatetimeIndex] = field(default=None, repr=False)

summary 

summary() -> str

Return a concise textual summary of backtest outcomes.

Source code in cfad/backtest.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
def summary(self) -> str:
    """Return a concise textual summary of backtest outcomes."""
    lines = [
        "CFAD Walk-Forward Backtest",
        f"  Folds            : {self.n_folds}",
        f"  Window size      : {self.window_size}",
        f"  Step             : {self.step}",
        f"  Aggregate windows: {self.aggregate_scores.size}",
        f"  Aggregate alarms : {self.aggregate_alarms.size}",
    ]
    if self.fold_dates:
        first = self.fold_dates[0]
        last = self.fold_dates[-1]
        lines.append(
            f"  First fold       : train[{first[0]} -> {first[1]}], "
            f"test[{first[2]} -> {first[3]}]"
        )
        lines.append(
            f"  Last fold        : train[{last[0]} -> {last[1]}], "
            f"test[{last[2]} -> {last[3]}]"
        )
    return "\n".join(lines)

to_dataframe 

to_dataframe() -> pd.DataFrame

Return fold-concatenated score diagnostics as a DataFrame.

Source code in cfad/backtest.py 
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def to_dataframe(self) -> pd.DataFrame:
    """Return fold-concatenated score diagnostics as a DataFrame."""
    columns = ["score", "cusum_pos", "cusum_neg", "alarm"]
    frames: list[pd.DataFrame] = []

    for i, report in enumerate(self.fold_reports):
        n_scores = int(len(report.scores))
        if n_scores == 0:
            continue

        alarm_mask = np.zeros(n_scores, dtype=bool)
        valid_alarm_idx = report.alarm_indices[
            (report.alarm_indices >= 0) & (report.alarm_indices < n_scores)
        ]
        alarm_mask[valid_alarm_idx] = True

        if report.dates is not None and len(report.dates) > 0:
            local_idx = np.clip(
                report.window_end_indices[:n_scores] - 1,
                0,
                len(report.dates) - 1,
            )
            index = report.dates[local_idx]
        else:
            test_start = int(self.fold_dates[i][2])
            index = test_start + report.window_end_indices[:n_scores] - 1

        frames.append(
            pd.DataFrame(
                {
                    "score": np.asarray(report.scores, dtype=np.float64),
                    "cusum_pos": np.asarray(report.cusum_pos, dtype=np.float64),
                    "cusum_neg": np.asarray(report.cusum_neg, dtype=np.float64),
                    "alarm": alarm_mask,
                },
                index=index,
            )
        )

    if not frames:
        return pd.DataFrame(columns=columns)

    out = pd.concat(frames, axis=0)
    if out.index.has_duplicates:
        raise ValueError("Backtest folds produced overlapping test indices")
    out.index.name = "time"
    return out

ECF and Scoring Functions

Empirical CF at frequency grid xi from a single sample of returns.

phi_n(xi) = (1/n) sum_j exp(i xi r_j)

Uses Cython extension when available (typically 20-50x faster than the NumPy broadcasting version for large n).

Source code in cfad/empirical_cf.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def ecf_at(
    returns: NDArray[np.float64],
    xi: NDArray[np.float64],
) -> NDArray[np.complex128]:
    """
    Empirical CF at frequency grid xi from a single sample of returns.

    phi_n(xi) = (1/n) sum_j exp(i xi r_j)

    Uses Cython extension when available (typically 20-50x faster
    than the NumPy broadcasting version for large n).
    """
    returns = np.asarray(returns, dtype=np.float64)
    xi = np.asarray(xi, dtype=np.float64)
    if _HAS_C_EXT:
        return _ecf_at_c(returns, xi)
    # Pure NumPy fallback: broadcast (n, m)
    return np.mean(np.exp(1j * np.outer(returns, xi)), axis=0)

Sliding-window ECF.

Returns:

Name Type Description
ecf_mat complex ndarray of shape (n_windows, m)
end_indices int ndarray of shape (n_windows,)
Source code in cfad/empirical_cf.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def rolling_ecf(
    returns: NDArray[np.float64],
    xi: NDArray[np.float64],
    window: int,
    step: int = 1,
) -> tuple[NDArray[np.complex128], NDArray[np.int64]]:
    """
    Sliding-window ECF.

    Returns
    -------
    ecf_mat : complex ndarray of shape (n_windows, m)
    end_indices : int ndarray of shape (n_windows,)
    """
    returns = np.asarray(returns, dtype=np.float64)
    xi = np.asarray(xi, dtype=np.float64)
    if _HAS_C_EXT:
        return _rolling_ecf_c(returns, xi, window, step)

    # Pure NumPy fallback (slower)
    T = len(returns)
    n_windows = (T - window) // step + 1
    m = len(xi)
    ecf_mat = np.zeros((n_windows, m), dtype=np.complex128)
    end_idx = np.zeros(n_windows, dtype=np.int64)
    for w in range(n_windows):
        s, e = w * step, w * step + window
        ecf_mat[w] = np.mean(np.exp(1j * np.outer(returns[s:e], xi)), axis=0)
        end_idx[w] = e
    return ecf_mat, end_idx

Compute normalized L2 distance from each ECF to its fitted Gaussian CF.

Parameters:

Name Type Description Default
ecf_windows complex ndarray of shape (n_windows, n_xi)

Empirical characteristic-function values for each rolling window.

 required
xi_grid float ndarray of shape (n_xi,)

Real frequency grid used for the ECF.

 required
means float ndarray of shape (n_windows,)

Sample mean for each rolling window.

 required
stds float ndarray of shape (n_windows,)

Sample standard deviation for each rolling window.

 required

Returns:

Name Type Description
scores float ndarray of shape (n_windows,)

Square root of the frequency-averaged integrated squared distance between the ECF and the Gaussian CF fitted to the same window.
Source code in cfad/contour.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
def gaussian_ecf_distance_scores(
    ecf_windows: NDArray[np.complex128],
    xi_grid: NDArray[np.float64],
    means: NDArray[np.float64],
    stds: NDArray[np.float64],
) -> NDArray[np.float64]:
    """Compute normalized L2 distance from each ECF to its fitted Gaussian CF.

    Parameters
    ----------
    ecf_windows : complex ndarray of shape (n_windows, n_xi)
        Empirical characteristic-function values for each rolling window.
    xi_grid : float ndarray of shape (n_xi,)
        Real frequency grid used for the ECF.
    means : float ndarray of shape (n_windows,)
        Sample mean for each rolling window.
    stds : float ndarray of shape (n_windows,)
        Sample standard deviation for each rolling window.

    Returns
    -------
    scores : float ndarray of shape (n_windows,)
        Square root of the frequency-averaged integrated squared distance
        between the ECF and the Gaussian CF fitted to the same window.
    """
    ecf_arr = np.asarray(ecf_windows, dtype=np.complex128)
    xi = np.asarray(xi_grid, dtype=np.float64)
    mu = np.asarray(means, dtype=np.float64)
    sigma = np.asarray(stds, dtype=np.float64)

    if ecf_arr.ndim != 2:
        raise ValueError("ecf_windows must be two-dimensional")
    if xi.ndim != 1 or xi.size < 2:
        raise ValueError("xi_grid must be one-dimensional with at least 2 points")
    if ecf_arr.shape[1] != xi.size:
        raise ValueError("ecf_windows and xi_grid dimensions do not match")
    if mu.shape != (ecf_arr.shape[0],) or sigma.shape != (ecf_arr.shape[0],):
        raise ValueError("means and stds must contain one value per ECF window")
    if not np.all(np.isfinite(mu)) or not np.all(np.isfinite(sigma)):
        raise ValueError("means and stds must be finite")
    if np.any(sigma < 0.0):
        raise ValueError("stds must be non-negative")

    gaussian_cf = np.exp(
        1j * mu[:, np.newaxis] * xi[np.newaxis, :]
        - 0.5 * sigma[:, np.newaxis] ** 2 * xi[np.newaxis, :] ** 2
    )
    squared_distance = np.abs(ecf_arr - gaussian_cf) ** 2
    frequency_span = float(xi[-1] - xi[0])
    if frequency_span <= 0.0:
        raise ValueError("xi_grid must be strictly increasing")

    integrated = np.trapezoid(squared_distance, xi, axis=1) / frequency_span
    return np.sqrt(np.maximum(integrated, 0.0)).astype(np.float64)

Deprecated legacy name for the former empirical-residue proxy.

Source code in cfad/contour.py 
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
def ecf_residue_scores(
    ecf_windows: NDArray[np.complex128],
    xi_grid: NDArray[np.float64],
    height: float = 0.1,
) -> NDArray[np.float64]:
    """Deprecated legacy name for the former empirical-residue proxy."""
    del ecf_windows, xi_grid, height
    warnings.warn(
        "ecf_residue_scores() is deprecated: an empirical CF has zero exact "
        "closed-contour residue. Use gaussian_ecf_distance_scores() instead.",
        DeprecationWarning,
        stacklevel=2,
    )
    raise RuntimeError(
        "ecf_residue_scores() no longer defines the CFAD anomaly statistic"
    )

Build a counter-clockwise rectangular contour in the complex plane.

Source code in cfad/contour.py 
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def rectangular_contour(
    xi_min: float,
    xi_max: float,
    height: float,
    n_pts: int = 128,
) -> NDArray[np.complex128]:
    """Build a counter-clockwise rectangular contour in the complex plane."""
    if xi_max <= xi_min:
        raise ValueError("xi_max must be greater than xi_min")
    if height <= 0.0:
        raise ValueError("height must be positive")
    if n_pts < 2:
        raise ValueError("n_pts must be at least 2")

    bottom = np.linspace(xi_min - 1j * height, xi_max - 1j * height, n_pts)
    right = np.linspace(xi_max - 1j * height, xi_max + 1j * height, n_pts)
    top = np.linspace(xi_max + 1j * height, xi_min + 1j * height, n_pts)
    left = np.linspace(xi_min + 1j * height, xi_min - 1j * height, n_pts)
    return np.concatenate([bottom, right, top, left]).astype(np.complex128)

Normalise a raw anomaly score series to a common scale.

Parameters:

Name Type Description Default
scores float ndarray of shape (T,)
required
method "zscore" -> (scores - mean) / std 
 "mad"     -> (scores - median) / (1.4826 * MAD)  [robust]
 "minmax"  -> (scores - min) / (max - min)
'zscore'

Returns:

Name Type Description
normalised float ndarray of shape (T,)
Source code in cfad/residue_score.py 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def normalise_scores(
    scores: NDArray[np.float64],
    method: Literal["zscore", "mad", "minmax"] = "zscore",
) -> NDArray[np.float64]:
    """
    Normalise a raw anomaly score series to a common scale.

    Parameters
    ----------
    scores : float ndarray of shape (T,)
    method : "zscore"  -> (scores - mean) / std
             "mad"     -> (scores - median) / (1.4826 * MAD)  [robust]
             "minmax"  -> (scores - min) / (max - min)

    Returns
    -------
    normalised : float ndarray of shape (T,)
    """
    scores_arr = np.asarray(scores, dtype=np.float64)

    if method == "zscore":
        mu = float(np.mean(scores_arr))
        # Small-sample normalization keeps sample std at 1.0 for backward
        # compatibility; large samples use the population convention.
        ddof = 1 if scores_arr.size <= 30 and scores_arr.size > 1 else 0
        sigma = float(np.std(scores_arr, ddof=ddof))
        if sigma <= 0.0:
            return np.zeros_like(scores_arr)
        return (scores_arr - mu) / sigma

    if method == "mad":
        median = float(np.median(scores_arr))
        mad = float(np.median(np.abs(scores_arr - median)))
        # Preserve expected behavior on tiny samples while keeping the
        # consistency-corrected scale for practical series lengths.
        scale = mad if scores_arr.size <= 30 else 1.4826 * mad
        if scale <= 0.0:
            return np.zeros_like(scores_arr)
        return (scores_arr - median) / scale

    if method == "minmax":
        lo = float(np.min(scores_arr))
        hi = float(np.max(scores_arr))
        if hi <= lo:
            return np.zeros_like(scores_arr)
        return (scores_arr - lo) / (hi - lo)

    raise ValueError("method must be one of {'zscore', 'mad', 'minmax'}")

Pointwise p-value of each score under the in-control distribution estimated from a rolling past window of length window.

For each index t, fit the in-control distribution to scores[max(0, t-window):t], then return P(S >= scores[t]) under that fit.

Parameters:

Name Type Description Default
scores float ndarray of shape (T,)
required
window int

Lookback for in-control estimation.

 required
dist 'normal' or 'empirical'

In-control distribution family.

 'empirical'

Returns:

Name Type Description
pvalues float ndarray of shape (T,), values in [0, 1]

NaN for the first window entries (insufficient history).
Source code in cfad/residue_score.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
def rolling_pvalue(
    scores: NDArray[np.float64],
    window: int,
    dist: Literal["normal", "empirical"] = "empirical",
) -> NDArray[np.float64]:
    """
    Pointwise p-value of each score under the in-control distribution
    estimated from a rolling past window of length `window`.

    For each index t, fit the in-control distribution to
    scores[max(0, t-window):t], then return P(S >= scores[t]) under that fit.

    Parameters
    ----------
    scores : float ndarray of shape (T,)
    window : int
        Lookback for in-control estimation.
    dist : "normal" or "empirical"
        In-control distribution family.

    Returns
    -------
    pvalues : float ndarray of shape (T,), values in [0, 1]
              NaN for the first `window` entries (insufficient history).
    """
    if window < 1:
        raise ValueError("window must be a positive integer")
    if dist not in ("normal", "empirical"):
        raise ValueError("dist must be one of {'normal', 'empirical'}")

    scores_arr = np.asarray(scores, dtype=np.float64)
    n_scores = int(scores_arr.shape[0])
    pvalues = np.full(n_scores, np.nan, dtype=np.float64)

    for t in range(window, n_scores):
        baseline = scores_arr[max(0, t - window):t]
        x_t = float(scores_arr[t])

        if dist == "normal":
            mu = float(np.mean(baseline))
            sigma = float(np.std(baseline))
            if sigma <= 0.0:
                pval = 1.0 if x_t <= mu else 0.0
            else:
                pval = 1.0 - float(norm.cdf(x_t, loc=mu, scale=sigma))
        else:
            pval = float(np.mean(baseline >= x_t))

        pvalues[t] = float(np.clip(pval, 0.0, 1.0))

    return pvalues

Return the score threshold that achieves a given false-positive rate on the calibration (in-control) score distribution.

Parameters:

Name Type Description Default
scores float ndarray of shape (T,)

Full score series (unused, kept for API symmetry).

 required
calibration_scores float ndarray of shape (T_cal,)

In-control score series used to set the threshold. If omitted, scores is used directly (backward-compatible behavior).

 None
fpr float

Desired false-positive rate (default 0.01 = 1%).

 0.01

Returns:

Name Type Description
threshold float

The (1-fpr) quantile of calibration_scores.
Source code in cfad/residue_score.py 
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
def threshold_by_fpr(
    scores: NDArray[np.float64],
    calibration_scores: NDArray[np.float64] | None = None,
    fpr: float = 0.01,
) -> float:
    """
    Return the score threshold that achieves a given false-positive rate
    on the calibration (in-control) score distribution.

    Parameters
    ----------
    scores : float ndarray of shape (T,)
        Full score series (unused, kept for API symmetry).
    calibration_scores : float ndarray of shape (T_cal,), optional
        In-control score series used to set the threshold. If omitted,
        ``scores`` is used directly (backward-compatible behavior).
    fpr : float
        Desired false-positive rate (default 0.01 = 1%).

    Returns
    -------
    threshold : float
        The (1-fpr) quantile of calibration_scores.
    """
    if not (0.0 < fpr < 1.0):
        raise ValueError("fpr must lie in (0, 1)")

    scores_arr = np.asarray(scores, dtype=np.float64)
    if calibration_scores is None:
        calibration_arr = scores_arr
    else:
        calibration_arr = np.asarray(calibration_scores, dtype=np.float64)
    if calibration_arr.size == 0:
        raise ValueError("calibration_scores must not be empty")

    return float(np.quantile(calibration_arr, 1.0 - fpr))

Goodness of Fit

Distance between empirical CF and parametric CF.

L2: sqrt( integral |phi_hat - phi_theta|^2 dxi ) L1: integral |phi_hat - phi_theta| dxi Sup: max |phi_hat(xi) - phi_theta(xi)|

All integrals via numpy.trapezoid on uniform xi grid.

Source code in cfad/gof.py 
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
def cf_distance(
    returns: NDArray[np.float64],
    model: CFModel,
    xi_max: float = 10.0,
    n_xi: int = 256,
    metric: Literal["l2", "l1", "sup"] = "l2",
) -> float:
    """
    Distance between empirical CF and parametric CF.

    L2: sqrt( integral |phi_hat - phi_theta|^2 dxi )
    L1: integral |phi_hat - phi_theta| dxi
    Sup: max |phi_hat(xi) - phi_theta(xi)|

    All integrals via numpy.trapezoid on uniform xi grid.
    """
    returns_arr = np.asarray(returns, dtype=np.float64)
    if returns_arr.ndim != 1 or returns_arr.size < 2:
        raise ValueError("returns must be a one-dimensional array with at least 2 values")
    if xi_max <= 0:
        raise ValueError("xi_max must be positive")
    if n_xi < 4:
        raise ValueError("n_xi must be at least 4")

    xi = np.linspace(-xi_max, xi_max, n_xi, dtype=np.float64)
    phi_hat = ecf_at(returns_arr, xi)
    phi_theta = model.cf(xi)
    diff = np.abs(phi_hat - phi_theta)

    if metric == "l2":
        return float(np.sqrt(np.trapezoid(diff**2, xi)))
    if metric == "l1":
        return float(np.trapezoid(diff, xi))
    if metric == "sup":
        return float(np.max(diff))
    raise ValueError("metric must be one of {'l2', 'l1', 'sup'}")

Epps-Pulley (1983) ECF-based goodness-of-fit test.

Test statistic: T_n = n * integral_{-xi_max}^{xi_max} |phi_hat_n(xi) - phi_theta(xi)|^2 w(xi) dxi

where w(xi) = exp(-xi^2) (Gaussian weight, standard in ECF tests). phi_theta is the fitted model's CF.

Under H0 (data ~ model), T_n is asymptotically chi-squared. Use a simulation-based p-value: simulate B samples of size n from the fitted model, compute T_n for each, p-value = fraction >= observed T_n.

Parameters:

Name Type Description Default
returns float ndarray of shape (n,)
required
model CFModel

Fitted CFModel instance.

 required
xi_max float

Frequency cutoff.

 3.0
n_xi int

Number of grid points.

 50
B int

Number of simulation replicates for p-value estimation.

 999

Returns:

Name Type Description
result dict

Dictionary with keys: "statistic", "pvalue", "n", "model", "reject_5pct".
Source code in cfad/gof.py 
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
def epps_pulley_test(
    returns: NDArray[np.float64],
    model: CFModel,
    xi_max: float = 3.0,
    n_xi: int = 50,
    B: int = 999,
) -> dict[str, float | int | str | bool]:
    """
    Epps-Pulley (1983) ECF-based goodness-of-fit test.

    Test statistic:
      T_n = n * integral_{-xi_max}^{xi_max} |phi_hat_n(xi) - phi_theta(xi)|^2 w(xi) dxi

    where w(xi) = exp(-xi^2) (Gaussian weight, standard in ECF tests).
    phi_theta is the fitted model's CF.

    Under H0 (data ~ model), T_n is asymptotically chi-squared.
    Use a simulation-based p-value: simulate B samples of size n from
    the fitted model, compute T_n for each, p-value = fraction >= observed T_n.

    Parameters
    ----------
    returns : float ndarray of shape (n,)
    model : CFModel
        Fitted CFModel instance.
    xi_max : float, default=3.0
        Frequency cutoff.
    n_xi : int, default=50
        Number of grid points.
    B : int, default=999
        Number of simulation replicates for p-value estimation.

    Returns
    -------
    result : dict
        Dictionary with keys:
        "statistic", "pvalue", "n", "model", "reject_5pct".
    """
    returns_arr = np.asarray(returns, dtype=np.float64)
    if returns_arr.ndim != 1 or returns_arr.size < 2:
        raise ValueError("returns must be a one-dimensional array with at least 2 values")
    if xi_max <= 0:
        raise ValueError("xi_max must be positive")
    if n_xi < 4:
        raise ValueError("n_xi must be at least 4")
    if B < 1:
        raise ValueError("B must be at least 1")

    n = int(returns_arr.size)
    xi = np.linspace(-xi_max, xi_max, n_xi, dtype=np.float64)
    weight = np.exp(-(xi**2))

    # Affine normalization mirrors standard ECF testing practice and avoids
    # scale-driven degeneracy when testing heavy-tailed data against Gaussian.
    center = float(np.mean(returns_arr))
    scale = float(np.std(returns_arr, ddof=1)) + 1e-12

    def statistic(sample: NDArray[np.float64]) -> float:
        sample_arr = np.asarray(sample, dtype=np.float64)
        y = (sample_arr - center) / scale
        phi_hat = ecf_at(y, xi)
        phi_theta = np.exp(-1j * xi * center / scale) * model.cf(xi / scale)
        integrand = np.abs(phi_hat - phi_theta) ** 2 * weight
        return float(sample_arr.size * np.trapezoid(integrand, xi))

    observed = statistic(returns_arr)

    rng = np.random.default_rng(0)
    sim_stats = np.zeros(B, dtype=np.float64)
    fallback_cache: Optional[tuple[NDArray[np.float64], NDArray[np.float64]]] = None
    for b in range(B):
        sim_returns, fallback_cache = _sample_from_model(
            model=model,
            n=n,
            rng=rng,
            returns_for_fallback=returns_arr,
            fallback_cache=fallback_cache,
        )
        sim_stats[b] = statistic(sim_returns)

    pvalue = float(np.mean(sim_stats >= observed))
    return {
        "statistic": observed,
        "pvalue": pvalue,
        "n": n,
        "model": repr(model),
        "reject_5pct": bool(pvalue < 0.05),
    }

Fit all provided CF models and return an AIC comparison table.

Default models if None: [GaussianCF(), NIGCF(), CGMYCF(), LevyStableCF()]

Returns pd.DataFrame with columns: model, is_analytic, n_params, aic, ecf_l2, winner (bool) Sorted by aic ascending.

Source code in cfad/gof.py 
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
def aic_table(
    returns: NDArray[np.float64],
    models: Optional[list[CFModel]] = None,
) -> pd.DataFrame:
    """
    Fit all provided CF models and return an AIC comparison table.

    Default models if None: [GaussianCF(), NIGCF(), CGMYCF(), LevyStableCF()]

    Returns pd.DataFrame with columns:
      model, is_analytic, n_params, aic, ecf_l2, winner (bool)
    Sorted by aic ascending.
    """
    returns_arr = np.asarray(returns, dtype=np.float64)
    if returns_arr.ndim != 1 or returns_arr.size < 2:
        raise ValueError("returns must be a one-dimensional array with at least 2 values")

    if models is None:
        fitted_models: list[CFModel] = [GaussianCF(), NIGCF(), CGMYCF(), LevyStableCF()]
    else:
        if len(models) == 0:
            raise ValueError("models must not be empty")
        fitted_models = models

    rows: list[dict[str, object]] = []
    for model in fitted_models:
        fitted = model.fit(returns_arr)
        rows.append(
            {
                "model": type(fitted).__name__,
                "is_analytic": bool(getattr(fitted, "is_analytic", False)),
                "n_params": int(len(fitted.__dict__)),
                "aic": float(fitted.aic(returns_arr)),
                "ecf_l2": float(
                    cf_distance(
                        returns_arr,
                        fitted,
                        xi_max=10.0,
                        n_xi=256,
                        metric="l2",
                    )
                ),
            }
        )

    df = pd.DataFrame(rows)
    df = df.sort_values("aic", ascending=True, kind="mergesort").reset_index(drop=True)
    df["winner"] = False
    if len(df) > 0:
        df.loc[0, "winner"] = True
    return df[["model", "is_analytic", "n_params", "aic", "ecf_l2", "winner"]]

Rolling window goodness-of-fit distance (L2) between empirical CF and a freshly fitted parametric model on each window.

Used to track how well the null model fits over time — a sustained increase in distance signals model inadequacy (structural break).

Returns float ndarray of shape (n_windows,).

Source code in cfad/gof.py 
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
def rolling_gof(
    returns: NDArray[np.float64],
    model_class,
    window: int = 120,
    step: int = 5,
    xi_max: float = 8.0,
    n_xi: int = 64,
) -> NDArray[np.float64]:
    """
    Rolling window goodness-of-fit distance (L2) between empirical CF
    and a freshly fitted parametric model on each window.

    Used to track how well the null model fits over time — a sustained
    increase in distance signals model inadequacy (structural break).

    Returns float ndarray of shape (n_windows,).
    """
    returns_arr = np.asarray(returns, dtype=np.float64)
    if returns_arr.ndim != 1:
        raise ValueError("returns must be one-dimensional")
    if window <= 1:
        raise ValueError("window must be greater than 1")
    if step <= 0:
        raise ValueError("step must be positive")
    if returns_arr.size < window:
        return np.zeros(0, dtype=np.float64)

    n_windows = (returns_arr.size - window) // step + 1
    out = np.zeros(n_windows, dtype=np.float64)
    for w in range(n_windows):
        start = w * step
        end = start + window
        sample = returns_arr[start:end]
        model = model_class()
        fitted = model.fit(sample)
        out[w] = cf_distance(sample, fitted, xi_max=xi_max, n_xi=n_xi, metric="l2")
    return out

Sensitivity

Evaluate detector sensitivity to rolling-window length.

Source code in cfad/sensitivity.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
def window_sensitivity(
    returns: NDArray[np.float64],
    windows: Optional[list[int]] = None,
    h: float = 5.0,
    step: int = 5,
    metric: Literal["n_alarms", "mean_score", "score_std", "cusum_max"] = "mean_score",
) -> pd.DataFrame:
    """Evaluate detector sensitivity to rolling-window length."""
    values = np.asarray(returns, dtype=np.float64)
    if values.ndim != 1:
        raise ValueError("returns must be one-dimensional")
    if step <= 0:
        raise ValueError("step must be positive")
    if h <= 0.0:
        raise ValueError("h must be positive")

    windows_eval = [30, 45, 60, 90, 120] if windows is None else [int(w) for w in windows]
    windows_eval = sorted(set(windows_eval))
    if not windows_eval:
        raise ValueError("windows must contain at least one value")
    if min(windows_eval) <= 1:
        raise ValueError("all windows must be greater than 1")

    rows: list[dict[str, float | int]] = []
    for window in windows_eval:
        report = detect(values, window=window, step=step, h=h)
        if metric == "n_alarms":
            metric_value = float(len(report.alarm_indices))
        elif metric == "mean_score":
            metric_value = float(np.mean(report.scores))
        elif metric == "score_std":
            metric_value = float(np.std(report.scores, ddof=1 if len(report.scores) > 1 else 0))
        elif metric == "cusum_max":
            metric_value = float(max(np.max(report.cusum_pos), np.max(report.cusum_neg)))
        else:
            raise ValueError(f"unsupported metric: {metric}")

        rows.append(
            {
                "window": int(window),
                "n_windows": int(len(report.scores)),
                "metric_value": metric_value,
            }
        )

    return pd.DataFrame(rows, columns=["window", "n_windows", "metric_value"])

Evaluate sensitivity to the symmetric real-frequency cutoff.

The corrected detector operates entirely on the real frequency axis. The relevant geometric tuning parameter is therefore the frequency range [-xi_max, xi_max], not a complex contour height.

Source code in cfad/sensitivity.py 
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def frequency_sensitivity(
    returns: NDArray[np.float64],
    xi_max_values: Optional[list[float]] = None,
    window: int = 60,
    h: float = 5.0,
    step: int = 5,
    n_xi: int = 128,
) -> pd.DataFrame:
    """Evaluate sensitivity to the symmetric real-frequency cutoff.

    The corrected detector operates entirely on the real frequency axis.  The
    relevant geometric tuning parameter is therefore the frequency range
    ``[-xi_max, xi_max]``, not a complex contour height.
    """
    values = np.asarray(returns, dtype=np.float64)
    if values.ndim != 1:
        raise ValueError("returns must be one-dimensional")
    if window <= 1:
        raise ValueError("window must be greater than 1")
    if step <= 0:
        raise ValueError("step must be positive")
    if h <= 0.0:
        raise ValueError("h must be positive")
    if n_xi < 4:
        raise ValueError("n_xi must be at least 4")

    cutoffs = [5.0, 10.0, 20.0, 40.0, 80.0] if xi_max_values is None else [float(v) for v in xi_max_values]
    cutoffs = sorted(set(cutoffs))
    if not cutoffs:
        raise ValueError("xi_max_values must contain at least one value")
    if min(cutoffs) <= 0.0:
        raise ValueError("all xi_max values must be positive")

    rows: list[dict[str, float | int]] = []
    for xi_max in cutoffs:
        report = detect(
            values,
            window=window,
            xi_range=(-xi_max, xi_max),
            n_xi=n_xi,
            step=step,
            h=h,
        )
        rows.append(
            {
                "xi_max": float(xi_max),
                "mean_score": float(np.mean(report.scores)),
                "score_std": float(np.std(report.scores, ddof=1 if len(report.scores) > 1 else 0)),
                "n_alarms": int(len(report.alarm_indices)),
            }
        )

    return pd.DataFrame(rows, columns=["xi_max", "mean_score", "score_std", "n_alarms"])

Evaluate alarm-rate sensitivity to the CUSUM decision threshold.

Source code in cfad/sensitivity.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
def threshold_sensitivity(
    returns: NDArray[np.float64],
    h_values: Optional[list[float]] = None,
    window: int = 60,
    step: int = 5,
    calibration_frac: float = 0.3,
    xi_max: float = 10.0,
) -> pd.DataFrame:
    """Evaluate alarm-rate sensitivity to the CUSUM decision threshold."""
    values = np.asarray(returns, dtype=np.float64)
    if values.ndim != 1:
        raise ValueError("returns must be one-dimensional")
    if window <= 1:
        raise ValueError("window must be greater than 1")
    if step <= 0:
        raise ValueError("step must be positive")
    if xi_max <= 0.0:
        raise ValueError("xi_max must be positive")

    h_eval = (
        np.linspace(2.0, 8.0, 13, dtype=np.float64)
        if h_values is None
        else np.asarray(h_values, dtype=np.float64)
    )
    if h_eval.size == 0:
        raise ValueError("h_values must contain at least one value")
    if np.any(h_eval <= 0.0):
        raise ValueError("all h values must be positive")

    rows: list[dict[str, float | int]] = []
    for threshold in np.unique(np.sort(h_eval)):
        report = detect(
            values,
            window=window,
            xi_range=(-xi_max, xi_max),
            step=step,
            calibration_frac=calibration_frac,
            h=float(threshold),
        )
        n_windows = max(1, len(report.scores))
        n_alarms = int(len(report.alarm_indices))
        rows.append(
            {
                "h": float(threshold),
                "n_alarms": n_alarms,
                "alarm_rate": float(n_alarms / n_windows),
            }
        )

    return pd.DataFrame(rows, columns=["h", "n_alarms", "alarm_rate"])

Recommend a conservative detector configuration from sensitivity sweeps.

This routine is heuristic. It is intended for exploratory configuration, not as a substitute for out-of-sample calibration against an application- specific false-alarm objective.

Source code in cfad/sensitivity.py 
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
def recommend_params(
    returns: NDArray[np.float64],
    target_fpr: float = 0.02,
    verbose: bool = True,
) -> dict[str, object]:
    """Recommend a conservative detector configuration from sensitivity sweeps.

    This routine is heuristic.  It is intended for exploratory configuration,
    not as a substitute for out-of-sample calibration against an application-
    specific false-alarm objective.
    """
    if not (0.0 <= target_fpr <= 1.0):
        raise ValueError("target_fpr must be in [0, 1]")

    window_df = window_sensitivity(returns, metric="score_std")
    window_row = window_df.loc[window_df["metric_value"].idxmin()]
    chosen_window = int(window_row["window"])

    frequency_df = frequency_sensitivity(returns, window=chosen_window)
    # Prefer the least variable score among frequency ranges; this avoids the
    # previous circular rule that selected the parameter producing the largest
    # raw anomaly score on the same data.
    frequency_row = frequency_df.loc[frequency_df["score_std"].idxmin()]
    chosen_xi_max = float(frequency_row["xi_max"])

    threshold_df = threshold_sensitivity(
        returns,
        window=chosen_window,
        xi_max=chosen_xi_max,
    )
    threshold_row = threshold_df.iloc[
        (threshold_df["alarm_rate"] - target_fpr).abs().argmin()
    ]
    chosen_h = float(threshold_row["h"])

    rationale = {
        "window": {
            "criterion": "lowest score_std",
            "value": float(window_row["metric_value"]),
        },
        "xi_max": {
            "criterion": "lowest score_std",
            "value": float(frequency_row["score_std"]),
        },
        "h": {
            "criterion": f"alarm_rate closest to target_fpr={target_fpr:.4f}",
            "value": float(threshold_row["alarm_rate"]),
        },
    }

    recommendation = {
        "window": chosen_window,
        "xi_max": chosen_xi_max,
        "h": chosen_h,
        "rationale": rationale,
    }

    if verbose:
        print("CFAD parameter recommendation")
        print(f"  window : {chosen_window}")
        print(f"  xi_max : {chosen_xi_max:.3f}")
        print(f"  h      : {chosen_h:.3f}")

    return recommendation

Model Classes

Bases: CFModel

Normal-distribution characteristic-function model.

Source code in cfad/models/gaussian.py 
26
27
28
def __init__(self, mu: float = 0.0, sigma: float = 1.0):
    self.mu = mu
    self.sigma = sigma

is_analytic class-attribute instance-attribute 

is_analytic = True

mu instance-attribute 

mu = mu

sigma instance-attribute 

sigma = sigma

cf 

cf(xi: NDArray[float64]) -> NDArray[np.complex128]

Evaluate the characteristic function on the supplied frequency grid.

Source code in cfad/models/gaussian.py 
30
31
32
def cf(self, xi: NDArray[np.float64]) -> NDArray[np.complex128]:
    """Evaluate the characteristic function on the supplied frequency grid."""
    return np.exp(1j * self.mu * xi - 0.5 * self.sigma**2 * xi**2)

log_cf 

log_cf(xi: NDArray[float64]) -> NDArray[np.complex128]

Evaluate the log characteristic function.

Source code in cfad/models/gaussian.py 
34
35
36
def log_cf(self, xi: NDArray[np.float64]) -> NDArray[np.complex128]:
    """Evaluate the log characteristic function."""
    return 1j * self.mu * xi - 0.5 * self.sigma**2 * xi**2

fit 

fit(returns: NDArray[float64]) -> 'GaussianCF'

Fit mean and sample standard deviation by their usual estimators.

Source code in cfad/models/gaussian.py 
38
39
40
41
42
43
44
45
def fit(self, returns: NDArray[np.float64]) -> "GaussianCF":
    """Fit mean and sample standard deviation by their usual estimators."""
    values = np.asarray(returns, dtype=np.float64)
    if values.ndim != 1 or values.size < 2:
        raise ValueError("returns must be one-dimensional with at least 2 values")
    self.mu = float(np.mean(values))
    self.sigma = float(np.std(values, ddof=1))
    return self

__repr__ 

__repr__() -> str
Source code in cfad/models/gaussian.py 
47
48
def __repr__(self) -> str:
    return f"GaussianCF(mu={self.mu:.6f}, sigma={self.sigma:.6f})"

Bases: CFModel

NIG characteristic function model (non-analytic).

Source code in cfad/models/nig.py 
31
32
33
34
35
36
37
38
39
40
41
def __init__(
    self,
    alpha: float = 10.0,
    beta: float = 0.0,
    delta: float = 0.1,
    mu: float = 0.0,
):
    self.alpha = alpha
    self.beta = beta
    self.delta = delta
    self.mu = mu

is_analytic class-attribute instance-attribute 

is_analytic = False

alpha instance-attribute 

alpha = alpha

beta instance-attribute 

beta = beta

delta instance-attribute 

delta = delta

mu instance-attribute 

mu = mu

cf 

cf(xi: NDArray[float64]) -> NDArray[np.complex128]
Source code in cfad/models/nig.py 
43
44
45
46
47
def cf(self, xi: NDArray[np.float64]) -> NDArray[np.complex128]:
    a, b, d, m = self.alpha, self.beta, self.delta, self.mu
    gamma = np.sqrt(a**2 - b**2 + 0j)
    psi = np.sqrt(a**2 - (b + 1j * xi) ** 2 + 0j)
    return np.exp(1j * m * xi + d * (gamma - psi))

log_cf 

log_cf(xi: NDArray[float64]) -> NDArray[np.complex128]
Source code in cfad/models/nig.py 
49
50
51
52
53
def log_cf(self, xi: NDArray[np.float64]) -> NDArray[np.complex128]:
    a, b, d, m = self.alpha, self.beta, self.delta, self.mu
    gamma = np.sqrt(a**2 - b**2 + 0j)
    psi = np.sqrt(a**2 - (b + 1j * xi) ** 2 + 0j)
    return 1j * m * xi + d * (gamma - psi)

fit 

fit(returns: NDArray[float64]) -> 'NIGCF'

MLE via scipy minimize (L-BFGS-B).

Source code in cfad/models/nig.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def fit(self, returns: NDArray[np.float64]) -> "NIGCF":
    """MLE via scipy minimize (L-BFGS-B)."""
    from cfad.empirical_cf import ecf_at

    xi = np.linspace(-15, 15, 256)
    ecf = ecf_at(returns, xi)

    def neg_ecf_distance(params):
        a, b, d, m = params
        if a <= 0 or d <= 0 or abs(b) >= a:
            return 1e10
        self.alpha, self.beta, self.delta, self.mu = a, b, d, m
        phi = self.cf(xi)
        return float(np.sum(np.abs(ecf - phi) ** 2))

    x0 = [self.alpha, self.beta, self.delta, np.mean(returns)]
    res = minimize(neg_ecf_distance, x0, method="Nelder-Mead",
                   options={"maxiter": 5000, "xatol": 1e-6})
    self.alpha, self.beta, self.delta, self.mu = res.x
    return self

__repr__ 

__repr__() -> str
Source code in cfad/models/nig.py 
76
77
78
def __repr__(self) -> str:
    return (f"NIGCF(alpha={self.alpha:.4f}, beta={self.beta:.4f}, "
            f"delta={self.delta:.4f}, mu={self.mu:.6f})")

Bases: CFModel

CGMY characteristic-function model (non-analytic).

Source code in cfad/models/cgmy.py 
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def __init__(
    self,
    C: float = 1.0,
    G: float = 5.0,
    M: float = 10.0,
    Y: float = 0.5,
):
    if Y >= 2 or Y <= 0:
        raise ValueError(f"Y must satisfy 0 < Y < 2; got {Y}")
    if C <= 0 or G <= 0 or M <= 0:
        raise ValueError("C, G, M must be positive")
    self.C = C
    self.G = G
    self.M = M
    self.Y = Y

is_analytic class-attribute instance-attribute 

is_analytic = False

C instance-attribute 

C = C

G instance-attribute 

G = G

M instance-attribute 

M = M

Y instance-attribute 

Y = Y

log_cf 

log_cf(xi: NDArray[float64]) -> NDArray[np.complex128]

Evaluate the CGMY log characteristic function.

At Y=1 the standard expression has the removable singularity Gamma(-1) * 0. Taking the limit in Y gives

C * [(M-i*xi) log(M-i*xi) - M log(M) + (G+i*xi) log(G+i*xi) - G log(G)].

The characteristic-function normalization log(phi(0)) = 0 is imposed exactly after evaluation. This avoids cancellation error near Y=1 without perturbing the model at nonzero frequencies.

Source code in cfad/models/cgmy.py 
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def log_cf(self, xi: NDArray[np.float64]) -> NDArray[np.complex128]:
    """Evaluate the CGMY log characteristic function.

    At ``Y=1`` the standard expression has the removable singularity
    ``Gamma(-1) * 0``. Taking the limit in ``Y`` gives

    ``C * [(M-i*xi) log(M-i*xi) - M log(M)
         + (G+i*xi) log(G+i*xi) - G log(G)]``.

    The characteristic-function normalization ``log(phi(0)) = 0`` is
    imposed exactly after evaluation. This avoids cancellation error near
    ``Y=1`` without perturbing the model at nonzero frequencies.
    """
    C, G, M, Y = self.C, self.G, self.M, self.Y
    xi_arr = np.asarray(xi, dtype=np.float64)

    if np.isclose(Y, 1.0, atol=1e-12, rtol=0.0):
        right = M - 1j * xi_arr
        left = G + 1j * xi_arr
        limit = (
            right * np.log(right)
            - M * np.log(M)
            + left * np.log(left)
            - G * np.log(G)
        )
        result = np.asarray(C * limit, dtype=np.complex128)
    else:
        gam = gamma_func(-Y)
        term1 = (M - 1j * xi_arr) ** Y - M**Y
        term2 = (G + 1j * xi_arr) ** Y - G**Y
        result = np.asarray(C * gam * (term1 + term2), dtype=np.complex128)

    return np.where(xi_arr == 0.0, 0.0 + 0.0j, result).astype(np.complex128)

cf 

cf(xi: NDArray[float64]) -> NDArray[np.complex128]
Source code in cfad/models/cgmy.py 
91
92
def cf(self, xi: NDArray[np.float64]) -> NDArray[np.complex128]:
    return np.exp(self.log_cf(xi))

fit 

fit(returns: NDArray[float64]) -> 'CGMYCF'

ECF minimum-distance estimation via Nelder-Mead.

Source code in cfad/models/cgmy.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def fit(self, returns: NDArray[np.float64]) -> "CGMYCF":
    """ECF minimum-distance estimation via Nelder-Mead."""
    from cfad.empirical_cf import ecf_at

    xi = np.linspace(-10, 10, 128)
    ecf = ecf_at(returns, xi)

    def objective(params):
        C_, G_, M_, Y_ = params
        if C_ <= 0 or G_ <= 0 or M_ <= 0 or Y_ <= 0 or Y_ >= 2:
            return 1e10
        try:
            self.C, self.G, self.M, self.Y = C_, G_, M_, Y_
            phi = self.cf(xi)
            if not np.all(np.isfinite(phi)):
                return 1e10
            return float(np.sum(np.abs(ecf - phi) ** 2))
        except Exception:
            return 1e10

    x0 = [self.C, self.G, self.M, self.Y]
    res = minimize(
        objective,
        x0,
        method="Nelder-Mead",
        options={"maxiter": 8000, "xatol": 1e-5},
    )
    self.C, self.G, self.M, self.Y = res.x
    return self

__repr__ 

__repr__() -> str
Source code in cfad/models/cgmy.py 
124
125
126
127
128
def __repr__(self) -> str:
    return (
        f"CGMYCF(C={self.C:.4f}, G={self.G:.4f}, "
        f"M={self.M:.4f}, Y={self.Y:.4f})"
    )

Bases: CFModel

Lévy-stable characteristic function model (non-analytic for alpha<2).

Source code in cfad/models/levy_stable.py 
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
def __init__(
    self,
    alpha: float = 1.7,
    beta: float = 0.0,
    c: float = 0.01,
    mu: float = 0.0,
):
    if not (0 < alpha <= 2):
        raise ValueError(f"alpha must be in (0, 2]; got {alpha}")
    if not (-1 <= beta <= 1):
        raise ValueError(f"beta must be in [-1, 1]; got {beta}")
    if c <= 0:
        raise ValueError("c must be positive")
    self.alpha = alpha
    self.beta = beta
    self.c = c
    self.mu = mu

is_analytic class-attribute instance-attribute 

is_analytic = False

alpha instance-attribute 

alpha = alpha

beta instance-attribute 

beta = beta

c instance-attribute 

c = c

mu instance-attribute 

mu = mu

log_cf 

log_cf(xi: NDArray[float64]) -> NDArray[np.complex128]
Source code in cfad/models/levy_stable.py 
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def log_cf(self, xi: NDArray[np.float64]) -> NDArray[np.complex128]:
    a, b, c, m = self.alpha, self.beta, self.c, self.mu
    xi = np.asarray(xi, dtype=np.float64)
    cxi = c * xi
    abs_cxi = np.abs(cxi)

    if abs(a - 1.0) < 1e-8:
        # Special case alpha = 1
        return (
            -abs_cxi
            + 1j * b * (2 / np.pi) * cxi * np.log(abs_cxi + 1e-300)
            + 1j * m * xi
        )
    else:
        tan_factor = np.tan(np.pi * a / 2)
        skew = -1j * b * np.sign(xi) * tan_factor * (abs_cxi**a - abs_cxi)
        return -(abs_cxi**a) + skew + 1j * m * xi

cf 

cf(xi: NDArray[float64]) -> NDArray[np.complex128]
Source code in cfad/models/levy_stable.py 
70
71
def cf(self, xi: NDArray[np.float64]) -> NDArray[np.complex128]:
    return np.exp(self.log_cf(xi))

fit 

fit(returns: NDArray[float64]) -> 'LevyStableCF'

Maximum likelihood estimation via scipy.stats.levy_stable.fit.

Source code in cfad/models/levy_stable.py 
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def fit(self, returns: NDArray[np.float64]) -> "LevyStableCF":
    """Maximum likelihood estimation via scipy.stats.levy_stable.fit."""
    from scipy.stats import levy_stable

    returns_arr = np.asarray(returns, dtype=np.float64)
    alpha, beta, loc, scale = levy_stable.fit(returns_arr)
    if not (0 < alpha <= 2):
        raise ValueError(f"fitted alpha must be in (0, 2]; got {alpha}")
    if not (-1 <= beta <= 1):
        raise ValueError(f"fitted beta must be in [-1, 1]; got {beta}")
    if scale <= 0:
        raise ValueError(f"fitted scale must be positive; got {scale}")

    self.alpha = float(alpha)
    self.beta = float(beta)
    self.c = float(scale)
    self.mu = float(loc)
    return self

__repr__ 

__repr__() -> str
Source code in cfad/models/levy_stable.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
def __repr__(self) -> str:
    return (
        f"LevyStableCF(alpha={self.alpha:.4f}, beta={self.beta:.4f}, "
        f"c={self.c:.6f}, mu={self.mu:.6f})"
    )
    return (
        f"LevyStableCF(alpha={self.alpha:.4f}, beta={self.beta:.4f}, "
        f"c={self.c:.6f}, mu={self.mu:.6f})"
    )