ライブラリ: 統計関数 - ページ 3

新しいコメント
 

これらの関数の1つについて 質問があります。


定常性の存在を証明するために標準ディッキー・フラー検定が必要です。Dickey Fullerには標準検定と拡張検定があり、ドリフトや定数などを導入した検定もあります。

どの検定が使われているのか、正確にご存知の方がいらっしゃいましたら教えてください。

ありがとうございました。

  
//+------------------------------------------------------------------+
//|統計mqh
//|Copyright 2015, Herajika                  |
//|morinoherajika@gmail.com
//+------------------------------------------------------------------+
#property copyright "ヘラジカ／アダム・S³ッキ"
#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; // オープンバーのディレイ（分
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;           // スムージングの種類 
input ENUM_APPLIED_PRICE   applied_price=PRICE_CLOSE;    // EMA価格の種類

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

int hnd1, hnd2;
double arr1[], arr2[];
input int statSampleSize = 100; //サンプル・サイズを100に制限する。

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)

//+------------------------------------------------------------------+
//| エキスパートティック機能|
//+------------------------------------------------------------------+
void OnTick() {


 //************************************************************************
   //--- バーの本数は十分か？
   if(Bars(_Symbol,_Period) < minBars) {   // バーの合計がminBarsより小さい場合
      Alert(StringFormat("NOT ENOUGH DATA: Less Bars than %d", minBars));
      return;
   }

   //--- 最新2バーの詳細を取得
   if(CopyRates(_Symbol, _Period, 0, 2, mrate) < 2) {
      Alert("Error copying rates/history data - error:", GetLastError(), "!!");
      ResetLastError();
      return;
   }
     
   ZeroMemory(mrequest);      // mrequest構造体の初期化
   
   
   // 指標値を格納する動的配列はarr1[]、arr2[]。
   // すでに上記で宣言済み
   
   
   // 配列のインデックスを、時系列と同じように設定する。
   // indexには最後のバーの値が格納され、indexは最後から1番目などである。
   ArraySetAsSeries(arr1, true);
   ArraySetAsSeries(arr2, true);
   
     
   // インジケーターのハンドルを使って、インジケーターの値をコピーしてみましょう。
   // この目的のために特別に用意されたバッファを配列に変換する。
   CopyBuffer(hnd1,0,0,statSampleSize,arr1);
   CopyBuffer(hnd2,0,0,statSampleSize,arr2);
   
   checkBarOpen();
   
    if (_sqIsBarOpen == true) //CPUとメモリーの使用量が多くなるのを避けるため、バーが開いているときだけ計算するようにする。
   {
          
      // ステップ 1: サンプルデータからトレンドを取り除く (ルーチンに回帰がある) - 配列を返す
      detrend(arr1, detrendres);      
      
      // ステップ2 標本の定常性の検定 - 値を返す
      AugDF = dickeyFuller(detrendres);
      
      //ステップ3 サンプル・コインテグレーションの検定 - 値を返す
      egt =  engleGrangerTest(detrendres, arr1, cointegrationCoeff);
      
      //ステップ4 AR1サンプルの次の値を予測する（Lag = 1の場合） - 値を返す
      AR1forcast = AR1(detrendres);
      
      Print("adf = ", AugDF, "  egt = ", egt , "  ar1 = ",AR1forcast);
  }
   
   
}

//+------------------------------------------------------------------+
//| エキスパート初期化関数|
//+------------------------------------------------------------------+
int OnInit()
{
   
   // すべての配列インジケータ・ハンドルを初期化する
   if(!initIndicators()) return(INIT_FAILED);
   
     
   //コピーされた配列 統計計算用 
   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);
   
}

//+------------------------------------------------------------------+
|エキスパート初期化関数|
//+------------------------------------------------------------------+
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(indicatorHandles[a] == INVALID_HANDLE) { 
         //--- 失敗を伝え、エラーコードを出力する。 
         Print("Failed to create handle of the indicator, error code %d", GetLastError()); 
         //--- インジケータは早期に停止する 
         return(false); 
      }
   }
   
   return(true);
}             


// 注意： 
// 配列のサイズを i に変更した場合、配列には arr[0] から arr[i-1] が含まれる。
// arr[i]は存在しないので、array out of rangeエラーとなる。

//--- 時系列分解（トレンド除去）
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]; // 注意：ここで "array out of range "を得る。

     }
  }
  
 
 //===================================================== 
//--- 回帰線の値を配列で返す
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;
     }
  }
  
 //===================================================== 
//--- 回帰直線の値と回帰直線の係数（y = ax+b の a と 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検定
//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================================= 
//--- βパラメータも返す 
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);
  }
  
  
//===================================================== 
//--- ラグ1の自己回帰モデル；次の期間のリターン予測
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→下、b→上、n→多項式]。
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;
  }
  
//===================================================== 
  
//--- 積分で使用するように編集する
double foo(double x)
  {
   return x;
  }
  
//===================================================== 

//--- クレジット：シット・チー・キーン
double erf(double x)
  {
//--- 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) {
         // 実際のオープンからX分後にオープンするバーを設定する
         processBarOpen = false;

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

      if(processBarOpen) {
         _sqIsBarOpen = true;
         lastBarTime = currentBarTime;      
      } 
   }
}
ヘルプMQL5のstatistics.mqhの更新に合わせて以前のスクリプトを更新しようとしたところ、detrendメソッドで「Array out of range error」に遭遇したようです。何か思い当たることはありますか？
123
新しいコメント