645 lines
28 KiB
Python
645 lines
28 KiB
Python
"""
|
|
Wavelet band decomposition (CWT + synchrosqueezed variant), ported near-verbatim from
|
|
c:\\DataS\\InstrumentSimulator\\backend\\app\\wavelet.py (project "Macro Causal Lab").
|
|
|
|
No project-specific dependencies — pure numpy/ssqueezepy, safe to call from any service.
|
|
"""
|
|
import logging
|
|
import warnings
|
|
|
|
import numpy as np
|
|
from ssqueezepy import cwt, icwt
|
|
from ssqueezepy.utils.cwt_utils import center_frequency
|
|
|
|
logger = logging.getLogger(__name__)
|
|
|
|
MIN_SCALES_PER_BAND = 4
|
|
|
|
# Fixed period (days) lower-bounds for bands 0..5, independent of the analysis
|
|
# window's length. Without this, dividing [scales.min(), scales.max()] into
|
|
# num_levels *equal* log-spaced buckets means "band i" covers a totally
|
|
# different real-world period range depending on how much history is fed in
|
|
# (a longer window resolves longer scales, which shifts every intermediate
|
|
# boundary upward) - so the same band index is not comparable across window
|
|
# lengths. Anchoring to absolute day-periods keeps "band i" meaning the same
|
|
# oscillation whether the window is 3 months or 3 years; only the last band's
|
|
# upper bound stays open, since a longer window can genuinely resolve slower
|
|
# cycles a short one cannot.
|
|
CANONICAL_PERIOD_EDGES_DAYS = [0.3, 1.5, 4.0, 12.0, 35.0, 90.0]
|
|
|
|
|
|
def _period_to_log_scale_fn(wavelet, scales, n, num_calib_points=16):
|
|
"""Build an interpolation from log(period_days) to log(scale) by sampling
|
|
center_frequency at a handful of scales (not every scale - it's not
|
|
vectorized and can be slow over hundreds/thousands of scales).
|
|
"""
|
|
log_scales_full = np.log(scales)
|
|
calib_log_scales = np.linspace(log_scales_full.min(), log_scales_full.max(), min(num_calib_points, len(scales)))
|
|
calib_periods = np.array([1 / center_frequency(wavelet, scale=float(np.exp(s)), N=n) for s in calib_log_scales])
|
|
calib_log_periods = np.log(calib_periods)
|
|
|
|
def period_to_log_scale(period_days):
|
|
return float(np.interp(np.log(period_days), calib_log_periods, calib_log_scales))
|
|
|
|
return period_to_log_scale
|
|
|
|
|
|
def band_decompose(
|
|
values: list[float],
|
|
dates: list[str],
|
|
num_levels: int = 4,
|
|
wavelet: str = "gmw",
|
|
) -> dict:
|
|
x = np.asarray(values, dtype=float)
|
|
n = len(x)
|
|
if n < 32:
|
|
raise ValueError("Serie trop courte pour une analyse ondelette (32 points minimum).")
|
|
|
|
x_mean = float(x.mean())
|
|
xc = x - x_mean
|
|
|
|
Wx, scales = cwt(xc, wavelet=wavelet)
|
|
log_scales = np.log(scales)
|
|
|
|
period_to_log_scale = _period_to_log_scale_fn(wavelet, scales, n)
|
|
period_edges = list(CANONICAL_PERIOD_EDGES_DAYS[:num_levels])
|
|
while len(period_edges) < num_levels:
|
|
period_edges.append(period_edges[-1] * 2.5)
|
|
|
|
# scales[0] = smallest scale = highest frequency = shortest period (impulses/noise);
|
|
# scales[-1] = largest scale = lowest frequency = longest period (slow trend).
|
|
bands = []
|
|
for i in range(num_levels):
|
|
lo = period_to_log_scale(period_edges[i])
|
|
hi = period_to_log_scale(period_edges[i + 1]) if i + 1 < len(period_edges) else log_scales.max()
|
|
mask = (log_scales >= lo) & (log_scales <= hi if i == num_levels - 1 else log_scales < hi)
|
|
scales_band = scales[mask]
|
|
|
|
reconstruction_failed = False
|
|
if len(scales_band) >= MIN_SCALES_PER_BAND:
|
|
try:
|
|
recon = icwt(Wx[mask, :], wavelet=wavelet, scales=scales_band, x_len=n)
|
|
except Exception:
|
|
# Root-caused: icwt's internal scaletype inference recursively splits
|
|
# `scales` looking for a log/linear transition point (infer_scaletype ->
|
|
# logscale_transition_idx), and on a masked SUBSET of the original CWT
|
|
# scales that recursive sub-slice can come out too short, crashing with
|
|
# "attempt to get argmax of an empty sequence" instead of a clean error.
|
|
# Confirmed by reproducing it directly against ssqueezepy 0.6.6: retrying
|
|
# with a regenerated, exactly log-uniform array over the same range/count
|
|
# sidesteps the fragile inference entirely (trivially satisfies its "is
|
|
# this log-spaced" check) and reliably recovers a real reconstruction —
|
|
# verified on 5 independently reproduced failures, all recovered.
|
|
try:
|
|
clean_scales = np.geomspace(scales_band.min(), scales_band.max(), len(scales_band))
|
|
recon = icwt(Wx[mask, :], wavelet=wavelet, scales=clean_scales, x_len=n)
|
|
logger.warning(
|
|
"band_decompose: icwt failed for band %d with original scales (n=%d, wavelet=%s) — "
|
|
"recovered using a regenerated log-uniform scale array", i, len(scales_band), wavelet,
|
|
)
|
|
except Exception:
|
|
logger.warning(
|
|
"band_decompose: icwt failed for band %d (n=%d scales, wavelet=%s) even after "
|
|
"retry — falling back to a zero band", i, len(scales_band), wavelet, exc_info=True,
|
|
)
|
|
recon = np.zeros(n)
|
|
reconstruction_failed = True
|
|
scale_lo, scale_hi = float(scales_band.min()), float(scales_band.max())
|
|
else:
|
|
# Too few scales in this bucket for a stable reconstruction; report a
|
|
# flat zero band (its content ends up folded into the residual) so the
|
|
# band list always stays num_levels long.
|
|
logger.warning(
|
|
"band_decompose: only %d scale(s) in band %d (< %d minimum) — falling back to a zero band",
|
|
len(scales_band), i, MIN_SCALES_PER_BAND,
|
|
)
|
|
recon = np.zeros(n)
|
|
reconstruction_failed = True
|
|
scale_lo, scale_hi = float(np.exp(lo)), float(np.exp(hi))
|
|
|
|
freq_at_scale_hi = center_frequency(wavelet, scale=scale_hi, N=n)
|
|
freq_at_scale_lo = center_frequency(wavelet, scale=scale_lo, N=n)
|
|
period_low_days = round(1 / freq_at_scale_lo, 1) if freq_at_scale_lo else None
|
|
period_high_days = round(1 / freq_at_scale_hi, 1) if freq_at_scale_hi else None
|
|
label = (
|
|
f"{period_low_days}-{period_high_days}j"
|
|
if period_low_days is not None and period_high_days is not None
|
|
else f"bande {i + 1}"
|
|
)
|
|
bands.append(
|
|
{
|
|
"index": i,
|
|
"label": label,
|
|
"period_low_days": period_low_days,
|
|
"period_high_days": period_high_days,
|
|
"series": [round(float(v), 6) for v in recon],
|
|
"reconstruction_failed": reconstruction_failed,
|
|
}
|
|
)
|
|
|
|
reconstructed_total = x_mean + sum(np.array(band["series"]) for band in bands)
|
|
residual = x - reconstructed_total
|
|
|
|
return {
|
|
"dates": dates,
|
|
"original": [round(float(v), 6) for v in x],
|
|
"mean": round(x_mean, 6),
|
|
"bands": bands,
|
|
"residual": [round(float(v), 6) for v in residual],
|
|
"wavelet": wavelet,
|
|
}
|
|
|
|
|
|
def windowed_band_decompose(
|
|
values: list[float],
|
|
dates: list[str],
|
|
window_size: int,
|
|
num_levels: int = 4,
|
|
wavelet: str = "gmw",
|
|
) -> dict:
|
|
"""Decompose a long series by running band_decompose independently on
|
|
consecutive slices of `window_size` points, then concatenating the
|
|
results. Each slice is analyzed on its own terms (own local mean, own
|
|
CWT scale range), so precision near the reconstructed curve stays high
|
|
even over a long history; the trade-off is a visible discontinuity at
|
|
each slice boundary, which is acceptable for this diagnostic view.
|
|
"""
|
|
n = len(values)
|
|
window_size = max(32, min(window_size, n)) if window_size else n
|
|
if n <= window_size:
|
|
result = band_decompose(values, dates, num_levels, wavelet)
|
|
result["window_size"] = window_size
|
|
result["chunks"] = 1
|
|
return result
|
|
|
|
all_dates: list[str] = []
|
|
all_original: list[float] = []
|
|
all_residual: list[float] = []
|
|
band_series: list[list[float]] = [[] for _ in range(num_levels)]
|
|
band_meta: list[dict | None] = [None] * num_levels
|
|
band_failed: list[bool] = [False] * num_levels
|
|
|
|
chunk_count = 0
|
|
start = 0
|
|
while start < n:
|
|
end = min(start + window_size, n)
|
|
if 0 < n - end < 32:
|
|
# Don't leave a too-short trailing slice: fold the remainder in.
|
|
end = n
|
|
chunk = band_decompose(values[start:end], dates[start:end], num_levels, wavelet)
|
|
chunk_count += 1
|
|
all_dates.extend(chunk["dates"])
|
|
all_original.extend(chunk["original"])
|
|
all_residual.extend(chunk["residual"])
|
|
for band in chunk["bands"]:
|
|
band_series[band["index"]].extend(band["series"])
|
|
band_meta[band["index"]] = band
|
|
# OR across chunks: if any slice couldn't reconstruct this band, flag the
|
|
# whole concatenated series as unreliable rather than losing that signal
|
|
# when a later chunk happens to succeed and overwrites band_meta.
|
|
band_failed[band["index"]] = band_failed[band["index"]] or band.get("reconstruction_failed", False)
|
|
start = end
|
|
|
|
bands = []
|
|
for i in range(num_levels):
|
|
meta = band_meta[i] or {"label": f"bande {i + 1}", "period_low_days": None, "period_high_days": None}
|
|
bands.append(
|
|
{
|
|
"index": i,
|
|
"label": meta["label"],
|
|
"period_low_days": meta["period_low_days"],
|
|
"period_high_days": meta["period_high_days"],
|
|
"series": band_series[i],
|
|
"reconstruction_failed": band_failed[i],
|
|
}
|
|
)
|
|
|
|
return {
|
|
"dates": all_dates,
|
|
"original": all_original,
|
|
"mean": round(float(np.mean(all_original)), 6) if all_original else 0.0,
|
|
"bands": bands,
|
|
"residual": all_residual,
|
|
"wavelet": wavelet,
|
|
"window_size": window_size,
|
|
"chunks": chunk_count,
|
|
}
|
|
|
|
|
|
def rolling_causal_bands(
|
|
values: list[float],
|
|
dates: list[str],
|
|
start_idx: int,
|
|
lookback: int,
|
|
num_levels: int = 4,
|
|
wavelet: str = "gmw",
|
|
step: int = 1,
|
|
) -> dict:
|
|
"""Walk-forward band decomposition with no look-ahead: the band value
|
|
reported for day `t` is computed from a CWT run only on the trailing
|
|
`lookback` points ending at `t` (never anything after `t`). This is
|
|
deliberately much slower than `band_decompose` (one CWT per day instead
|
|
of one for the whole range) - that cost is the point: a single
|
|
whole-range CWT lets every point "see" the entire future through the
|
|
transform's global/symmetric support, which makes any backtest built on
|
|
it meaningless (the peaks/troughs it finds were computed with hindsight).
|
|
`step` > 1 re-runs the CWT every `step` days and holds the last computed
|
|
value in between, trading fidelity for speed on long ranges.
|
|
"""
|
|
n = len(values)
|
|
out_dates: list[str] = []
|
|
out_original: list[float] = []
|
|
band_series: list[list[float]] = [[] for _ in range(num_levels)]
|
|
band_labels = [f"bande {i + 1}" for i in range(num_levels)]
|
|
band_failed = [False] * num_levels
|
|
last_tip: list[float] | None = None
|
|
recomputations = 0
|
|
|
|
for t in range(start_idx, n):
|
|
window_start = t - lookback + 1
|
|
if window_start < 0:
|
|
continue # not enough trailing history yet for a full lookback window
|
|
if last_tip is None or (t - start_idx) % step == 0:
|
|
window_values = values[window_start:t + 1]
|
|
window_dates = dates[window_start:t + 1]
|
|
result = band_decompose(window_values, window_dates, num_levels, wavelet)
|
|
last_tip = [band["series"][-1] for band in result["bands"]]
|
|
band_labels = [band["label"] for band in result["bands"]]
|
|
for i, band in enumerate(result["bands"]):
|
|
band_failed[i] = band_failed[i] or band.get("reconstruction_failed", False)
|
|
recomputations += 1
|
|
out_dates.append(dates[t])
|
|
out_original.append(values[t])
|
|
for i in range(num_levels):
|
|
band_series[i].append(last_tip[i])
|
|
|
|
bands = [
|
|
{
|
|
"index": i,
|
|
"label": band_labels[i],
|
|
"period_low_days": None,
|
|
"period_high_days": None,
|
|
"series": band_series[i],
|
|
"reconstruction_failed": band_failed[i],
|
|
}
|
|
for i in range(num_levels)
|
|
]
|
|
|
|
return {
|
|
"dates": out_dates,
|
|
"original": out_original,
|
|
"bands": bands,
|
|
"wavelet": wavelet,
|
|
"lookback": lookback,
|
|
"step": step,
|
|
"recomputations": recomputations,
|
|
}
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Synchrosqueezed variant (opt-in, method="ssq"). Everything above this line
|
|
# is completely untouched by what follows: band_decompose, windowed_band_decompose
|
|
# and rolling_causal_bands are the exact same functions used when method="cwt"
|
|
# (the default), so any existing simulation re-run with its original config
|
|
# hits the same code path and produces byte-identical results as before.
|
|
#
|
|
# Rationale: plain cwt/icwt smears a single true frequency's energy across
|
|
# several neighboring scales (that's *why* band_decompose needs canonical period
|
|
# edges and a minimum-scales guard - the boundaries are fuzzy). Synchrosqueezing
|
|
# reassigns that smeared energy onto its estimated true instantaneous frequency
|
|
# before splitting into bands, which should reduce cross-band leakage. It also
|
|
# exposes two things plain cwt does not: per-band energy |Tx|^2 (how dominant
|
|
# that cycle is *right now*, independent of direction) and a dominant-frequency
|
|
# ridge (which cycle length is currently winning, tracked over time).
|
|
#
|
|
# Separate period edges from the cwt path's CANONICAL_PERIOD_EDGES_DAYS: with
|
|
# fs=1.0 (one sample = one day, matching our daily data), ssq_cwt has a hard
|
|
# Nyquist floor at exactly 2.0 days (max resolvable frequency = fs/2). The cwt
|
|
# path's edges start at 0.3 days - below that floor, so reusing them here would
|
|
# leave the fastest band permanently empty. This does mean "band 0" covers a
|
|
# different absolute range under ssq than under cwt; that's an inherent
|
|
# consequence of ssq's stricter frequency floor, not a rounding choice.
|
|
# ---------------------------------------------------------------------------
|
|
|
|
CANONICAL_PERIOD_EDGES_DAYS_SSQ = [2.0, 4.0, 10.0, 25.0, 60.0, 120.0]
|
|
|
|
|
|
def band_decompose_ssq(
|
|
values: list[float],
|
|
dates: list[str],
|
|
num_levels: int = 4,
|
|
wavelet: str = "gmw",
|
|
) -> dict:
|
|
from ssqueezepy import issq_cwt, ssq_cwt
|
|
from ssqueezepy.ridge_extraction import extract_ridges
|
|
|
|
x = np.asarray(values, dtype=float)
|
|
n = len(x)
|
|
if n < 32:
|
|
raise ValueError("Serie trop courte pour une analyse ondelette (32 points minimum).")
|
|
|
|
x_mean = float(x.mean())
|
|
xc = x - x_mean
|
|
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("ignore", RuntimeWarning)
|
|
# fs=1.0: one sample = one day, so ssq_freqs comes out directly in
|
|
# cycles/day and 1/ssq_freqs is a period in days with no extra
|
|
# calibration step (unlike band_decompose's scale->period conversion,
|
|
# which needs _period_to_log_scale_fn precisely because plain cwt's
|
|
# scale axis has no absolute unit by itself).
|
|
Tx, _Wx, ssq_freqs, _scales = ssq_cwt(xc, wavelet=wavelet, fs=1.0)
|
|
|
|
ssq_freqs = np.asarray(ssq_freqs)
|
|
periods = 1.0 / ssq_freqs
|
|
|
|
period_edges = list(CANONICAL_PERIOD_EDGES_DAYS_SSQ[:num_levels])
|
|
while len(period_edges) < num_levels:
|
|
period_edges.append(period_edges[-1] * 2.5)
|
|
|
|
bands = []
|
|
for i in range(num_levels):
|
|
lo_period = period_edges[i]
|
|
hi_period = period_edges[i + 1] if i + 1 < len(period_edges) else float(periods.max())
|
|
mask = (periods >= lo_period) & (periods <= hi_period if i == num_levels - 1 else periods < hi_period)
|
|
|
|
reconstruction_failed = False
|
|
if mask.any():
|
|
try:
|
|
Tx_band = np.zeros_like(Tx)
|
|
Tx_band[mask, :] = Tx[mask, :]
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("ignore", RuntimeWarning)
|
|
recon = np.real(issq_cwt(Tx_band, wavelet=wavelet))
|
|
except Exception:
|
|
# Defensive, mirroring band_decompose: never let one band's edge
|
|
# case fail the whole request — but log it and flag it (see
|
|
# band_decompose's identical fix for why this used to be silent).
|
|
logger.warning(
|
|
"band_decompose_ssq: issq_cwt failed for band %d (wavelet=%s) — "
|
|
"falling back to a zero band", i, wavelet, exc_info=True,
|
|
)
|
|
recon = np.zeros(n)
|
|
reconstruction_failed = True
|
|
energy = np.sum(np.abs(Tx[mask, :]) ** 2, axis=0)
|
|
else:
|
|
logger.warning("band_decompose_ssq: no frequencies fall in band %d — falling back to a zero band", i)
|
|
recon = np.zeros(n)
|
|
energy = np.zeros(n)
|
|
reconstruction_failed = True
|
|
|
|
period_low_days = round(float(lo_period), 1)
|
|
period_high_days = round(float(hi_period), 1)
|
|
bands.append(
|
|
{
|
|
"index": i,
|
|
"label": f"{period_low_days}-{period_high_days}j",
|
|
"period_low_days": period_low_days,
|
|
"period_high_days": period_high_days,
|
|
"series": [round(float(v), 6) for v in recon],
|
|
"energy": [round(float(v), 8) for v in energy],
|
|
"reconstruction_failed": reconstruction_failed,
|
|
}
|
|
)
|
|
|
|
reconstructed_total = x_mean + sum(np.array(band["series"]) for band in bands)
|
|
residual = x - reconstructed_total
|
|
|
|
try:
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("ignore", RuntimeWarning)
|
|
ridge_idxs = extract_ridges(Tx, ssq_freqs, n_ridges=1, bw=4)
|
|
ridge_periods = 1.0 / ssq_freqs[ridge_idxs.ravel()]
|
|
ridge_period_days = [round(float(v), 3) for v in ridge_periods]
|
|
except Exception:
|
|
ridge_period_days = [None] * n
|
|
|
|
return {
|
|
"dates": dates,
|
|
"original": [round(float(v), 6) for v in x],
|
|
"mean": round(x_mean, 6),
|
|
"bands": bands,
|
|
"residual": [round(float(v), 6) for v in residual],
|
|
"wavelet": wavelet,
|
|
"ridge_period_days": ridge_period_days,
|
|
}
|
|
|
|
|
|
def rolling_causal_bands_ssq(
|
|
values: list[float],
|
|
dates: list[str],
|
|
start_idx: int,
|
|
lookback: int,
|
|
num_levels: int = 4,
|
|
wavelet: str = "gmw",
|
|
step: int = 1,
|
|
) -> dict:
|
|
"""Synchrosqueezed counterpart to rolling_causal_bands - identical
|
|
walk-forward, no-look-ahead structure (same day-by-day trailing window,
|
|
same step logic), just calling band_decompose_ssq instead of
|
|
band_decompose. rolling_causal_bands itself is untouched.
|
|
"""
|
|
n = len(values)
|
|
out_dates: list[str] = []
|
|
out_original: list[float] = []
|
|
band_series: list[list[float]] = [[] for _ in range(num_levels)]
|
|
energy_series: list[list[float]] = [[] for _ in range(num_levels)]
|
|
ridge_series: list[float | None] = []
|
|
band_labels = [f"bande {i + 1}" for i in range(num_levels)]
|
|
band_failed = [False] * num_levels
|
|
last_tip: list[float] | None = None
|
|
last_energy_tip: list[float] | None = None
|
|
last_ridge_tip: float | None = None
|
|
recomputations = 0
|
|
|
|
for t in range(start_idx, n):
|
|
window_start = t - lookback + 1
|
|
if window_start < 0:
|
|
continue
|
|
if last_tip is None or (t - start_idx) % step == 0:
|
|
window_values = values[window_start:t + 1]
|
|
window_dates = dates[window_start:t + 1]
|
|
result = band_decompose_ssq(window_values, window_dates, num_levels, wavelet)
|
|
last_tip = [band["series"][-1] for band in result["bands"]]
|
|
last_energy_tip = [band["energy"][-1] for band in result["bands"]]
|
|
last_ridge_tip = result["ridge_period_days"][-1]
|
|
band_labels = [band["label"] for band in result["bands"]]
|
|
for i, band in enumerate(result["bands"]):
|
|
band_failed[i] = band_failed[i] or band.get("reconstruction_failed", False)
|
|
recomputations += 1
|
|
out_dates.append(dates[t])
|
|
out_original.append(values[t])
|
|
for i in range(num_levels):
|
|
band_series[i].append(last_tip[i])
|
|
energy_series[i].append(last_energy_tip[i])
|
|
ridge_series.append(last_ridge_tip)
|
|
|
|
bands = [
|
|
{
|
|
"index": i,
|
|
"label": band_labels[i],
|
|
"period_low_days": None,
|
|
"period_high_days": None,
|
|
"series": band_series[i],
|
|
"energy": energy_series[i],
|
|
"reconstruction_failed": band_failed[i],
|
|
}
|
|
for i in range(num_levels)
|
|
]
|
|
|
|
return {
|
|
"dates": out_dates,
|
|
"original": out_original,
|
|
"bands": bands,
|
|
"wavelet": wavelet,
|
|
"lookback": lookback,
|
|
"step": step,
|
|
"recomputations": recomputations,
|
|
"method": "ssq",
|
|
"ridge_period_days": ridge_series,
|
|
}
|
|
|
|
|
|
def _turning_points(series: list[float], smooth_days: int) -> list[tuple[int, int]]:
|
|
"""Indices (and new direction, +1/-1) where the sign of the `smooth_days`-lag slope of
|
|
`series` flips. Using the raw day-to-day slope would flag noise as a "reversal"; lagging
|
|
over a few days smooths that out."""
|
|
turns: list[tuple[int, int]] = []
|
|
prev_sign = 0
|
|
for t in range(smooth_days, len(series)):
|
|
slope = series[t] - series[t - smooth_days]
|
|
sign = 1 if slope > 0 else (-1 if slope < 0 else 0)
|
|
if sign != 0:
|
|
if prev_sign != 0 and sign != prev_sign:
|
|
turns.append((t, sign))
|
|
prev_sign = sign
|
|
return turns
|
|
|
|
|
|
def _average_cycle_days(turn_points: list[tuple[int, int]]) -> float | None:
|
|
"""Empirical average gap (in days) between consecutive SAME-direction turning points —
|
|
peak-to-peak or trough-to-trough, i.e. one full oscillation — measured from the actual
|
|
reconstructed series rather than assumed from the band's theoretical period_low/high
|
|
labels (which describe the wavelet scale bucket, not necessarily the real cycle length
|
|
this specific window's price action produced). None if there aren't at least 2
|
|
same-direction turns to measure a gap from."""
|
|
peaks = [t for t, d in turn_points if d == -1] # a downturn = just passed a peak
|
|
troughs = [t for t, d in turn_points if d == 1] # an upturn = just passed a trough
|
|
gaps = [seq[i + 1] - seq[i] for seq in (peaks, troughs) for i in range(len(seq) - 1)]
|
|
return sum(gaps) / len(gaps) if gaps else None
|
|
|
|
|
|
def wavelet_reliability(
|
|
values: list[float],
|
|
dates: list[str],
|
|
start_idx: int,
|
|
lookback: int,
|
|
num_levels: int = 4,
|
|
wavelet: str = "gmw",
|
|
method: str = "cwt",
|
|
step: int = 1,
|
|
smooth_days: int = 3,
|
|
tolerance_pct: float = 0.10,
|
|
min_confirm_horizon: int = 3,
|
|
) -> dict:
|
|
"""
|
|
For every turning point a CAUSAL (walk-forward, no look-ahead) decomposition flags —
|
|
i.e. what a live signal would have shown on that date, using only data available up to
|
|
then — checks whether redoing the decomposition later still shows a same-direction
|
|
reversal near the original date. The fraction confirmed is a per-band confidence score:
|
|
CWT/SSQ reconstructions are well known to be least reliable right at the edge of the
|
|
available data — exactly where a live "the band just turned" reading gets made.
|
|
|
|
A fixed confirm horizon (e.g. always +10 days) is meaningless across bands whose natural
|
|
oscillation period can range from ~1 day to several weeks: 10 days is many cycles for a
|
|
fast band (trivially "confirmed", inflating its score) but less than half a cycle for a
|
|
slow one (never has time to actually turn back, deflating its score) — the exact pattern
|
|
that shows up as fast bands scoring ~100% and the slow band scoring ~10% for no real
|
|
reliability reason. Instead, per band: measure its own historical average peak-to-peak
|
|
(or trough-to-trough) period from the causal series, and use `avg_cycle * (1 +
|
|
tolerance_pct)` as a MAXIMUM horizon — a same-direction reversal confirms the original one
|
|
if it shows up ANYWHERE in the forward window [original date, original date + that
|
|
horizon], not just near one specific point in it.
|
|
"""
|
|
rolling_decomposer = rolling_causal_bands_ssq if method == "ssq" else rolling_causal_bands
|
|
batch_decompose = band_decompose_ssq if method == "ssq" else band_decompose
|
|
|
|
causal = rolling_decomposer(
|
|
values, dates, start_idx=start_idx, lookback=lookback,
|
|
num_levels=num_levels, wavelet=wavelet, step=step,
|
|
)
|
|
causal_dates = causal["dates"]
|
|
date_to_global = {d: i for i, d in enumerate(dates)}
|
|
|
|
bands_out = []
|
|
for band in causal["bands"]:
|
|
turn_points = _turning_points(band["series"], smooth_days)
|
|
avg_cycle_days = _average_cycle_days(turn_points)
|
|
|
|
if avg_cycle_days is None:
|
|
bands_out.append({
|
|
"label": band["label"], "n_tested": 0, "n_confirmed": 0, "confidence_pct": None,
|
|
"avg_cycle_days": None, "confirm_horizon": None, "turns": [],
|
|
})
|
|
continue
|
|
|
|
confirm_horizon = max(min_confirm_horizon, round(avg_cycle_days * (1 + tolerance_pct)))
|
|
|
|
tested = []
|
|
for t, direction in turn_points:
|
|
d_date = causal_dates[t]
|
|
g_idx = date_to_global.get(d_date)
|
|
if g_idx is None or g_idx + confirm_horizon >= len(values):
|
|
continue # not enough real future data yet to judge this one
|
|
|
|
end_idx = g_idx + confirm_horizon
|
|
window_start = end_idx - lookback + 1
|
|
if window_start < 0:
|
|
continue
|
|
window_values = values[window_start:end_idx + 1]
|
|
window_dates = dates[window_start:end_idx + 1]
|
|
try:
|
|
hindsight = batch_decompose(window_values, window_dates, num_levels, wavelet)
|
|
except ValueError:
|
|
continue
|
|
h_series = hindsight["bands"][band["index"]]["series"]
|
|
|
|
# Forward-only window [original date, original date + horizon] — the reversal
|
|
# can be confirmed anywhere in it, not just near a specific point. window_values
|
|
# was built to end exactly at end_idx = g_idx + confirm_horizon, so h_series's
|
|
# last index already IS "original date + horizon"; no separate bound needed.
|
|
d_pos = g_idx - window_start
|
|
hindsight_turns = _turning_points(h_series, smooth_days)
|
|
matches = [idx for idx, new_dir in hindsight_turns if idx >= d_pos and new_dir == direction]
|
|
confirmed = bool(matches)
|
|
# Actual delay = how many days after the original date the confirming reversal
|
|
# showed up (earliest match) — compared against avg_cycle_days below to check
|
|
# whether "one cycle" is actually a well-calibrated horizon, or systematically
|
|
# too long/short relative to how quickly a real reversal actually confirms.
|
|
actual_delay = (min(matches) - d_pos) if confirmed else None
|
|
tested.append({"date": d_date, "confirmed": confirmed, "actual_delay_days": actual_delay})
|
|
|
|
n_tested = len(tested)
|
|
n_confirmed = sum(1 for t in tested if t["confirmed"])
|
|
delays = [t["actual_delay_days"] for t in tested if t["actual_delay_days"] is not None]
|
|
avg_actual_delay = round(sum(delays) / len(delays), 1) if delays else None
|
|
calibration_gap = round(sum(abs(avg_cycle_days - d) for d in delays) / len(delays), 1) if delays else None
|
|
bands_out.append({
|
|
"label": band["label"],
|
|
"n_tested": n_tested,
|
|
"n_confirmed": n_confirmed,
|
|
"confidence_pct": round(n_confirmed / n_tested * 100, 1) if n_tested else None,
|
|
"avg_cycle_days": round(avg_cycle_days, 1),
|
|
"confirm_horizon": confirm_horizon,
|
|
"avg_actual_delay_days": avg_actual_delay,
|
|
"calibration_gap_days": calibration_gap,
|
|
"turns": tested,
|
|
})
|
|
|
|
return {
|
|
"bands": bands_out,
|
|
"smooth_days": smooth_days,
|
|
"tolerance_pct": tolerance_pct,
|
|
"method": method,
|
|
}
|