Adaptivity: Measures of Dominant Cycles and Price Trend [Loxx]
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © loxx

//@version=5

indicator('Adaptivity: Measures of Dominant Cycles and Price Trend [Loxx]', shorttitle = "AMDCPT [Loxx]", overlay = false, max_bars_back = 2000)

greencolor = #2DD204
redcolor = #D2042D 
bluecolor = #042dd2
purplecolor = #a904d2


auto_src = input.source(close, title='Auto Source', group = "Dominant Cycle: Autocorrelation Periodogram")
auto_min = input.int(8, minval = 1, title='Auto Min Length', group = "Dominant Cycle: Autocorrelation Periodogram")
auto_max = input.int(48, minval = 1, title='Auto Max Length', group = "Dominant Cycle: Autocorrelation Periodogram")
auto_avg = input.int(3, minval = 1, title='Auto Average Length', group = "Dominant Cycle: Autocorrelation Periodogram")

hilbert_src = input.source(hl2, title='Hilbert Source', group = "Dominant Cycle: Hilbert Period")
hilbert_len = input.int(7, title='Hilbert Length', minval=1, group = "Dominant Cycle: Hilbert Period")
hilbert_alpha = input.float(0.07, title = "Hilbert Alpha", step = 0.01, group = "Dominant Cycle: Hilbert Period")

instant_src = input.source(hl2, title='Instant Source', group = "Dominant Cycle: Instantaneous")
instant_min_len = input.int(8, title='Instant Min Length', minval=1, group = "Dominant Cycle: Instantaneous")
instant_max_len = input.int(48, title='Instant Max Length', minval=1, group = "Dominant Cycle: Instantaneous")

bp_period = input.int(20, "Band-pass Period", minval = 1, group = "Dominant Cycle: Band-pass")
Bandwidth = input.float(0.70, "Band-pass Width", step = 0.1, group = "Dominant Cycle: Band-pass")

LPPeriod = input(20, title='LP Period', group = "Dominant Cycle: Hilbert Dual Differentiator")
min_len_input = input.int(8, title='Min Length', group = "Dominant Cycle: Hilbert Dual Differentiator")
max_len_input = input.int(48, title='Max Length', group = "Dominant Cycle: Hilbert Dual Differentiator")

cfb_src = input.source(hlcc4, "Source", group = "Price Trend Adaptive: Composite Fractal Behavior Settings")
cfb_len = input.int(10, "CFB Depth", maxval = 10, group = "Price Trend Adaptive: Composite Fractal Behavior Settings")
cfb_smth = input.int(8, "CFB Smooth Length", group = "Price Trend Adaptive: Composite Fractal Behavior Settings")
cfb_jper = input.int(5, "Juirk Fitler Smoothing Period", group = "Price Trend Adaptive: Composite Fractal Behavior Settings")
cfb_phs = input.float(0.0, "Jurik Filter Phase", group = "Price Trend Adaptive: Composite Fractal Behavior Settings")


_a_jurik_filt(src, len, phase) =>
    //static variales
    volty = 0.0, avolty = 0.0, vsum = 0.0, bsmax = src, bsmin = src
    len1 = math.max(math.log(math.sqrt(0.5 * (len-1))) / math.log(2.0) + 2.0, 0)
    len2 = math.sqrt(0.5 * (len - 1)) * len1
    pow1 = math.max(len1 - 2.0, 0.5)
    beta = 0.45 * (len - 1) / (0.45 * (len - 1) + 2)
    div = 1.0 / (10.0 + 10.0 * (math.min(math.max(len-10, 0), 100)) / 100)
    phaseRatio = phase < -100 ? 0.5 : phase > 100 ? 2.5 : 1.5 + phase * 0.01
    bet = len2 / (len2 + 1)

    //Price volatility
    del1 = src - nz(bsmax[1])
    del2 = src - nz(bsmin[1])
    volty := math.abs(del1) > math.abs(del2) ? math.abs(del1) : math.abs(del2)

    //Relative price volatility factor
    vsum := nz(vsum[1]) + div * (volty - nz(volty[10]))
    avolty := nz(avolty[1]) + (2.0 / (math.max(4.0 * len, 30) + 1.0)) * (vsum - nz(avolty[1]))
    dVolty = avolty > 0 ? volty / avolty : 0
    dVolty := math.max(1, math.min(math.pow(len1, 1.0/pow1), dVolty))
    
    //Jurik volatility bands
    pow2 = math.pow(dVolty, pow1)
    Kv = math.pow(bet, math.sqrt(pow2))
    bsmax := del1 > 0 ? src : src - Kv * del1
    bsmin := del2 < 0 ? src : src - Kv * del2
    
    //Jurik Dynamic Factor
    alpha = math.pow(beta, pow2)

    //1st stage - prelimimary smoothing by adaptive EMA
    jma = 0.0, ma1 = 0.0, det0 = 0.0, e2 = 0.0
    ma1 := (1 - alpha) * src + alpha * nz(ma1[1])
    
    //2nd stage - one more prelimimary smoothing by Kalman filter
    det0 := (src - ma1) * (1 - beta) + beta * nz(det0[1])
    ma2 = ma1 + phaseRatio * det0

    //3rd stage - final smoothing by unique Jurik adaptive filter
    e2 := (ma2 - nz(jma[1])) * math.pow(1 - alpha, 2) + math.pow(alpha, 2) * nz(e2[1])
    jma := e2 + nz(jma[1])
    jma

