Mood Test Change-Point Detector

//@version=6

// ============================================================================
// Mood Test Change-Point Detector
// ----------------------------------------------------------------------------
// Non-parametric, rank-based detector for volatility REGIME SHIFTS (not
// regime levels). Based on the sequential change-point framework in
// Ross, Tasoulis & Adams (2011) and applied to VIX/return series in
// Nystrup, Madsen & Lindstrom (2016) — see
// .claude/skills/quant-research-notes/references/change-point-detection-vix-nystrup.md
//
// What it does:
//   - Every bar, ranks the last N log-returns once.
//   - Because a point's rank inside the pooled window does not depend on
//     WHERE you split the window, the Mood scale statistic for every
//     candidate split point k=1..N-1 can be computed from a single
//     prefix-sum pass over those ranks (no re-ranking per split).
//   - D_max = the largest standardized Mood statistic across all candidate
//     splits. A large D_max means the recent segment of the window has a
//     significantly different variance (dispersion) than the older segment
//     — i.e., a volatility regime shift happened somewhere inside the window.
//   - This is deliberately about variance/scale shifts, not about price
//     direction. It answers "did the market's noise level just change?",
//     a question nothing else in this repo currently answers directly.
//
// Calibration note: true ARL (average run length) calibration requires
// offline Monte Carlo simulation per instrument/timeframe (as in the
// source papers). The threshold input here is a practical, empirically
// tuned stand-in — start around 3.0-3.5 and adjust to your instrument's
// alarm frequency rather than treating it as a calibrated p-value.
// ============================================================================

indicator(title = "Mood Test Change-Point Detector", shorttitle = "MoodCPD", overlay = false)

// ============================================================================
// INPUTS
// ============================================================================
srcInput           = input.source(close, "Price Source (used unless VIX enabled)", group = "Data")
useVixInput        = input.bool(false, "Use External VIX Series Instead of Price", group = "Data", tooltip = "Nystrup et al. (2016) found VIX-based change points outperform price-based ones for detecting equity volatility regime shifts. Only meaningful on instruments correlated with the VIX symbol below (SPX, ES, NQ, MNQ, etc.).")
vixSymbolInput     = input.symbol("CBOE:VIX", "VIX Symbol", group = "Data")

windowInput        = input.int(40, "Rolling Window Length (N)", minval = 20, maxval = 80, group = "Change-Point Test", tooltip = "Number of bars pooled for the rank test. Larger N = more stable ranks but slower to react and more compute per bar.")
minSplitFracInput  = input.float(0.20, "Minimum Split Fraction (edge guard)", minval = 0.05, maxval = 0.40, step = 0.05, group = "Change-Point Test", tooltip = "Candidate split points are restricted to [fraction*N, (1-fraction)*N] so both segments always have enough observations for a stable variance estimate.")
thresholdInput     = input.float(3.2, "Detection Threshold (standardized, empirical)", minval = 1.5, maxval = 6.0, step = 0.1, group = "Change-Point Test", tooltip = "D_max is a MAXIMUM over many correlated candidate splits, so a naive single-test 1.96/2.58 z-cutoff fires far too often. Start around 3.0-3.5 and raise if alarms feel too frequent for your instrument/timeframe.")

emaLenInput        = input.int(20, "EWMA Vol Re-estimation Length (bars)", minval = 5, maxval = 100, group = "Volatility Re-estimation", tooltip = "After an alarm, volatility is re-estimated from scratch with this EMA length, mirroring the source paper's re-estimation step following a detected change.")
baselineLenInput   = input.int(100, "Baseline (Long-Run) Vol EMA Length", minval = 20, maxval = 300, group = "Volatility Re-estimation")

showBgInput        = input.bool(true, "Highlight Alarm Bars", group = "Display")
showTableInput     = input.bool(true, "Show Info Table", group = "Display")

// ============================================================================
// RETURNS
// ============================================================================
float vixClose = useVixInput ? request.security(vixSymbolInput, timeframe.period, close, lookahead = barmerge.lookahead_off) : na
float src = useVixInput ? vixClose : srcInput
float r = math.log(src / src[1])

// ============================================================================
// ROLLING RANK-BASED MOOD STATISTIC
// ============================================================================
int N = windowInput
bool enoughBars = bar_index >= N

var float dMax = na
var int kStar = na
var float standardizedAtKStar = na

