Libraries: Statistical Functions - page 3

[Victor Adamuz]

Hi, I have a question about one of these functions.

I need the standard Dickey Fuller test to establish the existence of stationarity. I have seen that it is in this set, but I don't really know which test it uses, since Dickey Fuller has both the standard test and the augmented test and within them we have several, those that introduce drift, constant, etc.

I have seen the code but the formulas don't fit, I have also checked the results with other programs like R and I don't get the same, so I would like to know if anyone knows exactly which test is being used, so I know if I am using the same test in all programs, or see where the differences are etc.

Thank you very much

[Jacob James]

//+------------------------------------------------------------------+
//|                                                   statistics.mqh |
//|                        Copyright 2015, Herajika                  |
//|                                         [e-mail] |
//+------------------------------------------------------------------+
#property copyright "Herajika/Adam S³ucki"
#property link      "[e-mail]"

//+------------------------------------------------------------------+
#include <Expert/Expert.mqh>
#include <Trade/SymbolInfo.mqh>
#include <Trade/OrderInfo.mqh>
#include <Trade/HistoryOrderInfo.mqh>
#include <Trade/Trade.mqh>
#include <Trade/PositionInfo.mqh>
#include <Trade/DealInfo.mqh>
#include <Trade/TerminalInfo.mqh>  
#include <Trade/AccountInfo.mqh>  
#include <Object.mqh>
#include <MovingAverages.mqh>
#include <Arrays\ArrayObj.mqh>
#include <Math\Stat\Math.mqh>
#include <Math\Alglib\statistics.mqh>
//--------------------------------------------------------------

MqlTradeRequest mrequest;  // Used for sending our trade requests
MqlTradeResult mresult;    // Used to get our trade results
MqlRates mrate[];          // Used to store the prices, volumes and spread of each bar
CTerminalInfo terminalInfo;
                       
datetime startTime = TimeCurrent();
datetime lastCheckedTradeTime = startTime;
                       
bool sqDisplayInfoPanel = MQLInfoInteger(MQL_TESTER) != 0 && MQLInfoInteger(MQL_OPTIMIZATION) != 0;
int sqLabelCorner = 1;
int sqOffsetHorizontal = 5;
int sqOffsetVertical = 20;
color sqLabelColor = clrWhite;

int magicNumber;
int minBars = 30;
int sqMaxSlippage;
int sqVerboseMode;
int sqMaxRetries = 5;

double gPointCoef = 0;                  

int deviationPoints = 10;
int sqTicket = 1;
datetime lastBarTime;
bool mcond[100];

double initialBalance = 0;

string valueIdentificationSymbol = "";

int indicatorHandles[];
//+------------------------------------------------------------------+


bool _sqIsBarOpen;
input int OpenBarDelay = 0; // open bar delay in minutes
input int ma_period1  = 55;
input int ma_shift1  = 0 ;
input int ma_period2  = 144;
input int ma_shift2  = 0 ;
input ENUM_MA_METHOD       movavg_method=MODE_EMA;           // type of smoothing 
input ENUM_APPLIED_PRICE   applied_price=PRICE_CLOSE;    // type of EMA price

input string statistics = "-----------statistics analysis  -----------";

int hnd1, hnd2;
double arr1[], arr2[];
input int statSampleSize = 100; //Limit sample size to only 100 values

double detrendres[] , olsres[];
double aLinCoeff = 0.0, bLinCoeff = 0.0, cointegrationCoeff = 0.0;
double AugDF,egt,AR1forcast;

#define EMA_1 0                   //iMA(NULL,0,ma_period1,ma_shift1,ma_method,PRICE_CLOSE)
#define EMA_2 1                   //iMA(NULL,0,ma_period2,ma_shift2,ma_method,PRICE_CLOSE)

