HurstExponent
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © nolantait

//@version=5

// @description Library to calculate Hurst Exponent refactored from Hurst Exponent - Detrended Fluctuation Analysis [pig]
library("HurstExponent")
import nolantait/Moments/1 as moments

// @function Calculates a series subtracted from the series mean.
// @param src The series used to calculate the difference from the mean (e.g. log returns).
// @returns The series subtracted from the series mean
export demean(float src, int length) =>
    src - moments.mean(src, length)

// @function Calculates a cumulated sum from the series.
// @param src The series used to calculate the cumulative sum (e.g. demeaned log returns).
// @param length The length used to calculate the cumulative sum (e.g. 100).
// @returns The cumulative sum of the series as an array
export cumsum(float src, int length) =>
	csum = array.new_float(length)
	sum = 0.0
	for i=0 to length-1
	    sum := src[i] + sum
	    array.set(csum, i, sum)
	csum

// @function Calculates an aproximated log scale. Used to save sample size
// @param scale The scale to aproximate.
// @param length The length used to aproximate the expected scale.
// @returns The aproximated log scale of the value
export aproximateLogScale(int scale, int length, int minScale = 8, int maxScale = 2) =>
    aproximation = math.pow(length / (minScale * maxScale), 0.1111111111)
    math.round(minScale * math.pow(aproximation, scale))


// @function Calculates linear trend to determine error between linear trend and cumulative sum
// @param cumulativeSum The cumulative sum array to regress.
// @param barId The barId for the slice
// @param numberOfSegments The total number of segments used for the regression calculation
// @returns The error between linear trend and cumulative sum
export rootMeanSum(float[] cumulativeSum, int barId, int numberOfSegements = 10) =>
    sequence = array.new_float(numberOfSegements)
    count = 0
    for i=0 to numberOfSegements-1
        count := count + 1
        array.set(sequence, i, count)
    y = array.slice(cumulativeSum, barId, barId + numberOfSegements)
    xStandardDeviation = array.stdev(sequence) * math.sqrt(numberOfSegements / (numberOfSegements - 1))
    yStandardDeviation = array.stdev(y) * math.sqrt(numberOfSegements / (numberOfSegements - 1))
    covariance = array.covariance(sequence, y) * (numberOfSegements / (numberOfSegements - 1))
    regression = math.pow(covariance / (xStandardDeviation * yStandardDeviation), 2)
    math.sqrt(1 - regression) * yStandardDeviation


// @function Calculates the Root Mean Sum Measured for each block (e.g the aproximated log scale)
// @param cumulativeSum The cumulative sum array to regress and determine the average of.
// @param barId The barId for the slice
// @param length The length used for finding the average
// @returns The average root mean sum error of the cumulativeSum
export averageRootMeanSum(float[] cumulativeSum, int barId, int length) =>
    calculatedLength = math.floor(length / barId)
    sum = 0.0
    for i=0 to calculatedLength-1
        sum := sum + rootMeanSum(cumulativeSum, i * barId, barId)
    math.log10(sum / calculatedLength)


// @function Calculates the critical values for a hurst exponent for a given length
// @param length The length used for finding the average
// @returns The critical value, upper critical value and lower critical value for a hurst exponent
export criticalValues(int length) =>
    critical = 1.645 * (0.3912 / math.pow(length, 0.3))
    upper = 0.5 + critical
    lower = 0.5 - critical
    [critical, upper, lower]

// @function Calculates the hurst exponent slope measured from root mean sum, scaled to log log plot using linear regression
// @param cumulativeSum The cumulative sum array to regress and determine the average of.
// @param length The length used for the hurst exponent sample size
// @returns The slope of the hurst exponent
export slope(float[] cumulativeSum, int length) =>
    fluxes = array.new_float(10)
    scales = array.new_float(10)
    for i=0 to 9
        logScale = aproximateLogScale(i, length)
        array.set(fluxes, i, averageRootMeanSum(cumulativeSum, logScale, length))
        array.set(scales, i, math.log10(logScale))
    array.covariance(scales, fluxes) / array.variance(scales)

// @function Smooths input using advanced linear regression
// @param src The series to smooth (e.g. hurst exponent slope)
// @param length The length used to smooth
// @returns The src smoothed according to the given length
export smooth(float src, int length) =>
    pi = 2.0 * math.asin(1.0)
    sqrt2 = math.sqrt(2.0)
    lambda = pi * sqrt2 / length
    a1 = math.exp(-lambda)
    coeffecient2 = 2.0 * a1 * math.cos(lambda)
    coeffecient3 = -math.pow(a1, 2.0)
    coeffecient1 = 1.0 - coeffecient2 - coeffecient3
    filter = 0.0
    filter := coeffecient1 * (src + nz(src[1])) * 0.5 + coeffecient2 * nz(filter[1]) + coeffecient3 * nz(filter[2])
    filter

// @function Wrapper function to calculate the hurst exponent slope
// @param src The series used for returns calculation (e.g. close)
// @param hurstLength The length used to calculate the hurst exponent (should be greater than 50)
// @returns The src smoothed according to the given length
export exponent(float src, int hurstLength = 100) =>
    returns = moments.logReturns(src)
    demeanedReturns = demean(returns, hurstLength)
    cumulativeSum = cumsum(demeanedReturns, hurstLength)
    slope(cumulativeSum, hurstLength)

// Example
hurstExponent = exponent(close)
smoothedHurstExponent = smooth(hurstExponent, 18)
plot(hurstExponent, "Hurst Exponent")
plot(smoothedHurstExponent, "Smoothed Hurst Exponent", color=color.yellow)