2011/12/06

ウェーブレットを用いたパターン認識■■■

※2011/12/9追記(コードのバグ修正) しました。
ウェーブレットの勉強をしている際に、ふと思いついたことを形にしてみました。


パターン認識とは

今回ここで挙げるパターン認識とは、過去のチャートの形から、現在のチャートの形に近いものを選択し、その後の動きを元に、近未来のチャートの動きを予測しようというものです。

ここで問題になったのが、2点。
1.どのような形を用いるか。
※今回は、ウェーブレット変換した形を採用しました。
2.どのように予測するか。
※最も近い形3ヶ所をピックアップし、その2時間後までの高値、安値を導く。その線を表示させて、検証する。(実は、まだ検討中。)


ウェーブレット

今回使用したウェーブレットは、前に作ったHarrのウェーブレットをウェーブレット変換のみ使用しています。
前作のウェーブレット変換を完了させるとg[]配列に以下の様に変換された数値が、格納される。
(L=2^N:N=5の場合)
g[0]=レベルNのスケーリング係数
g[1]=レベルNのウェーブレット係数
g[2~3]=レベルN-1のウェーブレット係数
g[4~7]=レベルN-2のウェーブレット係数
g[8~15]=レベルN-3のウェーブレット係数
g[16~31]=レベルN-4のウェーブレット係数
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
ここで、最下層レベルN-4(高周波領域)をフィルタリングによりのウェーブレット係数を0にすると、格納されているデータは、半分の量になります。画像で見てみると、上がg[]配列のフィルタ前、下がフィルタ後です。さらに、レベルN-2までフィルタリングするとデータ量は、さらに半減します。
N0
このように圧縮(加工)されたデータをパターン認識の形として採用しました。
ちなみに、レベルN-4でフィルタリングしたデータを再変換をすると、以下の様な画像となります。
N2
圧縮(加工)したデータを用いることで、パターン認識する測点が減り、ノイズの削除効果もあるためウェーブレット変換を採用しました。


その他もろもろ

照合するパターンの抽出間隔を1日(1440分=86400秒)にしてあります。これは、抽出作業の際に特異な期間(短期間)に抽出箇所が集中するのを防ぐ役割と、前回のブログの話で出てきた時間単位の平均値(期待値)の差を低減させるためです。また、夏時間の影響をカバーするために夏時間の変換コードも追加しておきました。
※5分足様に作成されています。
採用した値は、すべてOpenを採用しています。これは、演算負荷を低減させるために採用した値です。お好みで調整ください。


コード

//+------------------------------------------------------------------+
//|                                                         DWTP.mq4 |
//|                                        Copyright ゥ 2011, bighope |
//|                       http://expertadviser-bighope.blogspot.com/ |
//+------------------------------------------------------------------+
#property copyright "Copyright ゥ 2011, bighope"
#property link      "http://expertadviser-bighope.blogspot.com/"

#property indicator_chart_window
#property indicator_buffers 6
#property indicator_color1 Red
#property indicator_color2 Red
#property indicator_color3 Orange
#property indicator_color4 Orange
#property indicator_color5 Yellow
#property indicator_color6 Yellow

#define sqr 0.70710678

//---- input parameters
extern int       N      =   7;//対象データ数(2のN乗)
extern int       HBs    =   5;//高周波カット位置
extern int       LBs    =   7;//低周波カット位置
extern int       Shift  =   0;//データのシフト
extern int       loop   = 500;//検索回数
int    DWTPeriod ;
//---- buffers
double Hi0Band[];
double Lw0Band[];
double Hi1Band[];
double Lw1Band[];
double Hi2Band[];
double Lw2Band[];
double g[];
double gs[];
double now[];
int CusPeriod;
//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int init()
  {
//---- indicators
   SetIndexStyle(0,DRAW_LINE);
   SetIndexStyle(1,DRAW_LINE);
   SetIndexStyle(2,DRAW_LINE);
   SetIndexStyle(3,DRAW_LINE);
   SetIndexStyle(4,DRAW_LINE);
   SetIndexStyle(5,DRAW_LINE);
   SetIndexShift(0,24);
   SetIndexShift(1,24);
   SetIndexShift(2,24);
   SetIndexShift(3,24);
   SetIndexShift(4,24);
   SetIndexShift(5,24);
   SetIndexBuffer(0,Hi0Band);
   SetIndexBuffer(1,Lw0Band);
   SetIndexBuffer(2,Hi1Band);
   SetIndexBuffer(3,Lw1Band);
   SetIndexBuffer(4,Hi2Band);
   SetIndexBuffer(5,Lw2Band);
   
   CusPeriod = MathPow(2,N-HBs+1);
   DWTPeriod = MathPow(2,N);
   ArrayResize(g,DWTPeriod);
   ArrayResize(gs,DWTPeriod);
   ArrayResize(now,CusPeriod);
//----
   return(0);
  }
