Bibliotecas: Statistical Functions - página 3

Novo comentário
 

Olá, tenho uma dúvida sobre uma dessas funções.


Preciso do teste padrão de Dickey Fuller para estabelecer a existência de estacionariedade. Vi que ele está nesse conjunto, mas não sei realmente qual teste ele usa, já que o Dickey Fuller tem tanto o teste padrão quanto o teste aumentado e dentro deles temos vários, os que introduzem drift, constante, etc.

Vi o código mas as fórmulas não se encaixam, também verifiquei os resultados com outros programas como o R e não obtive o mesmo, então gostaria de saber se alguém sabe exatamente qual teste está sendo utilizado, para eu saber se estou utilizando o mesmo teste em todos os programas, ou ver onde estão as diferenças etc.

Muito obrigado

  
//+------------------------------------------------------------------+
//|estatísticas.mqh
//|Direitos autorais 2015, Herajika |
//|morinoherajika@gmail.com |
//+------------------------------------------------------------------+
#property copyright "Herajika/Adam S³ucki"
#property link      "morinoherajika@gmail.com"

//+------------------------------------------------------------------+
#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; // atraso da barra aberta em minutos
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;           // tipo de suavização 
input ENUM_APPLIED_PRICE   applied_price=PRICE_CLOSE;    // tipo de preço da EMA

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

int hnd1, hnd2;
double arr1[], arr2[];
input int statSampleSize = 100; //Limitar o tamanho da amostra a apenas 100 valores

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)

//+------------------------------------------------------------------+
//| Função de tique de especialista|
//+------------------------------------------------------------------+
void OnTick() {


 //************************************************************************
   //--- Temos barras suficientes para trabalhar?
   if(Bars(_Symbol,_Period) < minBars) {   // se o total de barras for menor que minBars
      Alert(StringFormat("NOT ENOUGH DATA: Less Bars than %d", minBars));
      return;
   }

   //--- Obtenha os detalhes das últimas 2 barras
   if(CopyRates(_Symbol, _Period, 0, 2, mrate) < 2) {
      Alert("Error copying rates/history data - error:", GetLastError(), "!!");
      ResetLastError();
      return;
   }
     
   ZeroMemory(mrequest);      // Inicialização da estrutura mrequest
   
   
   // As matrizes dinâmicas para armazenar os valores dos indicadores são arr1[], arr2[]
   // Já declarado acima
   
   
   // Definir a indexação em matrizes da mesma forma que em séries temporais, ou seja, elemento de matriz com zero
   // índice armazenará os valores da última barra, com o i-ésimo índice - o penúltimo, etc.
   ArraySetAsSeries(arr1, true);
   ArraySetAsSeries(arr2, true);
   
     
   // Usando os identificadores de indicadores, vamos copiar os valores do indicador
   // buffers para matrizes, especialmente preparados para essa finalidade
   CopyBuffer(hnd1,0,0,statSampleSize,arr1);
   CopyBuffer(hnd2,0,0,statSampleSize,arr2);
   
   checkBarOpen();
   
    if (_sqIsBarOpen == true) /Para evitar o alto uso de CPU e memória, vamos fazer o cálculo somente no Bar Open
   {
          
      //Etapa 1: remover a tendência dos dados de amostra (tem regressão na rotina) - retorna uma matriz
      detrend(arr1, detrendres);      
      
      //Etapa 2 Teste de estacionariedade da amostra - Retorna um valor
      AugDF = dickeyFuller(detrendres);
      
      //Etapa 3 Teste de cointegração de amostra - Retorna um valor
      egt =  engleGrangerTest(detrendres, arr1, cointegrationCoeff);
      
      //Etapa 4 Previsão do próximo valor da amostra AR1 se Lag = 1 - Retorna um valor
      AR1forcast = AR1(detrendres);
      
      Print("adf = ", AugDF, "  egt = ", egt , "  ar1 = ",AR1forcast);
  }
   
   
}

//+------------------------------------------------------------------+
//| Função de inicialização especializada|
//+------------------------------------------------------------------+
int OnInit()
{
   
   //Inicializar todos os identificadores do indicador de matriz
   if(!initIndicators()) return(INIT_FAILED);
   
     
   //cópias de matrizes para cálculos de estatísticas matemáticas 
   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);
   
}

//+------------------------------------------------------------------+
//| Função de desinicialização de especialista|
//+------------------------------------------------------------------+
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++){
      //--- se o identificador não for criado 
      if(indicatorHandles[a] == INVALID_HANDLE) { 
         //--- informar sobre a falha e emitir o código de erro 
         Print("Failed to create handle of the indicator, error code %d", GetLastError()); 
         //--- o indicador é interrompido antecipadamente 
         return(false); 
      }
   }
   
   return(true);
}             


// OBSERVAÇÃO: 
// Se você redimensionar uma matriz para i, a matriz conterá arr[0] a arr[i-1]
// arr[i] não existe e, portanto, você receberá o erro de matriz fora do intervalo.

//--- decomposição de séries temporais (remover tendência)
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]; // OBSERVAÇÃO: obtendo "matriz fora do intervalo" AQUI

     }
  }
  
 
 //===================================================== 
//--- retorna uma matriz com valores de linha de regressão
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;
     }
  }
  
 //===================================================== 
//--- retorna a matriz com os valores da linha de regressão, mais os coeficientes da linha de regressão (a e b de 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=============================== 
//--- Teste de Dickey Fuller aumentado para estacionariedade
//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);
  }
  
  
//===================================================== 
//--- Modelo autorregressivo com defasagem 1; previsão de retorno para o próximo período
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 -> para baixo, b -> para cima, n -> grau polinominal]
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;
  }
  
//===================================================== 
  
//--- editar para usar com integral
double foo(double x)
  {
   return x;
  }
  
//===================================================== 

//--- créditos a Sitt Chee Keen
double erf(double x)
  {
//--- Fórmula A&S 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) {
         // definir barra para abrir após X minutos da abertura real
         processBarOpen = false;

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

      if(processBarOpen) {
         _sqIsBarOpen = true;
         lastBarTime = currentBarTime;      
      } 
   }
}
AJUDA: Parece que estou encontrando o erro " Array out of range error" no método detrend depois de tentar atualizar o script anterior para que ele esteja em conformidade com o MQL5 statistics.mqh atualizado. Alguma ideia? O que estou perdendo? Qualquer orientação será muito bem-vinda
123
Novo comentário