// Licensed under the Apache License, Version 2.0
// © mihakralj
//@version=6
indicator("Normal Distribution CDF (NORMDIST)", "NORMDIST", overlay=false, precision=6)

//@function Computes Normal Distribution CDF for a normalized price series
//@param source Series to transform
//@param period Lookback period for z-score normalization
//@param mu Mean parameter (0.0 for standard normal after z-score)
//@param sigma Standard deviation parameter (1.0 for standard normal after z-score)
//@returns CDF value in [0,1]: Φ(z) = 0.5 × (1 + erf(z / √2))
//@optimized O(period) per bar for mean/variance scan; CDF itself is O(1)
normdist(series float source, simple int period, simple float mu, simple float sigma) =>
    if period <= 0
        runtime.error("Period must be greater than 0")
    if sigma <= 0.0
        runtime.error("Sigma must be greater than 0")

    // Compute rolling mean and standard deviation over lookback
    float sum = 0.0
    float sumSq = 0.0
    int count = 0
    for i = 0 to period - 1
        float v = source[i]
        if not na(v)
            sum += v
            sumSq += v * v
            count += 1

    float result = 0.5
    if count >= 2
        float mean = sum / count
        float variance = (sumSq / count) - (mean * mean)
        float stddev = variance > 0.0 ? math.sqrt(variance) : 0.0

        // Z-score: normalize source relative to its own rolling distribution
        float z = stddev > 0.0 ? (source - mean) / stddev : 0.0

        // Apply user-specified mu/sigma shift: z_final = (z - mu) / sigma
        float z_final = (z - mu) / sigma

        // Approximate erf via Abramowitz & Stegun (max error < 1.5e-7)
        // erf(x) = 1 - (a1*t + a2*t^2 + a3*t^3) * exp(-x^2)
        // where t = 1 / (1 + 0.47047 * |x|)
        float x = z_final / math.sqrt(2.0)
        float ax = math.abs(x)
        float t = 1.0 / (1.0 + 0.47047 * ax)
        float t2 = t * t
        float t3 = t2 * t

        float a1 = 0.3480242
        float a2 = -0.0958798
        float a3 = 0.7478556

        float erfApprox = 1.0 - (a1 * t + a2 * t2 + a3 * t3) * math.exp(-(ax * ax))
        float erf = x >= 0.0 ? erfApprox : -erfApprox

        // CDF: Φ(z) = 0.5 * (1 + erf(z / sqrt(2)))
        result := 0.5 * (1.0 + erf)

    result

// ---------- Main loop ----------

// Inputs
i_source = input.source(close, "Source")
i_period = input.int(50, "Lookback Period", minval=2, maxval=5000, tooltip="Rolling window for z-score normalization")
i_mu     = input.float(0.0, "Mu (μ)", step=0.1, tooltip="Mean shift parameter (0 = standard normal)")
i_sigma  = input.float(1.0, "Sigma (σ)", minval=0.01, step=0.1, tooltip="Scale parameter (1 = standard normal)")

// Calculation
float result = normdist(i_source, i_period, i_mu, i_sigma)

// Plot
plot(result, "NORMDIST", color=color.yellow, linewidth=2)
hline(0.5, "Midline", color=color.gray, linestyle=hline.style_dotted)
hline(0.975, "Upper 2σ", color=color.red, linestyle=hline.style_dashed)
hline(0.025, "Lower 2σ", color=color.green, linestyle=hline.style_dashed)
hline(0.841, "Upper 1σ", color=color.orange, linestyle=hline.style_dashed)
hline(0.159, "Lower 1σ", color=color.teal, linestyle=hline.style_dashed)