// Licensed under the Apache License, Version 2.0
// © mihakralj
//@version=6
indicator("Partial Autocorrelation Function (PACF)", "PACF", overlay=false)

//@function Calculates partial autocorrelation at a specified lag using Durbin-Levinson recursion
//@param src series float Input data series
//@param len simple int Lookback period for calculation
//@param lag simple int Lag order for partial autocorrelation (default 1)
//@returns float Partial autocorrelation coefficient between -1 and 1
//@optimized using Durbin-Levinson recursion to remove shorter-lag effects
pacf(series float src, simple int len, simple int lag) =>
    if len <= 0
        runtime.error("Period must be greater than 0")
    if lag < 1
        runtime.error("Lag must be at least 1")
    if len <= lag + 1
        runtime.error("Period must be greater than lag + 1")
    var int p = math.max(1, len)
    var array<float> buffer = array.new_float(p, na)
    var int head = 0, var int count = 0
    float oldest = array.get(buffer, head)
    if not na(oldest)
        count -= 1
    if not na(src)
        array.set(buffer, head, src)
        count += 1
    else
        array.set(buffer, head, na)
    head := (head + 1) % p
    if count <= lag
        na
    else
        // Calculate mean
        float sum = 0.0
        int validN = 0
        for i = 0 to p - 1
            float val = array.get(buffer, i)
            if not na(val)
                sum += val
                validN += 1
        if validN <= lag
            na
        else
            float mean = sum / validN

            // Calculate variance (population): Σ(x - mean)² / n
            float variance = 0.0
            for i = 0 to p - 1
                float val = array.get(buffer, i)
                if not na(val)
                    float diff = val - mean
                    variance += diff * diff
            variance /= validN

            if variance <= 0
                0.0
            else
                // Calculate ACF for lags 0 to target lag
                int startIdx = count < p ? 0 : head
                array<float> acfValues = array.new_float(lag + 1, 0.0)
                array.set(acfValues, 0, 1.0)

                for k = 1 to lag
                    float autocovariance = 0.0
                    for t = k to count - 1
                        int currentIdx = (startIdx + t) % p
                        int laggedIdx = (startIdx + t - k) % p
                        float xt = array.get(buffer, currentIdx)
                        float xtk = array.get(buffer, laggedIdx)
                        if not na(xt) and not na(xtk)
                            autocovariance += (xt - mean) * (xtk - mean)
                    autocovariance /= validN
                    array.set(acfValues, k, autocovariance / variance)

                // Durbin-Levinson recursion
                float pacfResult = na
                if lag == 1
                    // PACF at lag 1 equals ACF at lag 1
                    pacfResult := array.get(acfValues, 1)
                else
                    array<float> phi = array.new_float(lag + 1, 0.0)
                    array<float> phiPrev = array.new_float(lag + 1, 0.0)

                    // Initialize: φ_11 = r_1
                    array.set(phi, 1, array.get(acfValues, 1))

                    // Iterate for k = 2 to target lag
                    for k = 2 to lag
                        // Copy phi to phiPrev
                        for j = 0 to lag
                            array.set(phiPrev, j, array.get(phi, j))

                        // numerator = r_k - Σ(φ_{k-1,j} * r_{k-j}) for j=1..k-1
                        float numerator = array.get(acfValues, k)
                        for j = 1 to k - 1
                            numerator -= array.get(phiPrev, j) * array.get(acfValues, k - j)

                        // denominator = 1 - Σ(φ_{k-1,j} * r_j) for j=1..k-1
                        float denominator = 1.0
                        for j = 1 to k - 1
                            denominator -= array.get(phiPrev, j) * array.get(acfValues, j)

                        if math.abs(denominator) < 1e-15
                            pacfResult := 0.0
                            break

                        // φ_kk = numerator / denominator
                        array.set(phi, k, numerator / denominator)

                        // Update: φ_kj = φ_{k-1,j} - φ_kk * φ_{k-1,k-j}
                        for j = 1 to k - 1
                            array.set(phi, j, array.get(phiPrev, j) - array.get(phi, k) * array.get(phiPrev, k - j))

                    if na(pacfResult)
                        pacfResult := array.get(phi, lag)

                // Clamp to [-1, 1]
                math.max(-1.0, math.min(1.0, pacfResult))

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

// Inputs
i_source = input.source(close, "Source")
i_period = input.int(20, "Period", minval=3)
i_lag = input.int(1, "Lag", minval=1)

// Calculation
pacf_value = pacf(i_source, i_period, i_lag)

// Plot
plot(pacf_value, "PACF", color=color.yellow, linewidth=2)