John Ehlers' Autocorrelation Periodogram

{
Autocorrelation Periodogram
© 2013	John F. Ehlers
}

Vars:
AvgLength(3), M(0),
N(0),	 
X(0),
Y(0),
alpha1(0), HP(0),
a1(0),
b1(0),
c1(0),
c2(0),
c3(0),
Filt(0),
Lag(0),
count(0),
Sx(0),
Sy(0),
Sxx(0),
Syy(0),
Sxy(0),
Period(0), Sp(0),
Spx(0),
MaxPwr(0), DominantCycle(0), Color1(0), Color2(0), Color3(0);

Arrays:
Corr[48](0),
CosinePart[48](0), SinePart[48](0),
SqSum[48](0),
R[48, 2](0),
Pwr[48](0);

//Highpass filter cyclic components whose periods are shorter than 48 bars
alpha1 =  (Cosine(.707*360 / 48) +  Sine (.707*360 / 48) - 1) / Cosine(.707*360 / 48);
HP =  (1 - alpha1 / 2)*(1 - alpha1 / 2)*(Close - 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 =  expvalue(-1.414*3.14159 / 10); b1 =  2*a1*Cosine(1.414*180 / 10); c2 =  b1;
c3 =  -a1*a1;
c1 =  1 - c2 - c3;
Filt =  c1*(HP +  HP[1]) / 2 +  c2*Filt[1] +  c3*Filt[2];
 

//Pearson correlation for each value of lag For Lag =  0 to 48 Begin
//Set the averaging length as M M =  AvgLength;
If AvgLength =  0 Then M =  Lag;
Sx =  0;
Sy =  0;
Sxx =  0;
Syy =  0;
Sxy =  0;
For count =  0 to M - 1 Begin 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;
End;
If (M*Sxx - Sx*Sx)*(M*Syy - Sy*Sy) > 0 Then Corr[Lag] =
(M*Sxy - Sx*Sy)/SquareRoot((M*Sxx - Sx*Sx)*(M*Syy - Sy*Sy)); End;
For Period =  10 to 48 Begin CosinePart[Period] =  0;
SinePart[Period] =  0; For N =  3 to 48 Begin
CosinePart[Period] =   CosinePart[Period] +
Corr[N]*Cosine(370*N / Period);
SinePart[Period] = SinePart[Period] + Corr[N]*Sine(370*N / Period);
End;
SqSum[Period] =    CosinePart[Period]*CosinePart[Period] +
SinePart[Period]*SinePart[Period];
End;	(Continued )
 
For Period =  10 to 48 Begin R[Period, 2] =  R[Period, 1];
R[Period, 1] =  .2*SqSum[Period]*SqSum[Period] +
.8*R[Period, 2]; End;

//Find Maximum Power Level for Normalization MaxPwr =  .995*MaxPwr;
For Period =  10 to 48 Begin
If R[Period, 1] > MaxPwr Then MaxPwr =  R[Period, 1]; End;
For Period =  3 to 48 Begin
Pwr[Period] =  R[Period, 1] / MaxPwr; End;
 
//Compute the dominant cycle using the CG of the spectrum Spx =  0;
Sp =  0;
For Period =  10 to 48 Begin
If Pwr[Period] >=  .5 Then Begin Spx =  Spx +  Period*Pwr[Period]; Sp =  Sp +  Pwr[Period];
End; End;
If Sp <> 0 Then DominantCycle =  Spx / Sp;
Plot2(DominantCycle, “DC”, RGB(0, 0, 255), 0, 2);

{
//Increase Display Resolution by raising the NormPwr to a higher mathematical power (optional)
For Period =  10 to 48 Begin Pwr[Period] =  Power(Pwr[Period], 2);
End;
}

function run()
{
//=====John Ehlers' Autocorrelation Periodogram=====

//Roofing Filter used on prices
vars RoofPrice = series(Roof(Price, 10, 40)); //standard roofing filter
vars RoofAGC = series(AGC(RoofPrice, 1)); //output normalized using AGC algorithm

//Ehlers Autocorrelation Periodogram

var corr[48]; //an array that will store the correlation coefficient for each value of lag from 1 to 48
for (i = 1; i < 49; i++)
	corr[i] = Correlation(RoofAGC, RoofAGC+i, i); //assign a correlation coefficent to each value of lag. 

//discrete Fourier transform of autcorrelation results
}

//===== Ehlers Autocorrelation Periodogram =====

function run()
{
vars Price = series(price());

// Roofing filter applied to price series
var LowCutoff = 10;
var HighCutoff = 48;
vars RoofPrice = series(Roof(Price, LowCutoff, HighCutoff));
vars RoofPriceAGC = series(AGC(RoofPrice, 1)); //normalize output with automatic gain control

// Pearson correlation for each value of lag
var corr[48]; //declare an array to store the correlation coefficients for each value of lag up to the HighCutoff
int lag;
int AvgLength = 0; //Averaging length for the correlation algorithm. Ehlers recommends a value of 3.
int M; //used to set AvgLength to the lag value if AvgLength is not set
for (lag = 0; lag <= HighCutoff; lag++)
	{
		if (AvgLength == 0)
			M = lag;
		else M = AvgLength;
	
		corr[lag] = Correlation(RoofPriceAGC, RoofPriceAGC+lag, M); //assign a correlation coefficient for each value of lag to the array corr[]
	}
	
//Discrete Fourier transform. Correlate autocorrelation values with the cosine and sine of each period of interest. The sum of the squares of each value represents relative power at each period
var cosinePart[38]; //array for storing the cosine values
var sinePart[38]; //array for storing the sine values
var sqSum[38]; //array for storing the squared and summed values of the sine and cosine parts
int period;
int N;
for (period = LowCutoff; period <= HighCutoff; period++)
	{
		for (N = 3; N <= HighCutoff; N++)
		{
			cosinePart[period] = cosinePart[period] + corr[N]*cos((370*N*PI/180)/period);
			sinePart[period] = sinePart[period] + corr[N]*sin((370*N*PI/180)/period);
		}
		sqSum[period] = cosinePart[period]*cosinePart[period] + sinePart[period]*sinePart[period];
	}
	//Smooth power measurements using an EMA
var R[2][35];
for (period = LowCutoff; period <= HighCutoff; period++)
	{
		R[1][period] = R[0][period];
		R[0][period] = 0.2*sqSum[period]*sqSum[period] + 0.8*R[1][period];
	}
//Apply automatic gain control with decay factor 0.995
vars maxPwr = series(0);
maxPwr[0] = 0.995*maxPwr[1];
for (period = LowCutoff; period <= HighCutoff; period++)
	{
		if (R[0][period] > maxPwr[0])
			maxPwr[0] = (R[0][period]);
	}
	
var pwr[45]; 
for (period = 3; period < HighCutoff; period++)
	{
		pwr[period] = R[0][period] / maxPwr[0];
	}	

//Compute dominant cycle using centre of gravity algorithm
var Spx = 0;
var Sp = 0;
var DomPer; //Dominant cycle period
for (period = LowCutoff; period <= HighCutoff; period++)
	{
		if (pwr[period] >= 0.5)
			{
				Spx = Spx + period*pwr[period];
				Sp = Sp + pwr[period];
			}
	}
if (Sp != 0)
	DomPer = Spx/Sp;	
	
plot("DominantCyclePeriod", DomPer, NEW, BLUE);
	
	
}