Hurst Exponent (Dubuc's variation method)
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// Copyright ©2022 and_then_it_CRASHED.

// This is an implementation of Dubuc's variation method for estimating the Hurst exponent.
// Dubuc B, Quiniou JF, Roques-Carmes C, Tricot C. Evaluating the fractal dimension of profiles. Physical Review A. 1989;39(3):1500-1512. DOI: 10.1103/PhysRevA.39.1500
// https://www.intechopen.com/chapters/64463

//@version=5

library("Hurst")

import lejmer/OrdinaryLeastSquares/2

// @function Estimate the Hurst Exponent using Dubuc's variation method
// @param length The length of the history window to use. Large values do not cause lag.
// @param samples The number of scale samples to take within the window. These samples are then used for regression. The minimum value is 2 but 3+ is recommended. Large values give more accurate results but suffer from a performance penalty.
// @param hi The high value of the series to analyze.
// @param lo The low value of the series to analyze.
export hurst(int length=100, int samples = 5, float hi=high, float lo=low) =>

	if length < 2
		runtime.error('Dubuc.hurst() window length must be at least 2.')
	if samples < 2
		runtime.error('Dubuc.hurst() must have at least two samples.')

	var x = matrix.new<float>(samples, 2, 1) // the second column is full of 1's for the intercept
	var y = array.new<float>(samples)
	min_iep = math.exp(-1)  // inverse epsilon
	max_iep = math.exp(-length)
	float base_count = na
	string msg = ''
	for s = 0 to samples-1
		iep = math.log(length) * s / (samples-1)
		ep = int(math.exp(iep))
		t = ta.highest(hi, 2*ep+1) // neighborhood top
		b = ta.lowest(lo, 2*ep+1)  // neighborhood bottom
		count = 0
		area = 0.
		for i = 0 to length
			if not na(t) and not na(b)
				count += 1
				area += t - b
		if na(base_count) // this will be set on the first window, the smallest, which will have the highest count
			base_count := count
		area *= base_count / count // normalize areas due to different sample sizes
		msg := str.format('{0} {1}=>{2}*{3}/{4}',msg,ep,area,count,base_count)
		matrix.set(x, s, 0, math.log(1/ep)) // x value
		array.set(y, s, math.log(area/ep/ep))

	fit = OrdinaryLeastSquares.solve(x,y)
	slope = array.get(fit, 0)
	hurst = slope - 1
	hurst

length = input.int(100, 'Length', minval = 2, tooltip = 'The window size for analysis. Using a large window does not cause lag and is not detrimental to getting a good measurement.')
samples = input.int(5, 'Samples', minval = 2, tooltip = 'Within the window, this many samples of different "scales" are used for regression. Higher numbers give more accurate results but cause performance delays. At least 3 samples are recommended.')

ind = hurst(length, samples=samples)
plot(ind)
hline(0.5, color=color.new(color.white,70))