//+------------------------------------------------------------------+
//| Custom indicator deinitialization function                       |
//+------------------------------------------------------------------+
int deinit()
  {
//----
   Comment("");
//----
   return(0);
  }
  
//+------------------------------------------------------------------+
//|  function                                                        |
//+------------------------------------------------------------------+ 
bool str(int index,double B,int &No[],double &BL[]){
   int Hindex;
   Hindex=ArrayMaximum(BL);
   
   if(BL[Hindex]>B){
      BL[Hindex] = B;
      No[Hindex]  = index;
      return(true);
   }
   return(false);
}

//+------------------------------------------------------------------+
//| ウェーブレット変換                                               |
//+------------------------------------------------------------------+ 
 void DWT(double &input[],double &out[],int Ns,int DWTf,int Hibs,int Lwbs){
   int z,k;
    double difference;
   for(z=1 ; z<=Ns ;z++){
      DWTf =DWTf/2;
      for(k = 0; k < DWTf; k++){
         out[k]         = (input[k*2]+input[k*2+1])*sqr;
         difference = (input[k*2]-input[k*2+1])*sqr;
         if((Lwbs>=z)&&(z>=Hibs)){out[DWTf+k] = difference;}else{out[DWTf+k] = 0;}//フィルタリング作業
       }
      if(Ns>z)ArrayCopy(input,out);  
   }
 }

//+------------------------------------------------------------------+
//|夏時間関数                                                        |
//+------------------------------------------------------------------+ 
bool  DST(datetime Now)//夏時間にtrue
{
 bool Check = false;
 int Start_Month    =              3;//開始月
 int Start_Week     =              2;//第○週
 int Start_DayofWeek=              0;//0:日曜日
 int End_Month      =             11;//終了月
 int End_Week       =              1;//第○週
 int End_DayofWeek  =              0;//0:日曜日
 int Year_Now       =  TimeYear(Now);
 datetime DST_Start,DST_End ;
 string MakeTime;
 int Day_S,Day_E;
 
 if(Year_Now < 2000)return(0);
 switch(Year_Now)
 {
  case 2000:
            DST_Start = D'02.04.2000';
            DST_End   = D'29.10.2000';
            break;
  case 2001:
            DST_Start = D'01.04.2001';
            DST_End   = D'28.10.2001';
            break;
  case 2002:
            DST_Start = D'07.04.2002';
            DST_End   = D'27.10.2002';
            break;
  case 2003:
            DST_Start = D'06.04.2003';
            DST_End   = D'26.10.2003';
            break;
  case 2004:
            DST_Start = D'04.04.2004';
            DST_End   = D'31.10.2004';
            break;
  case 2005:
            DST_Start = D'03.04.2005';
            DST_End   = D'30.10.2005';
            break;
  case 2006:
            DST_Start = D'02.04.2006';
            DST_End   = D'29.10.2006';
            break;
  default:
            MakeTime =Year_Now + "."+Start_Month + ".01"; 
            DST_Start = StrToTime(MakeTime) - 86400.0;
            Day_S = Start_Week*7 - TimeDayOfWeek(DST_Start) + Start_DayofWeek;
            MakeTime = Year_Now + "."+ Start_Month  + "." + Day_S;
            DST_Start = StrToTime(MakeTime);
            
            MakeTime =Year_Now + "."+End_Month + ".01"; 
            DST_End = StrToTime(MakeTime) - 86400.0;
            Day_E = End_Week*7 - TimeDayOfWeek(DST_End) + End_DayofWeek;
            MakeTime = Year_Now + "."+ End_Month  + "." + Day_E;
            DST_End = StrToTime(MakeTime);
            break;
 }
 
 if(DST_Start < Now &&  Now < DST_End)Check = true;
 
 return(Check);
 }
 
//+------------------------------------------------------------------+
//| Custom indicator iteration function                              |
//+------------------------------------------------------------------+

