//+--------------------------------------------------------------------------------------+
//| Extrapolator.mq4 |
//| Copyright © 2008, gpwr |
//| vlad1004@yahoo.com |
//+--------------------------------------------------------------------------------------+
#property copyright "Copyright © 2008, gpwr"
#property indicator_chart_window
#property indicator_buffers 2
#property indicator_color1 Blue
#property indicator_color2 Red
#property indicator_width1 2
#property indicator_width2 2
//Global constants
#define pi 3.141592653589793238462643383279502884197169399375105820974944592
/*
Method 1: Quinn-Fernandes Algorithm
Method 2: Autocorrelation Method
Method 3: Weighted Burg Method
Method 4: Burg Method with Helme-Nikias weighting function
Method 5: Itakura-Saito (geometric) method
Method 6: Modified covariance method
*/
//Input parameters
extern int Method =1; //Extrapolation method
extern int LastBar =0; //Last bar in the past data
extern int PastBars =300; //Number of past bars
extern double LPOrder =0.6; //Order of linear prediction model; 0 to 1
//LPOrder*PastBars specifies the number of prediction coefficients a[1..LPOrder*PastBars] where a[0]=1
extern int FutBars =100; //Number of bars to predict; for LP is set at PastBars-Order-1
extern int HarmNo =20; //Number of frequencies for Mathod 1; HarmNo=0 computes PastBars harmonics
extern double FreqTOL =0.0001;//Tolerance of frequency calculation for Method 1
//FreqTOL > 0.001 may not converge
extern int BurgWin =0; //Windowing function for Weighted Burg Method; 0=no window 1=Hamming 2=Parabolic
//Indicator buffers
double pv[];
double fv[];
//Global variables
double ETA,INFTY,SMNO;
int np,nf,lb,no,it;
int init()
{
lb=LastBar;
np=PastBars;
no=LPOrder*PastBars;
nf=FutBars;
if(Method>1) nf=np-no-1;
if(HarmNo==0) HarmNo=np;
IndicatorBuffers(4);
SetIndexBuffer(0,pv);
SetIndexBuffer(1,fv);
SetIndexStyle(0,DRAW_LINE,STYLE_SOLID,2);
SetIndexStyle(1,DRAW_LINE,STYLE_SOLID,2);
SetIndexShift(0,-lb);//past data vector 0..np-1; 0 corresponds to bar=lb
SetIndexShift(1,nf-lb);//future data vector i=0..nf; nf corresponds to bar=lb
IndicatorShortName("Extrapolator");
return(0);
}
int deinit(){return(0);}
int start()
{
ArrayInitialize(pv,EMPTY_VALUE);
ArrayInitialize(fv,EMPTY_VALUE);
//Find average of past values
double av=0.0;
for(int i=0;i<np;i++) av+=Open[i+lb];
av/=np;
//Prepare data
if(Method==1)
{
for(i=0;i<np;i++)
{
pv[i]=av;
if(i<=nf) fv[i]=av;
}
}
else
{
double a[],x[];
ArrayResize(a,no+1);
ArrayResize(x,np);
for(i=0;i<np;i++) x[np-1-i]=Open[i+lb]-av;
}
//Select extrapolation method
switch(Method)
{
case 1: //Fit trigomometric series
double w,m,c,s;
for(int harm=0;harm<HarmNo;harm++)
{
Freq(w,m,c,s);
for(i=0;i<np;i++)
{
pv[i]+=m+c*MathCos(w*i)+s*MathSin(w*i);
if(i<=nf) fv[i]+=m+c*MathCos(w*i)-s*MathSin(w*i);
}
}
break;
//Use linear prediction
case 2: ACF(x,no,a);break;
case 3: WBurg(x,no,BurgWin,a);break;
case 4: HNBurg(x,no,a);break;
case 5: Geom(x,no,a);break;
case 6:
bool stop=0;
MCov(x,no,a,stop);
if(stop==1)
{
Print("Terminated prematurely");
return(0);
}
for(i=no;i>=1;i--) a[i]=a[i-1];
}
//Calculate linear predictions
if(Method>1)
{
for(int n=no;n<np+nf;n++)
{
double sum=0.0;
for(i=1;i<=no;i++)
if(n-i<np) sum-=a[i]*x[n-i];
else sum-=a[i]*fv[n-i-np+1];
if(n<np) pv[np-1-n]=sum;
else fv[n-np+1]=sum;
}
fv[0]=pv[0];
for(i=0;i<np-no;i++)
{
pv[i]+=av;
fv[i]+=av;
}
}
//Reorder the predicted vector
for(i=0;i<=(nf-1)/2;i++)
{
double tmp=fv[i];
fv[i]=fv[nf-i];
fv[nf-i]=tmp;
}
return(0);
}
//+--------------------------------------------------------------------------------------+
//Quinn and Fernandes algorithm
void Freq(double& w, double& m, double& c, double& s)
{
double z[],num,den;
ArrayResize(z,np);
double a=0.0;
double b=2.0;
z[0]=Open[lb]-pv[0];
while(MathAbs(a-b)>FreqTOL)
{
a=b;
z[1]=Open[1+lb]-pv[1]+a*z[0];
num=z[0]*z[1];
den=z[0]*z[0];
for(int i=2;i<np;i++)
{
z[i]=Open[i+lb]-pv[i]+a*z[i-1]-z[i-2];
num+=z[i-1]*(z[i]+z[i-2]);
den+=z[i-1]*z[i-1];
}
b=num/den;
}
w=MathArccos(b/2.0);
Fit(w,m,c,s);
return;
}
//+-------------------------------------------------------------------------+
void Fit(double w, double& m, double& c, double& s)
{
double Sc=0.0;
double Ss=0.0;
double Scc=0.0;
double Sss=0.0;
double Scs=0.0;
double Sx=0.0;
double Sxx=0.0;
double Sxc=0.0;
double Sxs=0.0;
for(int i=0;i<np;i++)
{
double cos=MathCos(w*i);
double sin=MathSin(w*i);
Sc+=cos;
Ss+=sin;
Scc+=cos*cos;
Sss+=sin*sin;
Scs+=cos*sin;
Sx+=(Open[i+lb]-pv[i]);
Sxx+=MathPow(Open[i+lb]-pv[i],2);
Sxc+=(Open[i+lb]-pv[i])*cos;
Sxs+=(Open[i+lb]-pv[i])*sin;
}
Sc/=np;
Ss/=np;
Scc/=np;
Sss/=np;
Scs/=np;
Sx/=np;
Sxx/=np;
Sxc/=np;
Sxs/=np;
if(w==0.0)
{
m=Sx;
c=0.0;
s=0.0;
}
else
{
//calculating c, s, and m
double den=MathPow(Scs-Sc*Ss,2)-(Scc-Sc*Sc)*(Sss-Ss*Ss);
c=((Sxs-Sx*Ss)*(Scs-Sc*Ss)-(Sxc-Sx*Sc)*(Sss-Ss*Ss))/den;
s=((Sxc-Sx*Sc)*(Scs-Sc*Ss)-(Sxs-Sx*Ss)*(Scc-Sc*Sc))/den;
m=Sx-c*Sc-s*Ss;
}
return;
}
//+-------------------------------------------------------------------------+
void ACF(double x[], int p, double& a[])
{
int n=ArraySize(x);
double rxx[],r,E,tmp;
ArrayResize(rxx,p+1);
int i,j,k,kh,ki;
//Initialize
for(j=0;j<=p;j++)
{
rxx[j]=0.0;
for(i=j;i<n;i++) rxx[j]+=x[i]*x[i-j];
}
E=rxx[0];
//Main loop
for(k=1;k<=p;k++)
{
//Calculate reflection coefficient
r=-rxx[k];
for(i=1;i<k;i++) r-=a[i]*rxx[k-i];
r/=E;
//Calculate prediction coefficients
a[k]=r;
kh=k/2;
for(i=1;i<=kh;i++)
{
ki=k-i;
tmp=a[i];
a[i]+=r*a[ki];
if(i!=ki) a[ki]+=r*tmp;
}
//Calculate new residual energy
E*=(1-r*r);
}
}
//+-------------------------------------------------------------------------+
void WBurg(double x[], int p, int w, double& a[])
{
int n=ArraySize(x);
double df[],db[];
ArrayResize(df,n);
ArrayResize(db,n);
int i,k,kh,ki;
double tmp,num,den,r;
for(i=0;i<n;i++)
{
df[i]=x[i];
db[i]=x[i];
}
//Main loop
for(k=1;k<=p;k++)
{
//Calculate reflection coefficient
num=0.0;
den=0.0;
if(k==1)
{
for(i=2;i<n;i++)
{
num+=win(i,2,n,w)*x[i-1]*(x[i]+x[i-2]);
den+=win(i,2,n,w)*x[i-1]*x[i-1];
}
r=-num/den/2.0;
if(r>1) r=1.0;
if(r<-1.0) r=-1.0;
}
else
{
for(i=k;i<n;i++)
{
num+=win(i,k,n,w)*df[i]*db[i-1];
den+=win(i,k,n,w)*(df[i]*df[i]+db[i-1]*db[i-1]);
}
r=-2.0*num/den;
}
//Calculate prediction coefficients
a[k]=r;
kh=k/2;
for(i=1;i<=kh;i++)
{
ki=k-i;
tmp=a[i];
a[i]+=r*a[ki];
if(i!=ki) a[ki]+=r*tmp;
}
if(k<p)
//Calculate new residues
for(i=n-1;i>=k;i--)
{
tmp=df[i];
df[i]+=r*db[i-1];
db[i]=db[i-1]+r*tmp;
}
}
}
//+-------------------------------------------------------------------------+
double win(int i, int k, int n, int w)
{
if(w==0) return(1.0);
if(w==1) return(0.54-0.46*MathCos(pi*(2.0*(i-k)+1.0)/(n-k)));
if(w==2) return(6.0*(i-k+1.0)*(n-i)/(n-k)/(n-k+1.0)/(n-k+2.0));
}
//+-------------------------------------------------------------------------+
void HNBurg(double x[], int p, double& a[])
{
int n=ArraySize(x);
double df[],db[];
ArrayResize(df,n);
ArrayResize(db,n);
int i,k,kh,ki;
double w,tmp,num,den,r;
for(i=0;i<n;i++)
{
df[i]=x[i];
db[i]=x[i];
}
//Main loop
for(k=1;k<=p;k++)
{
//Calculate reflection coefficient
num=0.0;
den=0.0;
if(k==1)
{
for(i=2;i<n;i++)
{
w=x[i-1]*x[i-1];
num+=w*x[i-1]*(x[i]+x[i-2]);
den+=w*x[i-1]*x[i-1];
}
r=-num/den/2.0;
if(r>1) r=1.0;
if(r<-1.0) r=-1.0;
}
else
{
w=0.0;
for(i=1;i<k;i++) w+=x[i]*x[i];
for(i=k;i<n;i++)
{
num+=w*df[i]*db[i-1];
den+=w*(df[i]*df[i]+db[i-1]*db[i-1]);
w=w+x[i]*x[i]-x[i-k+1]*x[i-k+1];
}
r=-2.0*num/den;
}
//Calculate prediction coefficients
a[k]=r;
kh=k/2;
for(i=1;i<=kh;i++)
{
ki=k-i;
tmp=a[i];
a[i]+=r*a[ki];
if(i!=ki) a[ki]+=r*tmp;
}
if(k<p)
//Calculate new residues
for(i=n-1;i>=k;i--)
{
tmp=df[i];
df[i]+=r*db[i-1];
db[i]=db[i-1]+r*tmp;
}
}
}
//+-------------------------------------------------------------------------+
void Geom(double x[], int p, double& a[])
{
int n=ArraySize(x);
double df[],db[];
ArrayResize(df,n);
ArrayResize(db,n);
int i,k,kh,ki;
double tmp,num,denf,denb,r;
for(i=0;i<n;i++)
{
df[i]=x[i];
db[i]=x[i];
}
//Main loop
for(k=1;k<=p;k++)
{
//Calculate reflection coefficient
num=0.0;
denf=0.0;
denb=0.0;
for(i=k;i<n;i++)
{
num+=df[i]*db[i-1];
denf+=df[i]*df[i];
denb+=db[i-1]*db[i-1];
}
r=-num/MathSqrt(denf)/MathSqrt(denb);
//Calculate prediction coefficients
a[k]=r;
kh=k/2;
for(i=1;i<=kh;i++)
{
ki=k-i;
tmp=a[i];
a[i]+=r*a[ki];
if(i!=ki) a[ki]+=r*tmp;
}
if(k<p)
//Calculate new residues
for(i=n-1;i>=k;i--)
{
tmp=df[i];
df[i]+=r*db[i-1];
db[i]=db[i-1]+r*tmp;
}
}
}
//+-------------------------------------------------------------------------+
/* Modified Covariance method from Marple's book. It solves
e[j]=y[j]+SUM(a[i]*y[j-i-1],i=0...m-1), j=0..n-1
for a[] by the least-squares minimization of e[j] in both directions of j.
The code substitues y[j] with x[n-1-j] because, in the provided data,
x[0] is the latest data point. */
void MCov(double x[], int ip, double& a[], bool& stop)
{
//Initialization
int n=ArraySize(x);
double c[],d[],r[],v;
ArrayResize(c,ip+1);
ArrayResize(d,ip+1);
ArrayResize(r,ip);
int k,m,mk;
double r1,r2,r3,r4,r5,delta,gamma,lambda,theta,psi,xi;
double save1,save2,save3,save4,c1,c2,c3,c4,ef,eb;
r1=0.0;
for(k=1;k<n-1;k++) r1+=2.0*x[k]*x[k];
r2=x[n-1]*x[n-1];
r3=x[0]*x[0];
r4=1.0/(r1+2.0*(r2+r3));
v=r1+r2+r3;
delta=1.0-r2*r4;
gamma=1.0-r3*r4;
lambda=x[0]*x[n-1]*r4;
c[0]=x[0]*r4;
d[0]=x[n-1]*r4;
//Main loop
for(m=0;;m++)
{
save1=0.0;
for(k=m+1;k<n;k++) save1+=x[n-1-k]*x[n-k+m];
save1*=2.0;
r[m]=save1;
theta=x[0]*d[0];
psi=x[0]*c[0];
xi=x[n-1]*d[0];
if(m>0)
{
for(k=1;k<=m;k++)
{
theta+=x[k]*d[k];
psi+=x[k]*c[k];
xi+=x[n-1-k]*d[k];
r[k-1]-=(x[m]*x[m-k]+x[n-1-m]*x[n-1-m+k]);
save1+=r[k-1]*a[m-k];
}
}
//order update of a[]
c1=-save1/v;
a[m]=c1;
v*=(1.0-c1*c1);
if(m>0)
{
for(k=0;k<(m+1)/2;k++)
{
mk=m-k-1;
save1=a[k];
a[k]=save1+c1*a[mk];
if(k!=mk) a[mk]+=c1*save1;
}
}
if(m==ip-1)
{
v*=(0.5/(n-1-m));
break;
}
//time update of c[],d[],gamma,delta, and lambda
r1=1.0/(delta*gamma-lambda*lambda);
c1=(theta*lambda+psi*delta)*r1;
c2=(psi*lambda+theta*gamma)*r1;
c3=(xi*lambda+theta*delta)*r1;
c4=(theta*lambda+xi*gamma)*r1;
for(k=0;k<=m/2;k++)
{
mk=m-k;
save1=c[k];
save2=d[k];
save3=c[mk];
save4=d[mk];
c[k]+=(c1*save3+c2*save4);
d[k]+=(c3*save3+c4*save4);
if(k!=mk)
{
c[mk]+=(c1*save1+c2*save2);
d[mk]+=(c3*save1+c4*save2);
}
}
r2=psi*psi;
r3=theta*theta;
r4=xi*xi;
r5=gamma-(r2*delta+r3*gamma+2.0*psi*lambda*theta)*r1;
r2=delta-(r3*delta+r4*gamma+2.*theta*lambda*xi)*r1;
gamma=r5;
delta=r2;
lambda+=(c3*psi+c4*theta);
if(v<=0.0)
{
Print("Error: negative or zero variance in MCov");
stop=true;
return;
}
if(delta<=0.0 || delta>1.0 || gamma<=0.0 || gamma>1.0)
{
Print("Error: delta and gamma are outside (0,1) in MCov");
stop=true;
return;
}
//time update of a[]; order updates of c[],d[],gamma,delta, and lambda
r1=1.0/v;
r2=1.0/(delta*gamma-lambda*lambda);
ef=x[n-m-2];
eb=x[m+1];
for(k=0;k<=m;k++)
{
ef+=a[k]*x[n-1-m+k];
eb+=a[k]*x[m-k];
}
c1=eb*r1;
c2=ef*r1;
c3=(eb*delta+ef*lambda)*r2;
c4=(ef*gamma+eb*lambda)*r2;
for(k=m;k>=0;k--)
{
save1=a[k];
a[k]=save1+c3*c[k]+c4*d[k];
c[k+1]=c[k]+c1*save1;
d[k+1]=d[k]+c2*save1;
}
c[0]=c1;
d[0]=c2;
r3=eb*eb;
r4=ef*ef;
v-=(r3*delta+r4*gamma+2.0*ef*eb*lambda)*r2;
delta-=r4*r1;
gamma-=r3*r1;
lambda+=ef*eb*r1;
if(v<=0.0)
{
Print("Error: negative or zero variance in MCov");
stop=true;
return(0);
}
if(delta<=0.0 || delta>1.0 || gamma<=0.0 || gamma>1.0)
{
Print("Error: delta and gamma are outside (0,1) in MCov");
stop=true;
return(0);
}
}
//for(m=ip;m>=1;m--) a[m]=a[m-1];
//a[0]=1.0;
}
Sample
Analysis
Market Information Used:
Series array that contains open prices of each bar
Indicator Curves created:
Implements a curve of type DRAW_LINE
Indicators Used:
Custom Indicators Used:
Order Management characteristics:
Other Features: