# Future Prediction System - Amibroker AFL Code

### Click Image To Enlarge. Please Rate And Comment.

```EnableScript("VBScript");
<%
' Called by ARG_Burg AFL funct. to compute Ki in Double Precision
function Compute_KI(f, b, W, i, LongBar)
n = LongBar - 1

for j = i to n
Num = Num + cDBL(W(j))*( cDBL(f(j))*cDBL(b(j-1)) )
Den = Den + cDBL(W(j))*( cDBL(f(j))^2 + cDBL(b(j-1))^2 )
next

KI = 2*Num/Den
Compute_KI = KI

end function

function Gaussian_Elimination (GE_Order, GE_N, GE_SumXn, GE_SumYXn)
Dim b(10, 10)
Dim w(10)
Dim Coeff(10)

for i = 1 To 10
Coeff(i) = 0
next

n = GE_Order + 1

for i = 1 to n
for j = 1 to  n
if i = 1 AND j = 1 then
b(i, j) = cDBL(GE_N)
else
b(i, j) = cDbl(GE_SumXn(i + j - 2))
end if
next
w(i) = cDbl(GE_SumYXn(i))
next

n1 = n - 1
for i = 1 to n1
big = cDbl(abs(b(i, i)))
q = i
i1 = i + 1

for j = i1 to n
ab = cDbl(abs(b(j, i)))
if (ab >= big) then
big = ab
q = j
end if
next

if (big <> 0.0) then
if (q <> i) then
for j = 1 to n
Temp = cDbl(b(q, j))
b(q, j) = b(i, j)
b(i, j) = Temp
next
Temp = w(i)
w(i) = w(q)
w(q) = Temp
end if
end if

for j = i1 to n
t = cDbl(b(j, i) / b(i, i))
for k = i1 to n
b(j, k) = b(j, k) - t * b(i, k)
next
w(j) = w(j) - t * w(i)
next
next

if (b(n, n) <> 0.0) then

Coeff(n) = w(n) / b(n, n)
i = n - 1

while i > 0
SumY = cDbl(0)
i1 = i + 1
for j = i1 to n
SumY = SumY + b(i, j) * Coeff(j)
next
Coeff(i) = (w(i) - SumY) / b(i, i)
i = i - 1
wend

Gaussian_Elimination = Coeff

end if
end function
%>

function PolyFit(GE_Y, GE_BegBar, GE_EndBar, GE_Order, GE_ExtraB, GE_ExtraF)
{
BI        = BarIndex();

GE_N      = GE_EndBar - GE_BegBar + 1;
GE_XBegin = -(GE_N - 1) / 2;
GE_X      = IIf(BI < GE_BegBar, 0, IIf(BI > GE_EndBar, 0, (GE_XBegin + BI - GE_BegBar)));

GE_X_Max  = LastValue(Highest(GE_X));

GE_X      = GE_X / GE_X_Max;

X1 = GE_X;

GE_Y      = IIf(BI < GE_BegBar, 0, IIf(BI > GE_EndBar, 0, GE_Y));

GE_SumXn  = Cum(0);
GE_SumXn   = LastValue(Cum(GE_X));
GE_X2     = GE_X * GE_X;     GE_SumXn   = LastValue(Cum(GE_X2));
GE_X3     = GE_X * GE_X2;    GE_SumXn   = LastValue(Cum(GE_X3));
GE_X4     = GE_X * GE_X3;    GE_SumXn   = LastValue(Cum(GE_X4));
GE_X5     = GE_X * GE_X4;    GE_SumXn   = LastValue(Cum(GE_X5));
GE_X6     = GE_X * GE_X5;    GE_SumXn   = LastValue(Cum(GE_X6));
GE_X7     = GE_X * GE_X6;    GE_SumXn   = LastValue(Cum(GE_X7));
GE_X8     = GE_X * GE_X7;    GE_SumXn   = LastValue(Cum(GE_X8));
GE_X9     = GE_X * GE_X8;    GE_SumXn   = LastValue(Cum(GE_X9));
GE_X10    = GE_X * GE_X9;    GE_SumXn  = LastValue(Cum(GE_X10));
GE_X11    = GE_X * GE_X10;   GE_SumXn  = LastValue(Cum(GE_X11));
GE_X12    = GE_X * GE_X11;   GE_SumXn  = LastValue(Cum(GE_X12));
GE_X13    = GE_X * GE_X12;   GE_SumXn  = LastValue(Cum(GE_X13));
GE_X14    = GE_X * GE_X13;   GE_SumXn  = LastValue(Cum(GE_X14));
GE_X15    = GE_X * GE_X14;   GE_SumXn  = LastValue(Cum(GE_X15));
GE_X16    = GE_X * GE_X15;   GE_SumXn  = LastValue(Cum(GE_X16));

GE_SumYXn = Cum(0);
GE_SumYXn  = LastValue(Cum(GE_Y));
GE_YX     = GE_Y    * GE_X;  GE_SumYXn  = LastValue(Cum(GE_YX));
GE_YX2    = GE_YX   * GE_X;  GE_SumYXn  = LastValue(Cum(GE_YX2));
GE_YX3    = GE_YX2  * GE_X;  GE_SumYXn  = LastValue(Cum(GE_YX3));
GE_YX4    = GE_YX3  * GE_X;  GE_SumYXn  = LastValue(Cum(GE_YX4));
GE_YX5    = GE_YX4  * GE_X;  GE_SumYXn  = LastValue(Cum(GE_YX5));
GE_YX6    = GE_YX5  * GE_X;  GE_SumYXn  = LastValue(Cum(GE_YX6));
GE_YX7    = GE_YX6  * GE_X;  GE_SumYXn  = LastValue(Cum(GE_YX7));
GE_YX8    = GE_YX7  * GE_X;  GE_SumYXn  = LastValue(Cum(GE_YX8));

GE_Coeff  = Cum(0);

GE_VBS    = GetScriptObject();
GE_Coeff  = GE_VBS.Gaussian_Elimination(GE_Order, GE_N, GE_SumXn, GE_SumYXn);

for (i = 1; i <= GE_Order + 1; i++)
printf(NumToStr(i, 1.0) + " = " + NumToStr(GE_Coeff[i], 1.9) + "\n");

GE_X   = IIf(BI < GE_BegBar - GE_ExtraB - GE_ExtraF, 0, IIf(BI > GE_EndBar, 0, (GE_XBegin + BI - GE_BegBar + GE_ExtraF) / GE_X_Max));

GE_X2  = GE_X   * GE_X; GE_X3  = GE_X2  * GE_X; GE_X4  = GE_X3  * GE_X; GE_X5  = GE_X4  * GE_X; GE_X6  = GE_X5  * GE_X;
GE_X7  = GE_X6  * GE_X; GE_X8  = GE_X7  * GE_X; GE_X9  = GE_X8  * GE_X; GE_X10 = GE_X9  * GE_X; GE_X11 = GE_X10 * GE_X;
GE_X12 = GE_X11 * GE_X; GE_X13 = GE_X12 * GE_X; GE_X14 = GE_X13 * GE_X; GE_X15 = GE_X14 * GE_X; GE_X16 = GE_X15 * GE_X;

GE_Yn = IIf(BI < GE_BegBar - GE_ExtraB - GE_ExtraF, -1e10, IIf(BI > GE_EndBar, -1e10,
GE_Coeff  +
GE_Coeff  * GE_X   + GE_Coeff  * GE_X2  + GE_Coeff  * GE_X3  + GE_Coeff  * GE_X4  + GE_Coeff  * GE_X5  +
GE_Coeff  * GE_X6  + GE_Coeff  * GE_X7  + GE_Coeff  * GE_X8));

return GE_Yn;
}

function Time_weighting(data, window, BegBar, EndBar) {
BI = BarIndex();
if (window == "Log") alphatime = log(Cum(1));
if (window == "Ramp") alphatime = Cum(1);
if (window == "Exp") alphatime = exp(Cum(0.01));
alphatime = Ref(alphatime, -BegBar);
alphatime = alphatime - alphatime[BegBar];
alphatime = alphatime / alphatime[EndBar];
alphatime = IIf(BI < BegBar, 0, IIf(BI > EndBar, 0, alphatime));
return alphatime*data;
}

function Volume_HighLow_weighting(data, coef_vol, coefdim, period, scale, translation) {
// Coefficient Volume
highvol = V; highvol = HHV(Median(V, 5), 5000);
nbr_vol = coef_vol*highvol/(1 - coef_vol) + 1;
alpha_vol = V/nbr_vol;
alpha_vol = IIf(alpha_vol > 0.99, 0.99, alpha_vol);
alpha_vol = IIf(alpha_vol < 0.01, 0.01, alpha_vol);

// Coefficient High/Low (Adpated from FRAMA from John Ehlers)
coefdim = coefdim * 2 - 0.01;
N = 2*floor(period/2);
delta1 = HHV(Ref(High, -N/2), N/2) - LLV(Ref(Low, -N/2), N/2);
delta2 = HHV(High, N/2) - LLV(Low, N/2);
delta = HHV(High, N) - LLV(Low, N);
Dimen = log((delta1 + delta2)/delta)/log(2);

// Mix Volume and High/Low, and scale them with Inverse Fisher Tansform
data_av = scale*(alpha - translation);
alpha = (exp(2*data_av) - 1)/(exp(2*data_av) + 1);
alpha = Min(Max(alpha,0.01),0.99);
return alpha*data;
}

function Filter_GAP_EOD(data) {
time_frame = Minute() - Ref(Minute(), -1);
time_frame = IIf(time_frame < 0, time_frame + 60, time_frame);
time_frame = time_frame;
day_begin = 152900 + 100*time_frame;
day_end = 215900;
delta = 0;
for (i = 1; i < BarCount; i++) {
if (timequote[i] == Day_begin) delta = data[i] - data[i-1];
data[i] = data[i] - delta;
}
return data;
}

function T3(data, periods, slope) { // High Order Low-Pass filter
e1 = EMA(data, periods);
e2 = EMA(e1, periods);
e3 = EMA(e2, periods);
e4 = EMA(e3, periods);
e5 = EMA(e4, periods);
e6 = EMA(e5, periods);
c1 = -slope*slope*slope;
c2 = 3*slope*slope + 3*slope*slope*slope;
c3 = -6*slope*slope - 3*slope - 3*slope*slope*slope;
c4 = 1 + 3*slope+slope*slope*slope + 3*slope*slope;
result = c1*e6 + c2*e5 + c3*e4 + c4*e3;
return result;
}

function f_centeredT3(data, periods, slope) { // Centered T3 Moving Average
slide = floor(periods/2);
centeredT3 = data;
centeredT3 = Ref(T3(data,periods,slope),slide);
centeredT3 = IIf( IsNan(centeredT3) OR !IsFinite(centeredT3) OR IsNull(centeredT3), data, centeredT3);
return centeredT3;
}

function f_detrend(data, BegBar, EndBar, ExtraF, method_detrend, PF_order) { // Detrend data
BI = BarIndex();
LongBar = EndBar - BegBar + 1 ;
detrended = 0;

if (method_detrend == "Log10") detrended = log10(data);

if (method_detrend == "Derivation") detrended = data - Ref(data, -1);

if (method_detrend == "Linear Regression") {
x = Cum(1); x = Ref(x,-1); x = 0; x = Ref(x, -BegBar);
a = LinRegSlope(data, LongBar); a = a[EndBar];
b = LinRegIntercept(data, LongBar); b = b[EndBar];
y = a*x + b;

global detrending_parameters;
detrending_parameters = y;

y = IIf(BI < BegBar, -1e10, IIf(BI > EndBar + ExtraF, -1e10, y));
Plot( y, "Linear Regression", ParamColor( "Color Linear Regression", colorCycle ), ParamStyle("Style") );
detrended = data - y;
}

if (method_detrend == "PolyFit") {
Yn = PolyFit(data, BegBar, EndBar, PF_order, 0, ExtraF);
Yn = Ref(Yn, -ExtraF);

global detrending_parameters;
detrending_parameters = Yn;

Yn = IIf(BI < BegBar, -1e10, IIf(BI > EndBar + ExtraF, -1e10, Yn));
Plot( Yn, "PolyFit", ParamColor( "Color PolyFit", colorCycle ), ParamStyle("Style") );
detrended = data - Yn;
}

return detrended;
}

function f_retrend(data, Value_BegBar, BegBar, EndBar, method_detrend) {
BI = BarIndex();
retrended = -1e10;

if (method_detrend == "Log10") {
retrended = 10^(data);
}

if (method_detrend == "Derivation") {
retrended[BegBar] = Value_BegBar;
for (i = BegBar + 1; i < EndBar + 1; i++) retrended[i] = data[i] + retrended[i-1];
}

if (method_detrend == "Linear Regression") {
retrended = data + detrending_parameters;
}

if (method_detrend == "PolyFit") {
retrended = data + detrending_parameters;
}

return retrended;
}

function durbinlevison(Autocorr, OrderAR, LongBar, AutoOrder, AutoOrder_Criterion) { // for Yule Walker method only
// Initialization
AR_Coeff = 0;
alpha = noise_variance = Autocorr;
beta = Autocorr;
k = Autocorr / Autocorr;
AR_Coeff = k;

// Initialize recursive criterion for order AR selection
if (AutoOrder == 1) {
FPE = AIC = MDL = CAT = CIC = -1e10;
CAT_factor = 0;
CIC_product = (1 + 1/(LongBar + 1))/(1 - 1/(LongBar + 1));
CIC_add = 1 + 1/(LongBar + 1);
}

// Main iterative loop
for (n = 1; n < OrderAR; n++) {

// Compute last coefficient
for (i = 1; i < n + 1; i++) AR_Coeff_inv[n+1-i] =  AR_Coeff[i];

temp = 0;
for (i = 1; i < n + 1; i++) temp = temp + Autocorr[i] * AR_Coeff_inv[i];
beta[n+1] = Autocorr[n+1] - temp;

alpha[n+1] = alpha[n] * (1 - k[n]*k[n]);
k[n+1] = beta[n+1] / alpha[n+1];
AR_Coeff[n+1] = k[n+1];

// Compute other coefficients by recursion
for (i = 1; i < n + 1; i++) New_AR_Coeff[i] = AR_Coeff[i] - k[n+1] * AR_Coeff_inv[i];
New_AR_Coeff[n+1] =  AR_Coeff[n+1];

// Update
AR_Coeff = New_AR_Coeff;
noise_variance[n+1] = alpha[n+1];

if (AutoOrder == 1) {
i = n + 1;
// Compute criterions for order AR selection
FPE[i] = noise_variance[i]*(LongBar + (i + 1))/(LongBar - (i + 1)); // Final Prediction Error
AIC[i] = log(noise_variance[i])+2*i/LongBar; // Akaine Information Criterion
MDL[i] = log(noise_variance[i]) + i*log(LongBar)/LongBar; // Minimum Description Length
CAT_factor = CAT_factor + (LongBar - i)/(LongBar*noise_variance[i]);
CAT[i] = CAT_factor/LongBar - (LongBar - i)/(LongBar*noise_variance[i]); // Criterion Autoregressive Transfer
CIC_product = CIC_product * (1 + 1/(LongBar + 1 - i))/(1 - 1/(LongBar + 1 - i));
CIC[i] = log(noise_variance[i]) + Max(CIC_product - 1,3*CIC_add); // Combined Information Criterion
}

// End main iterative loop
}

if (AutoOrder == 1) {
// Clean data because of rounding number for very low or very high value, and discard small changes
FPE = IIf(abs(LinRegSlope(FPE, 2)/FPE) < 0.005 OR abs(FPE) > 1e7, -1e10, FPE);
AIC = IIf(log(abs(LinRegSlope(AIC, 2)/AIC)) < 0.005 OR abs(AIC) > 1e7, -1e10, AIC);
MDL = IIf(log(abs(LinRegSlope(MDL, 2)/MDL)) < 0.005 OR abs(MDL) > 1e7, -1e10, MDL);
CAT = IIf(log(abs(LinRegSlope(CAT, 2)/CAT)) < 0.005 OR abs(CAT) > 1e7, -1e10, CAT);
CIC = IIf(abs(LinRegSlope(CIC, 2)/CIC) < 0.005 OR abs(CIC) > 1e7, -1e10, CIC);

// Find lower value
Mixed_Order = 0;
FPE_Order = Mixed_Order = LastValue(Max(ValueWhen(FPE == Lowest(FPE) AND IsFinite(FPE), Cum(1), 1) - 1, 2));
AIC_Order = Mixed_Order = LastValue(Max(ValueWhen(AIC == Lowest(AIC) AND IsFinite(AIC), Cum(1), 1) - 1, 2));
MDL_Order = Mixed_Order = LastValue(Max(ValueWhen(MDL == Lowest(MDL) AND IsFinite(MDL), Cum(1), 1) - 1, 2));
CAT_Order = Mixed_Order = LastValue(Max(ValueWhen(CAT == Lowest(CAT) AND IsFinite(CAT), Cum(1), 1) - 1, 2));
CIC_Order = Mixed_Order = LastValue(Max(ValueWhen(CIC == Lowest(CIC) AND IsFinite(CIC), Cum(1), 1) - 1, 2));

// Mixed array is an averaging from all criterions taking in consideration order choose by criterions is often too low
Mixed_Order_array = Median(Mixed_Order, 5); Mixed_Order = 2*ceil(Mixed_Order_array);

// Print informations about best order for each criterion
printf("\nFPE Order : "+FPE_Order);
printf("\nAIC Order : "+AIC_Order);
printf("\nMDL Order : "+MDL_Order);
printf("\nCAT Order : "+CAT_Order);
printf("\nCIC Order : "+CIC_Order);
printf("\n*** Mixed Order *** : "+Mixed_Order);

// Mark best order in a wich is returned at the end of the function
if (AutoOrder_Criterion == "FPE") AR_Coeff = FPE_Order;
if (AutoOrder_Criterion == "AIC") AR_Coeff = AIC_Order;
if (AutoOrder_Criterion == "MDL") AR_Coeff = MDL_Order;
if (AutoOrder_Criterion == "CAT") AR_Coeff = CAT_Order;
if (AutoOrder_Criterion == "CIC (good)") AR_Coeff = CIC_Order;
if (AutoOrder_Criterion == "Mixed (best)") AR_Coeff = Mixed_Order;
}

AR_Coeff = alpha[OrderAR]; // noise variance estimator
return AR_Coeff;
}

function AR_YuleWalker(data, BegBar, EndBar, OrderAR, Biased, AutoOrder, AutoOrder_Criterion) {
BI = BarIndex();
Data_all = data;
Data = IIf(BI < BegBar, 0, IIf(BI > EndBar, 0, Data));

LongBar = EndBar - BegBar + 1;

// Window Hamming to make it more stable for non-biased estimator
if (Biased == 0) {
pi = atan(1) * 4;
x = Cum(1); x = Ref(x,-1); x = 0;
W = 0.54 - 0.46*cos(2*pi*x/(LongBar-1));
Wi = Ref(W,-BegBar);
data = Wi*data;
}

// Compute Autocorrelation function
for (i = 0; i < OrderAR + 1; i++) {
temp = 0;
for (j = BegBar; j < EndBar + 1 - i; j++) {
temp = temp + data[j]*data[j+i];
}
if (Biased == 1) Autocorr[i]=(1/(LongBar))*temp; // biased estomator (best, model is always stable), small variance but less power
if (Biased == 0) Autocorr[i]=(1/(LongBar-i))*temp; // non-biased estimator (model can be unstable), large variance
}
Autocorr=Autocorr/Autocorr; // normalization (don't change result, to be less exposed to bad rounding)

// Compute AR parameters with Durbin-Levison recursive algorithm for symmetric Toeplitz matrix (Hermitian matrix)
AR_Coeff = durbinlevison(Autocorr, OrderAR, LongBar, AutoOrder, AutoOrder_Criterion);

// End
AR_Coeff = IIf(BI > OrderAR, -1e10, AR_Coeff); // usefull for AR_Pred and AR_Freq
return AR_Coeff; // AR_Coeff is noise variance estimator
}

function AR_Burg(data, BegBar, EndBar, OrderAR, Window, AutoOrder, AutoOrder_Criterion, DoublePrecision) {
BI = BarIndex();
LongBar = EndBar - BegBar + 1; // LongBar = the number of data

// Initialize recursive criterion for order AR selection
if (AutoOrder == 1) {
FPE = AIC = MDL = CAT = CIC = -1e10;
CAT_factor = 0;
CIC_product = (1 + 1/(LongBar + 1))/(1 - 1/(LongBar + 1));
CIC_add = 1 + 1/(LongBar + 1);
}

// Hamming window if selected
if (Window == 1) {
pi = atan(1) * 4;
x = Cum(1); x = Ref(x,-1); x = 0;
W = 0.54 - 0.46*cos(2*pi*x/(LongBar-1));
}

// Create data, forward and backward coefficients arrays
y = Ref(data,BegBar); // y[0:LongBar-1] are the data
y = IIf(BI < LongBar, y, 0); // to be sure not to compute future data
f = b = y; // initialize forward and backward coefficients

// Initialize noise variance estimate with data variance (= data autocorrelation)
noise_variance = Cum(y^2);
noise_variance = noise_variance[LongBar-1]/LongBar;

// Main loop
for (i = 1; i < OrderAR + 1; i++) {

if (Window == 1) Wi = Ref(W,-i); // center Hamming Window
else Wi = 1;

b_shifted = Ref(b,-1);
if (DoublePrecision == 0) { // single precision
k_temp = Sum(Wi*f*b_shifted, LongBar - i)/Sum(Wi*(f^2 + b_shifted^2), LongBar - i);
k[i] = 2*k_temp[LongBar - 1]; // k is reflection coefficient computed from i to LongBar - 1, like Burg tell it (!!! on Numerical Recipes it is computed from 1 to LongBar - i !!!)
}

/* -------------------------------------------------------------------------------------------------------------------------------------------- */
/* --- PART JUST BELOW ARE ONLY FOR EXPERIMENTATION ABOUT ROUNDING PROBLEM, COMPUTATION IS SLOWER AND THIS IS NOT NEEDED FOR SMALL AR ORDER --- */

// Use for FOR loop instead SUM function give different result (neither good, neither bad, but just different) because of rounding single-precision in AFL. De-comment this part to use FOR loop.
/*
Num = 0;Den = 0;
for (j = i; j < LongBar; j++) { // compute from i to LongBar - 1, like Burg tell it (!!! on Numerical Recipes it is computed from 1 to LongBar - i !!!)
Num = Num + Wi[j]*( f[j]*b[j-1] );
Den = Den + Wi[j]*( f[j]^2 + b[j-1]^2 );
}
k[i] = 2*Num/Den; // k is reflection coefficient, a is AR coefficient
*/

/* --- END OF EXPERIMENTAL PART --- */
/* -------------------------------------------------------------------------------------------------------------------------------------------- */

// If you wish Double Precision, VBS script must be used (result are more precise than SUM or FOR loop, but computation is a lot more slower)
if (DoublePrecision == 1) { // double precision
KI_VBS = GetScriptObject();
KI  = KI_VBS.Compute_KI(f, b, Wi, i, LongBar);
k[i] = KI;
}

noise_variance[i] = (1 - k[i]^2)*noise_variance[i-1]; // update noise variance estimator

if (AutoOrder == 1) {
// Compute criterions for order AR selection
FPE[i] = noise_variance[i]*(LongBar + (i + 1))/(LongBar - (i + 1)); // Final Prediction Error
AIC[i] = log(noise_variance[i])+2*i/LongBar; // Akaine Information Criterion
MDL[i] = log(noise_variance[i]) + i*log(LongBar)/LongBar; // Minimum Description Length
CAT_factor = CAT_factor + (LongBar - i)/(LongBar*noise_variance[i]);
CAT[i] = CAT_factor/LongBar - (LongBar - i)/(LongBar*noise_variance[i]); // Criterion Autoregressive Transfer
CIC_product = CIC_product * (1 + 1/(LongBar + 1 - i))/(1 - 1/(LongBar + 1 - i));
CIC[i] = log(noise_variance[i]) + Max(CIC_product - 1,3*CIC_add); // Combined Information Criterion
}

// Update all AR coefficients
for (j = 1; j < i; j++) a[j] = a_old[j] - k[i]*a_old[i-j];
a[i] = k[i];

// Store them for next iteration
a_old = a;

// Update forward and backward reflection coefficients
f_temp = f - k[i]*b_shifted;
b = b_shifted - k[i]*f;
f = f_temp;

// End of the main loop
}

if (AutoOrder == 1) {
// Clean data because of rounding number for very low or very high value, and discard small changes
FPE = IIf(abs(LinRegSlope(FPE, 2)/FPE) < 0.005 OR abs(FPE) > 1e7, -1e10, FPE);
AIC = IIf(log(abs(LinRegSlope(AIC, 2)/AIC)) < 0.005 OR abs(AIC) > 1e7, -1e10, AIC);
MDL = IIf(log(abs(LinRegSlope(MDL, 2)/MDL)) < 0.005 OR abs(MDL) > 1e7, -1e10, MDL);
CAT = IIf(log(abs(LinRegSlope(CAT, 2)/CAT)) < 0.005 OR abs(CAT) > 1e7, -1e10, CAT);
CIC = IIf(abs(LinRegSlope(CIC, 2)/CIC) < 0.005 OR abs(CIC) > 1e7, -1e10, CIC);

// Find lower value
Mixed_Order = 0;
FPE_Order = Mixed_Order = LastValue(Max(ValueWhen(FPE == Lowest(FPE) AND IsFinite(FPE), Cum(1), 1) - 1, 2));
AIC_Order = Mixed_Order = LastValue(Max(ValueWhen(AIC == Lowest(AIC) AND IsFinite(AIC), Cum(1), 1) - 1, 2));
MDL_Order = Mixed_Order = LastValue(Max(ValueWhen(MDL == Lowest(MDL) AND IsFinite(MDL), Cum(1), 1) - 1, 2));
CAT_Order = Mixed_Order = LastValue(Max(ValueWhen(CAT == Lowest(CAT) AND IsFinite(CAT), Cum(1), 1) - 1, 2));
CIC_Order = Mixed_Order = LastValue(Max(ValueWhen(CIC == Lowest(CIC) AND IsFinite(CIC), Cum(1), 1) - 1, 2));

// Mixed array is an averaging from all criterions taking in consideration order choose by criterions is often too low
Mixed_Order_array = Median(Mixed_Order, 5); Mixed_Order = 2*ceil(Mixed_Order_array);

// Print informations about best order for each criterion
printf("\nFPE Order : "+FPE_Order);
printf("\nAIC Order : "+AIC_Order);
printf("\nMDL Order : "+MDL_Order);
printf("\nCAT Order : "+CAT_Order);
printf("\nCIC Order : "+CIC_Order);
printf("\n*** Mixed Order *** : "+Mixed_Order);

// Mark best order in a wich is returned at the end of the function
if (AutoOrder_Criterion == "FPE") a = FPE_Order;
if (AutoOrder_Criterion == "AIC") a = AIC_Order;
if (AutoOrder_Criterion == "MDL") a = MDL_Order;
if (AutoOrder_Criterion == "CAT") a = CAT_Order;
if (AutoOrder_Criterion == "CIC (good)") a = CIC_Order;
if (AutoOrder_Criterion == "Mixed (best)") a = Mixed_Order;
}

// End of the function, return AR coefficients and noise variance estimator for the current selected order in a
a = IIf(BI > OrderAR, -1e10, a); a = noise_variance[i-1]; // usefull for AR_Pred and AR_Freq
return a;
}

function AR_Predict(data, AR_Coeff, EndBar, ExtraF) {
BI = BarIndex();
OrderAR = LastValue(BarCount - BarsSince(AR_Coeff) - 1);

// Print some informations about AR coefficients and noise variance
sum_coeffs = Cum(AR_Coeff) - AR_Coeff; // ARCoeff store noise_variance
printf("\n\nNoise Variance Estimator : " + NumToStr(abs(AR_Coeff), 1.9));
printf("\n\nSum of AR coeffs : " + NumToStr(sum_coeffs[OrderAR], 1.9));
for (i = 1; i < OrderAR + 1; i++)  printf("\nCoeff AR " + NumToStr(i, 1.0) + " = " + NumToStr(AR_Coeff[i], 1.9));

// Pass the data into the AR filter
data_predict = IIf(BI > EndBar, -1e10, Data); // clear data after EndBar
for (i = Endbar + 1; i < EndBar + ExtraF + 1; i++) {
data_predict[i] = 0;
for (j = 1; j < OrderAR + 1; j++)  data_predict[i] = data_predict[i] + AR_Coeff[j] * data_predict[i-j];
}

/* -------------------------------------------------------------------------------------------------------------------------------------------- */
/* --- PART JUST BELOW IS ONLY IF YOU NEED RESIDUALS "u" FROM FILTERING PROCESS OR TO COMPARE IT WITH THE NOISE VARIANCE ESTIMATOR --- */
/*
global BegBar;
data_predict = IIf(BI > EndBar, -1e10, Data);
for (i = BegBar + OrderAR; i < EndBar + 1 + ExtraF; i++) {
data_predict[i] = 0;
for (j = 1; j < OrderAR + 1; j++) {
data_predict[i] = data_predict[i] + AR_Coeff[j] * data_predict[i-j];
}
if (i > EndBar) u[i] = 0; else u[i] = data[i] - data_predict[i];
data_predict[i] = data_predict[i] + u[i];
}
global noise_variance_filterburg;
noise_variance_filterburg = StDev(u, EndBar - (BegBar + OrderAR) + 1); noise_variance_filterburg = noise_variance_filterburg[EndBar]^2;
printf( "\n\nNoise variance from filtering process : " + NumToStr(noise_variance_filterburg, 1.9));
*/
/* --- END OF RESIDUAL COMPUTATION PART --- */
/* -------------------------------------------------------------------------------------------------------------------------------------------- */

// Return data and data predicted from (EndBar + 1) to (EndBar + ExtraF)
return data_predict;
}

function AR_Freq(AR_Coeff, step) {
BI = BarIndex();
pi2 = atan(1) * 8;
OrderAR = LastValue(BarCount - BarsSince(AR_Coeff) - 1);
AR_Coeff=-AR_Coeff;
noise_variance = abs(AR_Coeff);
if (noise_variance == 0) noise_variance = 0.000001; //seven-th digit to one

S = -1e10; // usefull to determine the number of data for AR_Sin function

f = Cum(1); f = Ref(f,-1); f = 0; f = IIf(BI > 0.5/step, -1e10, f);
f = step*f;

// Danielson-Lanczos algorithm to fast compute exponential complex multiplication, using vector formulation AFL for multiple frequency in one-time recursion
W_cos = cos(pi2*f); W_sin = sin(pi2*f); // initialize cos and sinus constant for each frequency f
W_Re = 1; W_Im = 0; // initialize for recursion real and imaginary weight for AR coeffs
Re = 1; Im = 0; // initialize real and imaginary part of the denominator from the spectrum function
for (j = 1; j < OrderAR + 1; j++) {
W_Re_temp = W_Re; // Begin Danielson-Lanczos part
W_Re = W_Re*W_cos - W_Im*W_sin;
W_Im = W_Im*W_cos + W_Re_temp*W_sin; // End Danielson-Lanczos part
Re = Re + AR_Coeff[j]*W_Re; // real part of the denominator
Im = Im + AR_Coeff[j]*W_Im; // imaginary part of the denominator
}

/* Just for experimentation about speed between classic and Danielson-Lanczos algorithm, here is basic computation */
//Re = 1; Im = 0;
//for (j = 1; j < OrderAR + 1; j++) { Re = Re + AR_Coeff[j]*cos(pi2*f*j); Im = Im + AR_Coeff[j]*sin(pi2*f*j); }

S = noise_variance/(Re^2+Im^2); // S the spectrum function
Sdb = 10*log10(abs(s/noise_variance)); // Sdb the spectrum function in dB (power) with noise_variance as reference (so if Sdb is negative, power at this frequency is less than noise power)
return Sdb;
}

function AR_Sin(S, nb_sinus) {
// Initialization
step = 0.5/nb_data;

// Process data to keep only signifiant peaks
slope = LinRegSlope(S,2); // derivation
peaks = Ref(Cross(0, slope), 1); // find maximums of the spectrum
value_peaks = IIf(peaks == 1 AND S > 0, S, 0); // erase maximum wich have less power than the noise

// Find the reduced frequency of the most powerfull cycles
sinus = -1e10; sinus = 0; // to detect error or no cycle
for (i = 0; i < nb_sinus; i++) { // save nb_sinus dominant cycles
sinus[i] = BarCount - LastValue(HighestBars(value_peaks)) - 1; // save higher peak position
value_peaks[sinus[i]] = 0; // erase new detected peak for next iteration
}

sinus = 1/(sinus*step); // convert reduced frequency in periods
sinus = Nz(sinus); // no cycle found

return sinus;
}

/* ---------------------------------------------------------------------- */
/* -------------------------------- DEMO -------------------------------- */
/* ---------------------------------------------------------------------- */

// Choice of the parameters
empty = ParamList("Main parameters", "", 0);
data_source = ParamField("Price field",-1);
ARmodel = ParamList("AR Model", "Burg (best)|Yule-Walker", 0);
length = Param("Length to look back", 200, 1, 5000, 1);
OrderAR  = Param("Order AR", 20, 0, 100, 1);
ExtraF = Param("Number of bar to predict",  50, 0, 200, 1);
AutoOrder = ParamToggle("Auto AR order", "No|Yes", 0);
AutoOrder_Criterion = ParamList("Auto AR criterion", "FPE|AIC|MDL|CAT|CIC (good)|Mixed (best)", 5);
HammingWindow = ParamToggle("Burg: Hamming Window", "No|Yes", 1);
DoublePrecision = ParamToggle("Burg: Double Precision", "No|Yes", 0);
Biased = ParamToggle("Y-W: Biased (Yes: best)", "No|Yes", 1);

empty = ParamList("*** PRE-PROCESSING : ***", "", 0);
empty = ParamList("1- Intraday EOD GAP Filter", "", 0);
FiltrageGAPEOD = ParamToggle("Filtrage GAP EOD", "No|Yes", 0);

empty = ParamList("2- Denoise", "", 0);
periodsT3 = Param("Periods T3", 5, 1, 200, 1);
slopeT3 = Param("Slope T3", 0.7, 0, 3, 0.01);

empty = ParamList("3- Detrend", "", 0);
method_detrend = ParamList("Detrending method", "None|Log10|Derivation|Linear Regression|PolyFit", 3);
PF_order = Param("PolyFit Order", 3, 1, 9, 1);

empty = ParamList("4- Center", "", 0);
PreCenter  = ParamToggle("Center data", "No|Yes", 1);

empty = ParamList("5- Normalize", "", 0);
PreNormalise = ParamToggle("Normalise data", "No|Yes", 1);

empty = ParamList("6- Volume weight on data", "", 0);
PonderVolHL  = ParamToggle("Weighting Volume and H/L(Fractal)", "No|Yes", 0);
coef_vol = Param( "Strength Volume (- +)", 0.5, 0.01, 0.99, 0.01);
coefdim = Param( "Strength Fractal (- +)", 0.5, 0.01, 0.99, 0.01);
period = Param( "Periods Fract (even)", 16, 2, 200, 2);
scale = Param( "Scale Final Weights", 1, 0, 10, 0.01);
translation = Param( "Translate Final Weights", 0.3, 0, 1, 0.01);

empty = ParamList("7- Time weight on data", "", 0);
PonderTime  = ParamToggle("Weighting Time", "No|Yes", 0);
PonderTimeWindow  = ParamList("Weighting Time Window", "Log|Ramp|Exp", 0);

empty = ParamList("*** POST-PROCESSING : ***", "", 0);
method_retrend = ParamList("Retrending method", "Add back trend|Prediction without trend", 0);
empty = ParamList("For Derivation or Log10 :", "add automatically the trend", 0);

// Initialization
SetBarsRequired(20000,20000);
Title = "";
BI = BarIndex();
current_pos = SelectedValue( BI ) - BI[ 0 ];//1500
printf( "Position: " + NumToStr(current_pos) + "\n" );

// Adjust users choice depending number of bars data available or detrending/retrending necessity for derivation method
if ( method_retrend == "Prediction without trend" AND (method_detrend == "Derivation" OR method_detrend == "Log10") ) method_retrend = "Add back trend";

BegBar = current_pos - length + 1;
slide = floor(periodsT3/2);

if (BegBar < 6*periodsT3) Title = Title + "\n\n\n\n\n\n!!! WARNING : YOU HAVE NOT ENOUGH DATA HISTORY FOR CORRECT T3 FILTERING - Reduce \"Periods T3\" parameter OR move Position Bar !!!\n\n\n\n\n\n";

if (BegBar < 1) {
BegBar = 1;
length = current_pos - slide;
Title = Title + "\n\n\n\n\n\n!!! WARNING : YOU HAVE NOT ENOUGH DATA HISTORY - Reduce \"Length to look back\" parameter OR move Position Bar !!!\n\n\n\n\n\n";
}

EndBar = current_pos - slide;
if ( EndBar + ExtraF > BarCount - 1) {
ExtraF = (BarCount - 1) - EndBar;
}

if (AutoOrder == 1) OrderAR = floor(length / 2);
if (OrderAR >= length) {
OrderAR = length - 1;
Title = Title + "\n\n\n\n\n\n!!! WARNING : ORDER AR IS MORE HIGH THAN (LENGTH TO LOOK BACK - 1) - Reduce \"Order AR\" OR increase \"Length to look back\" parameters !!!\n\n\n\n\n\n";
}
data = data_source;
data_filtred = f_centeredT3(data, periodsT3, slopeT3);

/* -------------------- BEGIN PROCESSING -------------------- */

// Filter GAP End of the days - only needed for intraday
if (FiltrageGAPEOD == 1) data = Filter_GAP_EOD(data);

// Denoise data
data = f_centeredT3(data, periodsT3, slopeT3);

// Detrend data
if (method_detrend != "None") data = f_detrend(data, BegBar, EndBar, ExtraF, method_detrend, PF_order);
data = IIf(BI < BegBar, 0, data); // fill with 0 before BegBar to be able to compute native MA Function AFL

// Center data
mean_data_before = MA(data, EndBar - BegBar + 1);
printf( "\nMean data not centered : " + WriteVal(mean_data_before[EndBar]) + "\n" );
if (PreCenter == 1) {
data = data - mean_data_before[EndBar];
mean_data_after = MA(data, EndBar - BegBar + 1);
printf( "\nMean data centered : " + WriteVal(mean_data_after[EndBar]) + "\n\n" );
}

// Normalize data
if (PreNormalise == 1) {
max_data = HHV(data, EndBar - BegBar + 1);
data = data / max_data[EndBar];
}

// Weigthing only for the data which will feed the AR_Burg function (modeling AR), not the AR_Pred function (prediction)
data_weighted = data;
// Volume and High/Low (Fractal) weighting
if (PonderVolHL == 1) data_weighted = Volume_HighLow_weighting(data_weighted, coef_vol, coefdim, period, scale, translation);
// Time weighting
if (PonderTime == 1) data_weighted = Time_weighting(data_weighted, PonderTimeWindow, BegBar, EndBar);

// Compute AR model
if (ARmodel == "Burg (best)") { // Burg model (MEM method) (better than Yule-Walker for sharp spectrum and short data history)
AR_Coeff = AR_Burg(data_weighted, BegBar, EndBar, OrderAR, HammingWindow, AutoOrder, AutoOrder_Criterion, DoublePrecision);
if (AutoOrder == 1) { // In case of AutoOrder, a new computation to find coefficients with the ideal order is needed
OrderAR = AR_Coeff;
AR_Coeff = AR_Burg(data_weighted, BegBar, EndBar, OrderAR, HammingWindow, 0, AutoOrder_Criterion, DoublePrecision); // Compute the AR from AutoOrder
}
}
if (ARmodel == "Yule-Walker") { // Yule-Walker model (Autocorrelation method)
AR_Coeff = AR_YuleWalker(data_weighted, BegBar, EndBar, OrderAR, Biased, AutoOrder, AutoOrder_Criterion);
if (AutoOrder == 1) { // In case of AutoOrder, a new computation to find coefficients with the ideal order is needed
OrderAR = AR_Coeff;
AR_Coeff = AR_YuleWalker(data_weighted, BegBar, EndBar, OrderAR, Biased, 0, AutoOrder_Criterion); // Compute the AR from AutoOrder
}
}
noise_variance = abs(AR_Coeff);

// Prediction based on the computed AR model
data_predict = AR_Predict(data, AR_Coeff, EndBar, ExtraF);

// De-normalize
if (PreNormalise == 1) {
data_predict = data_predict * max_data[EndBar];
noise_variance = noise_variance * max_data[EndBar];
}

// De-center
if (PreCenter == 1) data_predict = data_predict + mean_data_before[EndBar];

// Add the trend back again
if (method_detrend != "None" AND method_retrend == "Add back trend") data_predict = f_retrend(data_predict, data_filtred[EndBar], EndBar, EndBar + ExtraF, method_detrend);

// Translate along value axis filtered EOD GAP data or not retrended data, so no gap occurs
if (FiltrageGAPEOD == 1 OR method_retrend == "Prediction without trend") {
delta = data_predict[EndBar + 1] - data_filtred[EndBar];
data_predict = data_predict - delta;
data_predict = Ref(data_predict, 1); // so there is no discontinuity (but there is one data less for forward prediction)
}

// Plot AR Prediction and Burg predicted channels based on ATR and noise variance (if prediction is bad, channel will be larger)
dataATR = ATR(periodsT3); // no need to detrend or center ATR before prediction (considerer centred without any trend)
printf("\n\nATR AR model for Channel prediction :\n");
AR_Coeff_ATR = AR_Burg(dataATR, BegBar, EndBar, OrderAR, HammingWindow, OrderAR, AutoOrder_Criterion, DoublePrecision);
data_predict_ATR = AR_Predict(dataATR, AR_Coeff_ATR, EndBar, ExtraF);
DeltaBand = 2*(data_predict_ATR + noise_variance);

Data_reconstruct = -1e10;
Data_reconstruct = IIf( BI <= EndBar AND BI >= BegBar, data_filtred, IIf( BI > EndBar, data_predict, Data_reconstruct));
Plot(Data_reconstruct, "AR("+ NumToStr(OrderAR, 1.0) +") Prediction", IIf(BI > EndBar + slide, colorRed, IIf(BI > EndBar AND BI <= EndBar + slide, colorBlue, colorBrightGreen)), styleThick, Null, Null, 0);
Plot(Data_reconstruct + DeltaBand, "AR("+ NumToStr(OrderAR, 1.0) +") Prediction Upper Channel", IIf(BI > EndBar + slide, colorRed, IIf(BI > EndBar AND BI <= EndBar + slide, colorBlue, colorBrightGreen)), styleDashed, Null, Null, 0);
Plot(Data_reconstruct - DeltaBand, "AR("+ NumToStr(OrderAR, 1.0) +") Prediction Lower Channel", IIf(BI > EndBar + slide, colorRed, IIf(BI > EndBar AND BI <= EndBar + slide, colorBlue, colorBrightGreen)), styleDashed, Null, Null, 0);

/* -------------------------------------------------------------------------------------- */
/* Demo for density spectrum power (DSP) analyse based on AR modeling (parametric method) */

// Compute spectrum based directly on AR model. Decomment following line to plot spectrum begining in first Bar (Bar)
Sdb = AR_Freq(AR_Coeff, 0.001); // 0.001 is the step for frequency resolution
//Plot( Sdb, "Spectrum in decibels", ParamColor( "Color Spectrum in decibels", colorCycle ), ParamStyle("Style") );

// Look for the frequency which are the more powerfull (dominant cycles)
sinus = AR_Sin(Sdb, 10); // 10 is the maximum number of frequency to look for
Title = Title + "\nDominant periods: " + WriteIf(sinus >= 0 AND IsFinite(sinus), NumToStr(round(sinus),1.0), "");
for (i = 1; i < 10; i++) Title = Title + WriteIf(sinus[i] > 0 AND IsFinite(sinus[i]), ", "+NumToStr(round(sinus[i]),1.0), "");

/* End DSP demo - For more information look for MEM or MESA                               */
/* -------------------------------------------------------------------------------------- */

/* ---------------------------------------------------------------------- */
/* ------------------------------- END DEMO ----------------------------- */
/* ---------------------------------------------------------------------- */
```
Blogger Comment