_jcfbaux(src, depth)=>
    jrc04 = 0.0, jrc05 = 0.0, jrc06 = 0.0, jrc13 = 0.0, jrc03 = 0.0, jrc08 = 0.0
    if (bar_index >= bar_index - depth * 2) 
        for k = depth - 1 to 0
            jrc04 += math.abs(nz(src[k]) - nz(src[k+1]))
            jrc05 += (depth + k) * math.abs(nz(src[k]) - nz(src[k+1]))
            jrc06 += nz(src[k+1])
    if(bar_index < bar_index - depth * 2)
        jrc03 := math.abs(src - nz(src[1]))
        jrc13 := math.abs(nz(src[depth]) - nz(src[depth+1]))
        jrc04 := jrc04 - jrc13 + jrc03
        jrc05 := jrc05 - jrc04 + jrc03 * depth
        jrc06 := jrc06 - nz(src[depth+1]) + nz(src[1])
    jrc08 := math.abs(depth * src - jrc06)
    jcfbaux = jrc05 == 0.0 ? 0.0 : jrc08/jrc05
    jcfbaux

_jcfb(src, len, smth) =>
    a1 = array.new<float>(0, 0)
    c1 = array.new<float>(0, 0)
    d1 = array.new<float>(0, 0)
    cumsum = 2.0, result = 0.0

    //crete an array of 2^index, pushed onto itself
    //max values with injest of depth = 10: [1, 1, 2, 2, 4, 4, 8, 8, 16, 16, 32, 32, 64, 64, 128, 128, 256, 256, 512, 512]
    for i = 0 to len -1
        array.push(a1, math.pow(2, i))
        array.push(a1, math.pow(2, i))
        
    //cumulative sum of 2^index array
    //max values with injest of depth = 10: [2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536]
    for i = 0 to array.size(a1) - 1 
        array.push(c1, cumsum)
        cumsum += array.get(a1, i)
        
    //add values from jcfbaux calculation using the cumumulative sum arrary above
    for i = 0 to array.size(c1) - 1 
        array.push(d1, _jcfbaux(src, array.get(c1, i)))
    
    baux = array.copy(d1)
    cnt = array.size(baux)
    arrE = array.new<float>(cnt, 0)
    arrX = array.new<float>(cnt, 0)
    
    if bar_index <= smth
        for j = 0 to bar_index - 1
            for k = 0 to cnt - 1
                eget = array.get(arrE, k)
                array.set(arrE, k, eget + nz(array.get(baux, k)[bar_index - j]))
        for k = 0 to cnt - 1
            eget = nz(array.get(arrE, k))
            array.set(arrE, k, eget / bar_index)
    else
        for k = 0 to cnt - 1
            eget = nz(array.get(arrE, k))
            array.set(arrE, k, eget + (nz(array.get(baux, k)) - nz(array.get(baux, k)[smth])) / smth)
    if bar_index > 5
        a = 1.0, b = 1.0
        for k = 0 to cnt - 1
            if (k % 2 == 0)
                array.set(arrX, k, b * nz(array.get(arrE, k)))
                b := b * (1 - nz(array.get(arrX, k)))
            else
                array.set(arrX, k, a * nz(array.get(arrE, k)))
                a := a * (1 - nz(array.get(arrX, k)))
        sq = 0.0, sqw = 0.0 
        for i = 0 to cnt - 1
            sqw += nz(array.get(arrX, i)) * nz(array.get(arrX, i)) * nz(array.get(c1, i))
        for i = 0 to array.size(arrX) - 1
            sq += math.pow(nz(array.get(arrX, i)), 2)
        result := sq == 0.0 ? 0.0 : sqw / sq
        result




_hilber_dual(len_in, min_len, max_len) =>
    alpha1 = 0.00
    HP = 0.00
    a1 = 0.00
    b1 = 0.00
    c1 = 0.00
    c2 = 0.00
    c3 = 0.00
    Filt = 0.00
    QFilt = 0.00
    Real = 0.00
    Quad = 0.00
    Imag = 0.00
    IPeak = 0.00
    QPeak = 0.00
    IDot = 0.00
    QDot = 0.00
    Period = 0.00
    DomCycle = 0.00
    
    alpha1 := (math.cos(0.707 * 2 * math.pi / max_len) + math.sin(0.707 * 2 * math.pi / max_len) - 1) / math.cos(0.707 * 2 * math.pi / max_len)
    HP := (1 - alpha1 / 2) * (1 - alpha1 / 2) * (close - 2 * close[1] + close[2]) + 2 * (1 - alpha1) * nz(HP[1]) - (1 - alpha1) * (1 - alpha1) * nz(HP[2])
    
    a1 := math.exp(-math.sqrt(2) * math.pi / LPPeriod)
    b1 := 2 * a1 * math.cos(math.sqrt(2) * math.pi / LPPeriod)
    c2 := b1
    c3 := -a1 * a1
    c1 := 1 - c2 - c3
    
    Filt := c1 * (HP + nz(HP[1])) / 2 + c2 * nz(Filt[1]) + c3 * nz(Filt[2])
    
    IPeak := 0.991 * nz(IPeak[1])
    IPeak := math.abs(Filt) > IPeak ? math.abs(Filt) : IPeak
    
    Real := Filt / IPeak
    Quad := Real - nz(Real[1])
    QPeak := 0.991 * nz(QPeak[1])
    QPeak := math.abs(Quad) > QPeak ? math.abs(Quad) : QPeak
    Imag := Quad / QPeak
    IDot := Real - nz(Real[1])
    QDot := Imag - nz(Imag[1])
    
    Period := Real * QDot - Imag * IDot != 0 ? 2 * math.pi * (Real * Real + Imag * Imag) / (-Real * QDot + Imag * IDot) : Period
    
    Period := math.min(math.max(Period, min_len), max_len)

    DomCycle := c1 * (Period + nz(Period[1])) / 2 + c2 * nz(DomCycle[1]) + c3 * nz(DomCycle[2])
    DomCycle

_corrrelation(x, y, len) => 
    lenMinusOne = len - 1
    meanx = 0.0, meany = 0.0
    for i=0.0 to lenMinusOne
        meanx := meanx + nz(x[i])
        meany := meany + nz(y[i])
    meanx := meanx / len
    meany := meany / len
    sumxy=0.0, sumx=0.0, sumy=0.0
    for i=0  to lenMinusOne
        sumxy := sumxy +    (nz(x[i]) - meanx) * (nz(y[i]) - meany)
        sumx  := sumx  + math.pow(nz(x[i]) - meanx, 2)
        sumy  := sumy  + math.pow(nz(y[i]) - meany, 2)
    sumxy / math.sqrt(sumy * sumx)

_cycleit(src, len, alpha) =>
    instPeriod = 0.0
    smooth = (src + 2 * nz(src[1]) + 2 * nz(src[2]) + nz(src[3])) / 6
    cycle = 0.0
    cycle := bar_index < 7 ? (src - 2 * nz(src[1]) + nz(src[2])) / 4 : 
         (1 - .5*alpha)*(1 - .5*alpha)*(smooth - 2*smooth[1] + smooth[2]) + 2*(1 - alpha)*cycle[1] - (1 - alpha)*(1 - alpha)*cycle[2]
    q1 = (0.0962 * cycle + 0.5769 * nz(cycle[2]) - 0.5769 * nz(cycle[4]) - 0.0962 * nz(cycle[6])) * (0.5 + 0.08 * nz(instPeriod[1]))
    i1 = nz(cycle[3])
    deltaPhase = q1 != 0 and nz(q1[1]) != 0 ? (i1 / q1 - nz(i1[1]) / nz(q1[1])) / (1 + i1 * nz(i1[1]) / (q1 * nz(q1[1]))) : 0
    deltaPhase := math.min(math.max(deltaPhase, 0.1), 1.1)
    medianDelta = ta.percentile_nearest_rank(deltaPhase, len, 50)
    dc = medianDelta != 0 ? math.pi*2 / medianDelta + 0.5 : 15
    instPeriod := 0.33 * dc + 0.67 * nz(instPeriod[1])
    period = 0.0
    period := 0.15 * instPeriod + 0.85 * nz(period[1])
    number = math.round(math.floor(ta.wma(period, 4)))
    number
 