//+------------------------------------------------------------------+
//| Expert tick function                                             |
//+------------------------------------------------------------------+
void OnTick() {


 //************************************************************************
   //--- Do we have enough bars to work with?
   if(Bars(_Symbol,_Period) < minBars) {   // if total bars is less than minBars
      Alert(StringFormat("NOT ENOUGH DATA: Less Bars than %d", minBars));
      return;
   }

   //--- Get the details of the latest 2 bars
   if(CopyRates(_Symbol, _Period, 0, 2, mrate) < 2) {
      Alert("Error copying rates/history data - error:", GetLastError(), "!!");
      ResetLastError();
      return;
   }
     
   ZeroMemory(mrequest);      // Initialization of mrequest structure
   
   
   // Dynamic arrays to store indicators values  are arr1[], arr2[]
   // Already declared above
   
   
   // Setting the indexing in arrays the same as in timeseries, i.e. array element with zero
   // index will store the values of the last bar, with ith index - the last but one, etc.
   ArraySetAsSeries(arr1, true);
   ArraySetAsSeries(arr2, true);
   
     
   // Using indicators handles, let's copy the values of indicator
   // buffers to arrays, specially prepared for this purpose
   CopyBuffer(hnd1,0,0,statSampleSize,arr1);
   CopyBuffer(hnd2,0,0,statSampleSize,arr2);
   
   checkBarOpen();
   
    if (_sqIsBarOpen == true) //To avoid high CPU & memory usage let's do calculation only at Bar Open
   {
          
      //Step 1:  Remove Trend from the Sample Data (has regression in-routine) - returns an Array
      detrend(arr1, detrendres);      
      
      //Step 2 Test for Sample Stationarity   -  Returns a Value
      AugDF = dickeyFuller(detrendres);
      
      //Step 3 Test for Sample Cointegration  -  Returns a Value
      egt =  engleGrangerTest(detrendres, arr1, cointegrationCoeff);
      
      //Step 4 Forcast next value of sample AR1 if Lag = 1  - Returns a Value
      AR1forcast = AR1(detrendres);
      
      Print("adf = ", AugDF, "  egt = ", egt , "  ar1 = ",AR1forcast);
  }
   
   
}

//+------------------------------------------------------------------+
//| Expert initialization function                                   |
//+------------------------------------------------------------------+
int OnInit()
{
   
   //Initialize All Array Indicator handles
   if(!initIndicators()) return(INIT_FAILED);
   
     
   //copied Arrays For mathematical statistics calculations 
   hnd1 = indicatorHandles[EMA_1];     //iMA(NULL,0,ma_period1,ma_shift1,ma_method,applied_price)
   hnd2 = indicatorHandles[EMA_2];     //iMA(NULL,0,ma_period2,ma_shift2,ma_method,applied_price)
   
   return(INIT_SUCCEEDED);
   
}

//+------------------------------------------------------------------+
//| Expert deinitialization function                                 |
//+------------------------------------------------------------------+
void OnDeinit(const int reason)
  {
    writeReportFile(); 
  }


bool initIndicators(){
   
   ArrayResize(indicatorHandles, ArraySize(indicatorHandles) + 1, 10);
   indicatorHandles[EMA_1] = iMA(NULL,0,ma_period1,ma_shift1,movavg_method,applied_price) ;
   
   ArrayResize(indicatorHandles, ArraySize(indicatorHandles) + 1, 10);
   indicatorHandles[EMA_2] = iMA(NULL,0,ma_period2,ma_shift2,movavg_method,applied_price) ;
   
         
   for(int a=0; a<ArraySize(indicatorHandles); a++){
      //--- if the handle is not created 
      if(indicatorHandles[a] == INVALID_HANDLE) { 
         //--- tell about the failure and output the error code 
         Print("Failed to create handle of the indicator, error code %d", GetLastError()); 
         //--- the indicator is stopped early 
         return(false); 
      }
   }
   
   return(true);
}             


// NOTE: 
// If you resize an array to i the array contains arr[0] to arr[i-1]
// arr[i] does not exist and so you will get the array out of range error.

