#include <profile.c>

var rad(var degrees) { return degrees*PI/180; }
var rcos(var degrees) { return cos(rad(degrees)); }
var rsin(var degrees) { return sin(rad(degrees)); }

var AdaptiveStochastic(vars Close)
{
	var AvgLength = 3;
	var M;
	var N;
	var X;
	var Y;
	var alpha1;
	var* HP = series(0,3);
	var a1;
	var b1;
	var c1;
	var c2;
	var c3;
	var* Filt = series(0,3);
	var Lag;
	var count;
	var Sx;
	var Sy;
	var Sxx;
	var Syy;
	var Sxy;
	var Period;
	var Sp;
	var Spx;
	var* MaxPwr = series(0,3);
	var DominantCycle;
	
	var* Corr = series(0, 48);
	var* CosinePart = series(0, 48);
	var* SinePart = series(0, 48);
	var* SqSum = series(0, 48);
	var R[48][2];
	var Pwr[48];	
	
	//Highpass filter cyclic components whose periods are shorter than 48 bars
	alpha1 = (rcos(.707*360 / 48) + rsin(.707*360 / 48) - 1) / rcos(.707*360 / 48);
	
	HP[0] = (1 - alpha1 / 2)*(1 - alpha1 / 2)*
		(Close[0] - 2*Close[1] + Close[2]) + 
		2*(1 - alpha1)*HP[1] - 
		(1 - alpha1) * (1 - alpha1)*HP[2];
	
	
	//Smooth with a Super Smoother Filter from equation 3-3
	a1 = exp(-1.414*3.14159 / 10);
	b1 = 2*a1*rcos(1.414*180 / 10);
	c2 = b1;
	c3 = -a1*a1;
	c1 = 1 - c2 - c3;
	Filt[0] = c1*(HP[0] + HP[1]) / 2 + c2*Filt[1] + c3*Filt[2];
		
	//Pearson correlation for each value of lag
	for (Lag = 0; Lag < 48; Lag++) {		
		//Set the averaging length as M
		M = AvgLength;
		if (AvgLength == 0) M = Lag;
		Sx = 0;
		Sy = 0;
		Sxx = 0;
		Syy = 0;
		Sxy = 0;
		
		for (count = 0; count < M - 1; count++) {
			X = Filt[count];
			Y = Filt[Lag + count];
			Sx = Sx + X;
			Sy = Sy + Y;
			Sxx = Sxx + X*X;
			Sxy = Sxy + X*Y;
			Syy = Syy + Y*Y;
		}
		
		if ( (M*Sxx - Sx*Sx)*(M*Syy - Sy*Sy) > 0 ) {
			Corr[Lag] = (M*Sxy - Sx*Sy)/sqrt((M*Sxx-Sx*Sx)*(M*Syy-Sy*Sy));
		}
	}
	
	for (Period = 0; Period < 48; Period++) {
		CosinePart[Period] = 0;
		SinePart[Period] = 0;
		
		for(N = 3; N < 48; N++) {
			CosinePart[Period] = CosinePart[Period] + Corr[N]*rcos(360*N / Period);
			SinePart[Period] = SinePart[Period] + Corr[N]*rsin(360*N / Period);
		}
		SqSum[Period] = CosinePart[Period]*CosinePart[Period] +
		SinePart[Period]*SinePart[Period];
	}
	
	for (Period = 0; Period < 48; Period++) {
		R[Period][2] = R[Period][1];
		R[Period][1] = .2*SqSum[Period]*SqSum[Period] +.8*R[Period][2];
	}
	
	// Find Maximum Power Level for Normalization
	MaxPwr[0] = .995*MaxPwr[1];
	for (Period = 10; Period < 48; Period++) {
		if (R[Period][1] > MaxPwr[0]) MaxPwr[0] = R[Period][1];
	}
	
	for (Period = 3; Period < 48; Period++) {
		Pwr[Period] = R[Period][1] / MaxPwr[0];
	}
	
	//Compute the dominant cycle using the CG of the spectrum
	Spx = 0;
	Sp = 0;
	for(Period = 10; Period < 48; Period++) {
		if (Pwr[Period] >= .5) {
			Spx = Spx + Period*Pwr[Period];
			Sp = Sp + Pwr[Period];
		}
	}
	
	if (Sp != 0) DominantCycle = Spx / Sp;
	if (DominantCycle < 10) DominantCycle = 10;
	if (DominantCycle > 48) DominantCycle = 48;
	
	//Stochastic Computation starts here
	var* HighestC = series(0, 3);
	var* LowestC = series(0, 3);
	var* Stoc = series(0, 3);
	var* SmoothNum = series(0, 3);
	var* SmoothDenom = series(0, 3);
	var* AdaptiveStochastic = series(0, 3);
	
	HighestC[0] = Filt[0];
	LowestC[0] = Filt[0];
	
	for (count = 0; count < DominantCycle - 1; count++ ) {
		if (Filt[count] > HighestC[0]) HighestC[0] = Filt[count];
		if (Filt[count] < LowestC[0]) LowestC[0] = Filt[count];	
	}
	
	Stoc[0] = (Filt[0] - LowestC[0]) / (HighestC[0] - LowestC[0]);
	return AdaptiveStochastic[0] = 
		c1*(Stoc[0] + Stoc[1]) / 2 + c2*AdaptiveStochastic[1] + c3*AdaptiveStochastic[2];
}

function run()
{
	set(PLOTNOW);
	BarPeriod = 1;
	LookBack = 500;
	StartDate = 20210809;
	EndDate = 20210810;

	vars closes = series(priceClose(), 48);
	var astoch = AdaptiveStochastic(closes);
	
	var original = StochEhlers(closes, 48, 10, 48);
	
	plot("ASTOCH", astoch, NEW, RED);
        plot("Zorro`s EStoch", original, NEW, BLUE);
}