_f_hp(_src, max_len) =>
    var c = 360 * math.pi / 180
    _alpha = (1 - math.sin(c / max_len)) / math.cos(c / max_len)
    _hp = 0.0
    _hp := 0.5 * (1 + _alpha) * (_src - nz(_src[1])) + _alpha * nz(_hp[1])
    _hp

_f_ess(_src, _len) =>
    var s = math.sqrt(2)
    _a = math.exp(-s * math.pi / _len)
    _b = 2 * _a * math.cos(s *  math.pi  / _len)
    _c2 = _b
    _c3 = -_a * _a
    _c1 = 1 - _c2 - _c3
    _out = 0.0
    _out := _c1 * (_src + nz(_src[1])) / 2 + _c2 * nz(_out[1], nz(_src[1], _src)) + _c3 * nz(_out[2], nz(_src[2], nz(_src[1], _src)))
    _out

_auto_dom(src, min_len, max_len, ave_len) =>
    var c = 2 * math.pi
    var s = math.sqrt(2)
    avglen = ave_len

    filt = _f_ess(_f_hp(src, max_len), min_len)
    arr_size = max_len * 2
    
    var corr = array.new_float(arr_size, initial_value=0)
    var cospart = array.new_float(arr_size, initial_value=0)
    var sinpart = array.new_float(arr_size, initial_value=0)
    var sqsum = array.new_float(arr_size, initial_value=0)
    var r1 = array.new_float(arr_size, initial_value=0)
    var r2 = array.new_float(arr_size, initial_value=0)
    var pwr = array.new_float(arr_size, initial_value=0)
    
    for lag = 0 to max_len by 1
        m = avglen == 0 ? lag : avglen
        Sx = 0.0
        Sy = 0.0
        Sxx = 0.0
        Syy = 0.0
        Sxy = 0.0
        for i = 0 to m - 1 by 1
            x = nz(filt[i])
            y = nz(filt[lag + i])
            Sx += x
            Sy += y
            Sxx += x * x
            Sxy += x * y
            Syy += y * y
            Syy
        if (m * Sxx - Sx * Sx) * (m * Syy - Sy * Sy) > 0
            array.set(corr, lag, (m * Sxy - Sx * Sy) / math.sqrt((m * Sxx - Sx * Sx) * (m * Syy - Sy * Sy)))
    
    for period = min_len to max_len by 1
        array.set(cospart, period, 0)
        array.set(sinpart, period, 0)
        for n = ave_len to max_len by 1
            array.set(cospart, period, nz(array.get(cospart, period)) + nz(array.get(corr, n)) * math.cos(c * n / period))
            array.set(sinpart, period, nz(array.get(sinpart, period)) + nz(array.get(corr, n)) * math.sin(c * n / period))
        array.set(sqsum, period, math.pow(nz(array.get(cospart, period)), 2) + math.pow(nz(array.get(sinpart, period)), 2))
    
    for period = min_len to max_len by 1
        array.set(r2, period, nz(array.get(r1, period)))
        array.set(r1, period, 0.2 * math.pow(nz(array.get(sqsum, period)), 2) + 0.8 * nz(array.get(r2, period)))

    maxpwr = 0.0
    
    for period = min_len to max_len by 1
        if nz(array.get(r1, period)) > maxpwr
            maxpwr := nz(array.get(r1, period))

    for period = ave_len to max_len by 1
        array.set(pwr, period, nz(array.get(r1, period)) / maxpwr)
    
    dominantcycle = 0.0
    peakpwr = 0.0
    
    for period = min_len to max_len by 1
        if nz(array.get(pwr, period)) > peakpwr 
            peakpwr := nz(array.get(pwr, period))

    spx = 0.0
    sp = 0.0

    for period = min_len to max_len by 1
        if peakpwr >= 0.25 and nz(array.get(pwr, period)) >= 0.25
            spx += period * nz(array.get(pwr, period))
            sp += nz(array.get(pwr, period))

    if sp != 0 
        dominantcycle := spx / sp
        
    if sp < 0.25 
        dominantcycle := dominantcycle[1]
        
    if dominantcycle < 1 
        dominantcycle := 1

    dom_in = math.min(math.max(dominantcycle, min_len), max_len)
    dom_in

    
