Autocorrelation Periodogram Dominant Cycle

I think I fixed bugs. At least most of them.


Spikes looks suspicious.

And I still don't understand what's that MaxPwr decay line is supposed to do.
Perhaps some odd EasyLanguage behavior?

I mean - decay doesn't get applied. MaxPwr will always get initialized as zero.

Some room for refactoring: Correlation, EMA, AGC should be able to replace some lines.

Still a bit sluggish but might be actually usable.

Seems to be more accurate than calculating dominant period via Hilbert transform (assuming that it works correctly).
Combined with BandPass, some curve fitting and decreased smoothing - it seems like a reasonable strategy foundation.

Code
#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);
}


[Linked Image]

Quote

Monte Carlo Analysis... Median AR 435%
Win 0.23$ MI 0.69$ DD 0.15$ Capital 1.61$
Trades 394 Win 47.5% Avg +5.8p Bars 28
AR 512% PF 1.12 SR 0.00 UI 0% R2 1.00


^ 10 days



no moar spikey:


Code
	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);
	}


Code

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;
}


[Linked Image]

Failed to use Correlation function successfully.
I mean - I can pinpoint what should be replaced yet failing to do that cause of LiteC and my stupidity.

Last edited by Lapsa; 09/08/21 10:21.