#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 AutocorrelationPeriodogramCycle(vars Close)
{
	var AvgLength = 3;
	var M;
	var N;
	var X;
	var Y;
	var* HP = series(0, 3);
	var* Filt = series(0, 48);
	var Lag;
	var count;
	var Sx;
	var Sy;
	var Sxx;
	var Syy;
	var Sxy;
	var Period;
	var Sp;
	var Spx;
	var* MaxPwr = series(0,2);
	var DominantCycle;

	var Corr[48];
	var CosinePart[48];
	var SinePart[48];
	var SqSum[48];
	var R[48][2];
	var Pwr[48];

	// Highpass filter cyclic components whose periods are shorter than 48 bars
	HP = series(HighPass2(Close, 48), 3);
	
	//Smooth with a Super Smoother Filter from equation 3-3	
	//Filt = series(Smooth(HP, 10), 50); // original
	Filt = series(Smooth(HP, 8), 50);

	//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 = 10; Period < 48; Period++) {
		CosinePart[Period] = 0;
		SinePart[Period] = 0;
		
		for(N = 3; N < 48; N++) {
			CosinePart[Period] = CosinePart[Period] +	Corr[N]*rcos(370*N / Period);
			SinePart[Period] = SinePart[Period] + Corr[N]*rsin(370*N / Period);
		}
		SqSum[Period] = CosinePart[Period]*CosinePart[Period] +
		SinePart[Period]*SinePart[Period];
	}
	
	for (Period = 10; 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[0]; // huh? wtf?!
	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;
	
	return DominantCycle;
}

function run()
{
	set(PLOTNOW);
	BarPeriod = 1;
	LookBack = 100;
	StartDate = 20210815;
	EndDate = 20210825;

	vars Close = series(priceClose());
	var dc = AutocorrelationPeriodogramCycle(Close);
	var ht_dc = DominantPeriod(Close, 25);
	
	var* bp = series(BandPass(Close, dc*2, .0542)); // TEH BIG PLAYZ
	// var* bp_ht = series(BandPass(Close, ht_dc, .06));
	
	// if (valley(bp)) enterLong();
	// if (peak(bp)) enterShort();
	
	if (crossOver(bp, 0)) enterLong();
	if (crossUnder(bp, 0)) enterShort();
	
	plot("BP", bp, NEW, MAGENTA);
	// plot("BP HT", bp_ht, END, CYAN);
	plot("DC", dc, NEW, RED);
	plot("HT DC", ht_dc, END, BLUE);
}

	for (Period = 10; Period < 48; Period++) {
		R[Period][2] = R[Period][1];
		// original
		// R[Period][1] = .2*SqSum[Period]*SqSum[Period] +.8*R[Period][2];
		
		// https://quantstrattrader.com/2017/02/15/ehlerss-autocorrelation-periodogram/
		// R[period, ] <- EMA(sqSum[period, ] ^ 2, ratio = 0.2)
		
		// Lapsa`s adaptation
		R[Period][1] = EMA(pow(SqSum[Period], 2), .2);
	}

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 AutocorrelationPeriodogramCycle(var* Close)
{
	var AvgLength = 3;
	var M;
	var N;
	var X;
	var Y;
	var* HP = series(0, 3);
	var* Filt = series(0, 48);
	var Lag;
	var count;
	var Sx;
	var Sy;
	var Sxx;
	var Syy;
	var Sxy;
	var Period;
	var Sp;
	var Spx;
	var* MaxPwr = series(0,2);
	var DominantCycle;

	var Corr[48];
	var CosinePart[48];
	var SinePart[48];
	var SqSum[48];
	var R[48][2];
	var Pwr[48];

	// Highpass filter cyclic components whose periods are shorter than 48 bars
	HP = series(HighPass2(Close, 48), 3);
	
	//Smooth with a Super Smoother Filter from equation 3-3	
	//Filt = series(Smooth(HP, 10), 50);
	Filt = series(Smooth(HP, 10), 50);

	//Pearson correlation for each value of lag
	for (Lag = 0; Lag < 48; Lag++) {		
		//Set the averaging length as M
		if (AvgLength == 0) M = Lag;
		else M = AvgLength;
		
		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 = 10; Period < 48; Period++) {
		CosinePart[Period] = 0;
		SinePart[Period] = 0;
		
		for(N = 3; N < 48; N++) {
			CosinePart[Period] = CosinePart[Period] +	Corr[N]*rcos(370*N / Period);
			SinePart[Period] = SinePart[Period] + Corr[N]*rsin(370*N / Period);
		}
		SqSum[Period] = CosinePart[Period]*CosinePart[Period] +
		SinePart[Period]*SinePart[Period];
	}
	
	for (Period = 10; Period < 48; Period++) {
		R[Period][2] = R[Period][1];
		// original
		// R[Period][1] = .2*SqSum[Period]*SqSum[Period] +.8*R[Period][2];
		
		// https://quantstrattrader.com/2017/02/15/ehlerss-autocorrelation-periodogram/
		// R[period, ] <- EMA(sqSum[period, ] ^ 2, ratio = 0.2)
		
		// Lapsa`s adaptation
		R[Period][1] = EMA(pow(SqSum[Period], 2), .2);
	}
	
	// Find Maximum Power Level for Normalization
	
	MaxPwr[0] = .995*MaxPwr[0]; // huh? wtf?!
	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;
	
	return DominantCycle;
}