if enoughBars
    array<float> win = array.new<float>(N)
    for i = 0 to N - 1
        array.set(win, i, r[N - 1 - i])

    array<float> rnk = array.new<float>(N)
    for i = 0 to N - 1
        float vi = array.get(win, i)
        float lessCount = 0.0
        float eqCount = 0.0
        for j = 0 to N - 1
            float vj = array.get(win, j)
            if vj < vi
                lessCount += 1.0
            if vj == vi
                eqCount += 1.0
        array.set(rnk, i, lessCount + 1.0 + (eqCount - 1.0) / 2.0)

    array<float> dev2 = array.new<float>(N)
    float centerRank = (N + 1) / 2.0
    for i = 0 to N - 1
        float d = array.get(rnk, i) - centerRank
        array.set(dev2, i, d * d)

    int kMin = int(math.round(N * minSplitFracInput))
    int kMax = N - kMin
    float cumsum = 0.0
    float bestAbs = -1.0
    int bestK = na
    float bestStandardized = na
    for k = 1 to N - 1
        cumsum += array.get(dev2, k - 1)
        if k >= kMin and k <= kMax
            int m = k
            int n = N - k
            float expM = m * (N * N - 1) / 12.0
            float varM = m * n * (N + 1) * (N * N - 4) / 180.0
            float standardized = varM > 0 ? (cumsum - expM) / math.sqrt(varM) : 0.0
            if math.abs(standardized) > bestAbs
                bestAbs := math.abs(standardized)
                bestK := k
                bestStandardized := standardized

    dMax := bestAbs
    kStar := bestK
    standardizedAtKStar := bestStandardized

bool alarm = enoughBars and not na(dMax) and dMax > thresholdInput
bool volIncreasing = alarm and not na(standardizedAtKStar) and standardizedAtKStar < 0
bool volDecreasing = alarm and not na(standardizedAtKStar) and standardizedAtKStar > 0

// ============================================================================
// EWMA VOLATILITY RE-ESTIMATION (RESETS ON EACH ALARM)
// ============================================================================
var float emaSinceAlarm = na
float alpha = 2.0 / (emaLenInput + 1)
if na(emaSinceAlarm) or alarm
    emaSinceAlarm := r * r
else
    emaSinceAlarm := emaSinceAlarm + alpha * (r * r - emaSinceAlarm)

float baselineVol = ta.ema(r * r, baselineLenInput)

var int barsSinceAlarm = 0
if alarm
    barsSinceAlarm := 0
else
    barsSinceAlarm += 1

// ============================================================================
// PLOTS
// ============================================================================
plot(dMax, title = "D_max (Mood Statistic)", color = color.aqua, linewidth = 2)
hline(thresholdInput, title = "Threshold", color = color.red, linestyle = hline.style_dashed)

bgcolor(showBgInput and alarm ? color.new(color.yellow, 70) : na, title = "Alarm Highlight")

plotshape(volIncreasing, title = "Vol Increasing Alarm", style = shape.triangleup, location = location.top, color = color.orange, size = size.small)
plotshape(volDecreasing, title = "Vol Decreasing Alarm", style = shape.triangledown, location = location.top, color = color.blue, size = size.small)

plot(math.sqrt(emaSinceAlarm), title = "Vol Since Last Alarm (EWMA)", color = color.fuchsia, display = display.data_window)
plot(math.sqrt(baselineVol), title = "Baseline Vol (Long-Run EMA)", color = color.gray, display = display.data_window)

// ============================================================================
// INFO TABLE
// ============================================================================
if showTableInput and barstate.islast
    var table infoTable = table.new(position.bottom_right, 2, 5, bgcolor = color.new(color.rgb(15, 20, 27), 10), border_width = 1)
    table.cell(infoTable, 0, 0, "D_max", text_color = color.white)
    table.cell(infoTable, 1, 0, na(dMax) ? "N/A" : str.tostring(dMax, "#.##"), text_color = color.white)
    table.cell(infoTable, 0, 1, "Threshold", text_color = color.white)
    table.cell(infoTable, 1, 1, str.tostring(thresholdInput, "#.##"), text_color = color.white)
    table.cell(infoTable, 0, 2, "Direction", text_color = color.white)
    string dirWord = na(standardizedAtKStar) ? "N/A" : standardizedAtKStar > 0 ? "Vol Decreasing" : "Vol Increasing"
    table.cell(infoTable, 1, 2, dirWord, text_color = na(standardizedAtKStar) ? color.gray : standardizedAtKStar > 0 ? color.blue : color.orange)
    float volRatio = baselineVol > 0 ? math.sqrt(emaSinceAlarm) / math.sqrt(baselineVol) : na
    table.cell(infoTable, 0, 3, "Vol Ratio (since/baseline)", text_color = color.white)
    table.cell(infoTable, 1, 3, na(volRatio) ? "N/A" : str.tostring(volRatio, "#.##"), text_color = color.white)
    table.cell(infoTable, 0, 4, "Bars Since Alarm", text_color = color.white)
    table.cell(infoTable, 1, 4, str.tostring(barsSinceAlarm), text_color = color.white)