Ehlers Autocorrelation Periodogram

// The MIT License (MIT)
// © mihakralj
//@version=6
indicator("EACP: Ehlers Autocorrelation Periodogram","EACP",overlay=false)

//@function Autocorrelation periodogram dominant cycle estimator
//@param source Price input series
//@param minPeriod Minimum period to evaluate
//@param maxPeriod Maximum period to evaluate
//@param avgLength Averaging length for Pearson correlation (0 uses lag length)
//@param enhance Apply cubic emphasis to highlight dominant peaks
//@returns Smoothed dominant cycle estimate
//@optimized Removed buffer complexity, uses native PineScript historical operator for O(n) correlation
//@validation wolfram:"Wiener-Khinchin theorem","Pearson correlation coefficient" external:"TradingView TASC 2025.02 Autocorrelation","ImmortalFreedom Ehlers ACP","QuantStrat autocorrPeriodogram"
eacp(series float source,simple int minPeriod,simple int maxPeriod,simple int avgLength,simple bool enhance)=>
    if minPeriod<3
        runtime.error("Min period must be at least 3")
    if maxPeriod<=minPeriod
        runtime.error("Max period must be greater than min period")
    if avgLength<0
        runtime.error("Average length must be non-negative")
    int size=maxPeriod+1
    var array<float> corr=array.new_float(0)
    var array<float> power=array.new_float(0)
    var array<float> smooth=array.new_float(0)
    var int storedSize=0
    var int storedMin=0
    var int storedMax=0
    var bool configured=false
    var float hp=0.0
    var float filt=0.0
    var float dom=0.0
    var float domPower=0.0
    var float maxPwr=0.0
    var float e=1.0
    var bool warmup=true
    if not configured or storedSize!=size or storedMin!=minPeriod or storedMax!=maxPeriod
        corr:=array.new_float(size,0.0)
        power:=array.new_float(size,0.0)
        smooth:=array.new_float(size,0.0)
        storedSize:=size
        storedMin:=minPeriod
        storedMax:=maxPeriod
        configured:=true
        hp:=0.0
        filt:=0.0
        dom:=(minPeriod+maxPeriod)*0.5
        domPower:=0.0
        maxPwr:=0.0
        e:=1.0
        warmup:=true
    float price=nz(source)
    float alphaHP=(math.cos(math.sqrt(2.0)*math.pi/float(maxPeriod))+math.sin(math.sqrt(2.0)*math.pi/float(maxPeriod))-1.0)/math.cos(math.sqrt(2.0)*math.pi/float(maxPeriod))
    hp:=math.pow(1.0-alphaHP/2.0,2.0)*(price-2.0*nz(price[1])+nz(price[2]))+2.0*(1.0-alphaHP)*nz(hp[1])-math.pow(1.0-alphaHP,2.0)*nz(hp[2])
    float a1=math.exp(-math.sqrt(2.0)*math.pi/float(minPeriod))
    float b1=2.0*a1*math.cos(math.sqrt(2.0)*math.pi/float(minPeriod))
    float c2=b1
    float c3=-(a1*a1)
    float c1=1.0-c2-c3
    filt:=c1*(hp+nz(hp[1]))*0.5+c2*nz(filt[1])+c3*nz(filt[2])
    for lag=0 to maxPeriod
        if lag<2
            array.set(corr,lag,0.0)
        else
            int window=avgLength==0?lag:avgLength
            if window<2
                window:=2
            float sx=0.0
            float sy=0.0
            float sxx=0.0
            float syy=0.0
            float sxy=0.0
            int valid=0
            for k=0 to window-1
                float x=nz(filt[k])
                float y=nz(filt[lag+k])
                sx+=x
                sy+=y
                sxx+=x*x
                syy+=y*y
                sxy+=x*y
                valid+=1
            float corrVal=0.0
            if valid>1
                float denomX=float(valid)*sxx-sx*sx
                float denomY=float(valid)*syy-sy*sy
                float denom=denomX*denomY
                corrVal:=denom>0.0?(float(valid)*sxy-sx*sy)/math.sqrt(denom):0.0
            array.set(corr,lag,corrVal)
    for period=minPeriod to maxPeriod
        float cosAcc=0.0
        float sinAcc=0.0
        for n=2 to maxPeriod
            float corrVal=array.get(corr,n)
            float angle=2.0*math.pi*float(n)/float(period)
            cosAcc+=corrVal*math.cos(angle)
            sinAcc+=corrVal*math.sin(angle)
        float sq=cosAcc*cosAcc+sinAcc*sinAcc
        array.set(smooth,period,0.2*sq*sq+0.8*array.get(smooth,period))
    float localMaxPwr=0.0
    for period=minPeriod to maxPeriod
        float smoothVal=array.get(smooth,period)
        if smoothVal>localMaxPwr
            localMaxPwr:=smoothVal
    float diff=float(maxPeriod-minPeriod)
    float K=diff>0?math.pow(10.0,-0.15/diff):1.0
    if localMaxPwr>maxPwr
        maxPwr:=localMaxPwr
    else
        maxPwr:=K*maxPwr
    float weighted=0.0
    float sumWeight=0.0
    float peakPwr=0.0
    for period=minPeriod to maxPeriod
        float smoothVal=array.get(smooth,period)
        float pwr=maxPwr>0.0?smoothVal/maxPwr:0.0
        if enhance
            pwr:=math.pow(pwr,3.0)
        array.set(power,period,pwr)
        if pwr>peakPwr
            peakPwr:=pwr
        if pwr>=0.5
            weighted+=float(period)*pwr
            sumWeight+=pwr
    float base=sumWeight>=0.25?weighted/sumWeight:dom
    float alpha=0.2
    float beta=1.0-alpha
    dom:=alpha*(base-dom)+dom
    if warmup
        e*=beta
        float c=1.0/(1.0-e)
        dom:=c*dom
        warmup:=e>1e-10
    int domIdx=math.min(math.max(int(math.round(dom)),minPeriod),maxPeriod)
    domPower:=array.get(power,domIdx)
    [dom,domPower]

// ---------- Main loop ----------
i_source=input.source(close,"Source")
i_minPeriod=input.int(8,"Min Period",minval=3,maxval=500)
i_maxPeriod=input.int(48,"Max Period",minval=4,maxval=500)
i_avgLength=input.int(3,"Autocorrelation Length",minval=0,maxval=500)
i_enhance=input.bool(true,"Enhance Resolution")
[dominantCycle,normalizedPower]=eacp(i_source,i_minPeriod,i_maxPeriod,i_avgLength,i_enhance)
plot(dominantCycle,"Dominant Cycle",color=color.yellow,linewidth=2)
plot(normalizedPower,"Normalized Power",color=color.orange,linewidth=2)