//--- time series decomposition (remove trend)
void detrend(double &arr[],double &res[])
  {
   double regRes[];
   double iterator[];
   int length=ArraySize(arr);
   ArrayResize(regRes,length);
   ArrayResize(iterator,length);
//---
   for(int i=0; i<length; i++)
     {
      iterator[i]=i+1;
     }
//---
   regression(iterator,arr,regRes);
//---

   for(int i=0; i<length;i++)
     {
      Print("   arr = ", arr[i], "   regRes  = ", regRes[i], "    SumXY =", arr[i] - regRes[i] );
      res[i] = arr[i] - regRes[i]; // NOTE: getting "array out of range"  HERE

     }
  }
  
 
 //=====================================================  
//--- return array with regression line values
void regression(double &tsarr1[],double &tsarr2[],double &res[])
  {
//---
   double a = 0;
   double b = 0;
//---
   int length=ArraySize(tsarr1);
//---
   ArrayResize(res,length);
//---
   double meanX = MathMean(tsarr1);
   double meanY = MathMean(tsarr2);
//---
   double sumXY=0;
   double sqSumX=0;
//---
   for(int i=0; i<length; i++)
     {
      Print("   tsarr1 = ", tsarr1[i], "   tsarr2  = ", tsarr2[i], "    SumXY =", tsarr1[i]*tsarr2[i] );
      sumXY+=tsarr1[i]*tsarr2[i];
      sqSumX+=MathPow(tsarr1[i],2);
     }
//---
   a = (sumXY - length*meanX*meanY)/(sqSumX - length*MathPow(meanX,2));
   b = meanY - a*meanX;
//---
   for(int i=0; i<length; i++)
     {
      res[i]=a*tsarr1[i]+b;
     }
  }
  
 //=====================================================   
//--- return array with regression line values, plus regression line coefficients (a and b from y = ax+b)
void regression(double &tsarr1[],double &tsarr2[],double &res[],double &aCoeff,double &bCoeff)
  {
//---
   double a = 0;
   double b = 0;
//---
   int length=ArraySize(arr1);
//---
   ArrayResize(res,length);
//---
   double meanX = MathMean(tsarr1);
   double meanY = MathMean(tsarr2);
//---
   double sumXY=0;
   double sqSumX=0;
//---
   for(int i=0; i<length; i++)
     {
      sumXY+=tsarr1[i]*tsarr2[i];
      sqSumX+=MathPow(tsarr1[i],2);
     }
//---
   a = (sumXY - length*meanX*meanY)/(sqSumX - length*MathPow(meanX,2));
   b = meanY - a*meanX;
//---
   for(int i=0; i<length; i++)
     {
      res[i]=a*tsarr1[i]+b;
     }
   aCoeff = a;
   bCoeff = b;
  }
  
//======================Test for STATIONARITY===============================   
//--- Augumented Dickey Fuller test for stationarity
//bool dickeyFuller(double &arr[])
double dickeyFuller(double &arr[])
  {
// n=25     50     100    250    500    >500
// {-2.62, -2.60, -2.58, -2.57, -2.57, -2.57};
   double cVal;
   bool result;
   int n=ArraySize(arr);
   double tValue;
   double corrCoeff;
   double copyArr[];
   double difference[];
   ArrayResize(difference,n-1);
//---
   for(int i=0; i<n-1; i++)
     {
      difference[i]=arr[i+1]-arr[i];
     }
//---
   ArrayCopy(copyArr,arr,0,0,n-1);
   //corrCoeff=myStats.PearsonCorr2(copyArr,difference,statSampleSize);
   corrCoeff=myStats.SpearmanCorr2(copyArr,difference,statSampleSize);
   tValue=corrCoeff*MathSqrt((n-2)/(1-MathPow(corrCoeff,2)));
//---
   if(n<25)
     {
      cVal=-2.62;
        }else{
      if(n>=25 && n<50)
        {
         cVal=-2.60;
           }else{
         if(n>=50 && n<100)
           {
            cVal=-2.58;
              }else{
            cVal=-2.57;
           }
        }
     }
//Print(tValue);
   result=tValue>cVal;
   return(result);
   return (tValue);
  }
 
 