_instant_cycle(src, llen, hlen) =>
    v1 = src - nz(src[llen])
    ip = 0.0
    ip := 1.25 * (nz(v1[4]) - 0.635 * nz(v1[2])) + 0.635 * nz(ip[3])
    qu = 0.0
    qu := nz(v1[2]) - 0.338 * v1 + 0.338 * nz(qu[2])
    phase = math.abs(ip + nz(ip[1])) > 0 ? math.atan(math.abs((qu + nz(qu[1])) / (ip + nz(ip[1])))) * 180 / math.pi : 0
    phase := ip < 0 and qu > 0 ? 180 - phase : phase
    phase := ip < 0 and qu < 0 ? 180 + phase : phase
    phase := ip > 0 and qu < 0 ? 360 - phase : phase
    dPhase = nz(phase[1]) - phase
    dPhase := nz(phase[1]) < 90 and phase > 270 ? 360 + nz(phase[1]) - phase : dPhase
    dPhase := math.max(math.min(60, dPhase), 1)
    instPeriod = 0.0
    v4 = 0.0
    for i = 0 to hlen by 1
        v4 += nz(dPhase[i])
        if v4 > 360 and instPeriod == 0
            instPeriod := i
            instPeriod
    dcPeriod = 0.0
    dcPeriod := 0.25 * instPeriod + 0.75 * nz(dcPeriod[1])
    dom_in = math.min(math.max(dcPeriod, llen), hlen)
    dom_in

_bpDom(len, bpw) =>
    HP = 0.0
    BP = 0.0
    Peak = 0.0
    Real = 0.0
    counter = 0.0
    DC = 0.0
    alpha2 = (math.cos(0.25 * bpw * 2 * math.pi / len) + math.sin(0.25 * bpw * 2 * math.pi  / len) - 1) / math.cos(0.25 * bpw * 2 * math.pi / len)
    HP := (1 + alpha2 / 2) * (close - nz(close[1])) + (1 - alpha2) * nz(HP[1])
    beta1 = math.cos(2 * math.pi / len)
    gamma1 = 1 / math.cos(2 *  math.pi * bpw / len)
    alpha1 = gamma1 - math.sqrt(gamma1 * gamma1 - 1)
    BP := 0.5 * (1 - alpha1) * (HP - nz(HP[2])) + beta1 * (1 + alpha1) * nz(BP[1]) - alpha1 * nz(BP[2])
    BP := bar_index == 1 or bar_index == 2  ? 0 : BP
    Peak := 0.991 * Peak
    Peak := math.abs(BP) > Peak ? math.abs(BP) : Peak
    Real := Peak != 0 ?  BP / Peak : Real
    DC := nz(DC[1])
    DC := DC < 6 ? 6 : DC
    counter := counter[1] + 1
    if ta.crossover(Real, 0) or ta.crossunder(Real, 0)
        DC := 2 * counter
        if 2 * counter > 1.25 * nz(DC[1])
            DC := 1.25 * DC[1]
        if 2 * counter < 0.8 * nz(DC[1])
            DC := 0.8 * nz(DC[1])
        counter := 0
    DC



auto = _auto_dom(auto_src, auto_min, auto_max, auto_avg)
auto := math.floor(auto) < 1 ? 1 : math.floor(auto)

instant = _instant_cycle(instant_src, instant_min_len, instant_max_len)
instant := math.floor(instant) < 1 ? 1 : math.floor(instant) 

hilbert = _cycleit(hilbert_src, hilbert_len, hilbert_alpha)
hilbert := math.floor(hilbert) < 1 ? 1 : math.floor(hilbert) 

band = _bpDom(bp_period, Bandwidth)
band := math.floor(band) < 1 ? 1 : math.floor(band)

hilbert_dd = _hilber_dual(LPPeriod, min_len_input, max_len_input)
hilbert_dd := math.floor(hilbert_dd) < 1 ? 1 : math.floor(hilbert_dd)



cfb_out = _a_jurik_filt(_a_jurik_filt(_jcfb(cfb_src, cfb_len, cfb_smth), cfb_jper, cfb_phs), cfb_jper, cfb_phs)
cfb_out := math.ceil(cfb_out) < 1 ? 1 : math.ceil(cfb_out)


plot(auto, "Ehlers - Autocorrelation Periodogram Cycle", color = greencolor, linewidth = 2)
plot(instant, "Ehlers - Instantaneous Cycle", color = redcolor, linewidth = 2)
plot(hilbert, "Ehlers - Hilbert Cycle", color = bluecolor, linewidth = 2)
plot(band, "Ehlers - Band-pass Cycle", color = color.yellow, linewidth = 2)
plot(hilbert_dd, "Ehlers - Hilber Dual Differentiator Cycle", color = purplecolor, linewidth = 2)
plot(cfb_out, "Jurik - Composite Fractal Behavior PA Timeframe", color = color.orange, linewidth = 2)






