int start()
  {

   int counted_bars=IndicatorCounted();
   if(counted_bars>0&&Volume[0]>1)return(0);
   int y=0,k=0,index=0;
   datetime day=0;
   int    BarNo[3]={0,0,0};
   double     L[3]={1000.0,1000.0,1000.0};
   double    R=1000.0;
   double    Hi[3]={0,0,0};
   double    Lw[3]={0,0,0};
   
   for(int i=0;i<loop;i++){
      ArrayInitialize(g,0);
      ArrayInitialize(gs,0); 
//※修正箇所
      day= Time[Shift]-84400*i;
   //day=Time[Shift]-1440*i;訂正前
      if(DST(day))day-=3600;
  //if(DST(day))day--;訂正前
//※訂正箇所ここまで 
      if(i!=0)index = iBarShift(NULL,0,day,true);
      if(index<0){ Print("Orver data",i); break;}
      for(k=0;k<DWTPeriod;k++) g[k]=(Open[k+Shift+index]-Open[Shift+DWTPeriod+index])*MathPow(10,Digits-2);
      DWT(g,gs,N,DWTPeriod,HBs,LBs);
      
      if(i==0){
         ArrayCopy(now,gs,0,0,CusPeriod);
      }else{
             R=0;
             for(y=0;y<CusPeriod;y++) R += MathPow(now[y]-gs[y],2);
             str(Shift+index,R,BarNo,L);
           }
   }
    R=1000.0;
    for(k=0;k<3;k++){
      Hi[k] = High[iHighest(NULL,0,MODE_HIGH,24,BarNo[k]-24)]-Open[BarNo[k]];
      Lw[k] = Low[iLowest(NULL,0,MODE_LOW,24,BarNo[k]-24)]-Open[BarNo[k]];
      if(R>L[k]){R=L[k];index=k;}
    }

    for(k=0;k<24;k++){
      switch(index){
         case 0:
               Hi0Band[Shift+k] = Hi[0]+Open[Shift];
               Lw0Band[Shift+k] = Lw[0]+Open[Shift];
               Hi1Band[Shift+k] = Hi[1]+Open[Shift];
               Lw1Band[Shift+k] = Lw[1]+Open[Shift];
               Hi2Band[Shift+k] = Hi[2]+Open[Shift];
               Lw2Band[Shift+k] = Lw[2]+Open[Shift];
               break;
         case 1:
               Hi0Band[Shift+k] = Hi[1]+Open[Shift];
               Lw0Band[Shift+k] = Lw[1]+Open[Shift];
               Hi1Band[Shift+k] = Hi[0]+Open[Shift];
               Lw1Band[Shift+k] = Lw[0]+Open[Shift];
               Hi2Band[Shift+k] = Hi[2]+Open[Shift];
               Lw2Band[Shift+k] = Lw[2]+Open[Shift];
               break;
         case 2:
               Hi0Band[Shift+k] = Hi[2]+Open[Shift];
               Lw0Band[Shift+k] = Lw[2]+Open[Shift];
               Hi1Band[Shift+k] = Hi[0]+Open[Shift];
               Lw1Band[Shift+k] = Lw[0]+Open[Shift];
               Hi2Band[Shift+k] = Hi[1]+Open[Shift];
               Lw2Band[Shift+k] = Lw[1]+Open[Shift];
               break;
         default:
               Print("Err Buffer select");
               break;
      }
 
    }
   Comment(" R = ",L[index]," : ",L[0]," : ",L[1]," : ",L[2], "   No= ",BarNo[index]," : ",BarNo[0]," : ",BarNo[1]," : ",BarNo[2]);
  
//----
   return(0);
  }
//+------------------------------------------------------------------+

使いかた

パラメータは、ほぼ前回のウェーブレットの時に用いた物と同じですので省略して、追加したパラメータの説明をすると、【loop】は、検索するデータ量を設定するものです。デフォルトの500とは、過去500日分のデータを照合しなさいと言うことになります。(※設定したデータ量がない場合は、Printでお知らせし、そこでbreakします。)

表示内容の説明
インジケータを表示すると以下の様になると思います。
y0
各水平ランは、色ごとのペア(3種類のペア)で2時間後までの高値と安値を表示しています。ちなみに赤色のラインが最も近い形の予測値となります。
コメント欄(画像左上)の【R】は、検索された物がどれだけ近い形を示しているかを表す数値です。(小さい方が良い)【R=最も近かった数値:R1:R2:R3(3種類のR値)】を表示しています。【No】は、検索された値の日数を表します。【No=最も近かった数値の日数(●日前):R1の日数:R2の日数:R3の日数】
しばらく作動させると履歴が残りますのでどんな動きだったか検証することもできますが、パラメータ【Shift】の数値を変更しても確認できます。

まとめ

いろいろテストしてみましたが、完全に予測するというものではありません。ただし、Rの数値は小さく各3組の線がほぼ同じ状態になった時に、予測できているかな?と言う程度です。では!