//====================Test for COINTEGRATION=================================   
//--- return also beta parameter 
bool engleGrangerTest(double &tsarr1[],double &tsarr2[],double &cointCoeff)
  {
   bool result;
   int length=ArraySize(tsarr1);
   double regressionRes[];
   double residuals[];
   double copyArr1[],copyArr2[];
   double a;
   double b;
//---
   ArrayResize(regressionRes,length);
   ArrayResize(residuals,length);
   ArrayResize(copyArr1,length);
   ArrayResize(copyArr2,length);
//---
   ArrayCopy(copyArr1,tsarr1,0,0);
   ArrayCopy(copyArr2,tsarr2,0,0);
//---
   for(int i=0; i<length;i++)
     {
      copyArr1[i] = MathLog(copyArr1[i]);
      copyArr2[i] = MathLog(copyArr2[i]);
     }
//---
   regression(copyArr1,copyArr2,regressionRes,a,b);
   cointCoeff=a;
//---
   for(int i=0; i<length; i++)
     {
      residuals[i]=copyArr2[i]-regressionRes[i];
     }
//---
   result=dickeyFuller(residuals);
//---
   return(result);
  }
  
  
//=====================================================  
//--- Autoregressive model with lag 1; return forecast for next period
double AR1(double &arr[])
  {
   double arrY[],arrX[],regRes[];
   double a,b,fcst;
   int n=ArraySize(arr);
   ArrayResize(arrY,n-1);
   ArrayResize(arrX,n-1);
//---
   ArrayCopy(arrY,arr,0,1,n-1);
   ArrayCopy(arrX,arr,0,0,n-1);
//---
   regression(arrX,arrY,regRes,a,b);
   fcst=arr[n-1]*a+b;
//---
   return(fcst);
  }
  
//===================================================== 
  
//--- [a -> down, b -> up, n -> polynominal degree]
double signedIntegral(double a,double b,int n)
  {
//---
   int mult;
   double h=(b-a)/n;
   double integral=foo(a);
//---
   for(int i=1; i<n; i++)
     {
      if(i%2==0) mult=4; else mult=2;
      integral+=mult*(foo(a+i*h));
     }
//---
   integral += foo(a+n*h);
   integral *= h/3;
//---
   return integral;
  }
  
//===================================================== 
  
//--- edit to use with integral
double foo(double x)
  {
   return x;
  }
  
//===================================================== 

//--- credits to Sitt Chee Keen
double erf(double x)
  {
//--- A&S formula 7.1.26
   double a1 = 0.254829592;
   double a2 = -0.284496736;
   double a3 = 1.421413741;
   double a4 = -1.453152027;
   double a5 = 1.061405429;
   double p=0.3275911;
   x=MathAbs(x);
   double t=1/(1+p*x);
//---
   return 1 - ((((((a5 * t + a4) * t) + a3) * t + a2) * t) + a1) * t * MathExp(-1 * x * x);
  }
  
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
double normDistZ(double z)
  {
   double sign=1;
   if(z<0) sign=-1;
   return 0.5 * (1.0 + sign * erf(MathAbs(z)/MathSqrt(2)));
  }
//+------------------------------------------------------------------+

void checkBarOpen(){
   datetime currentBarTime = mrate[0].time;   
   _sqIsBarOpen = false;
   
   if(lastBarTime == 0){
      _sqIsBarOpen = true;          
      lastBarTime = currentBarTime;
   }
   else if(currentBarTime != lastBarTime){
      bool processBarOpen = true;

      if(OpenBarDelay > 0) {
         // set bar to open after X minutes from real open
         processBarOpen = false;

         int diffInSeconds = (int) (TimeCurrent() - currentBarTime);
         if(diffInSeconds >= OpenBarDelay * 60) {
            processBarOpen = true;
         }
      }

      if(processBarOpen) {
         _sqIsBarOpen = true;
         lastBarTime = currentBarTime;      
      } 
   }
}

HELP: I seem to run into " Array out of range error" in detrend method after trying to update the former script to comply with updated MQL5 statistics.mqh . Any idea ? what am i missing ? any guidance well appreciated