2013/11/17

f(x)=a+bX1+cX2+dX3+eX4+fX5+gX6 で重回帰分析

今回は、多項式回帰のコードを流用して複数(6つ以下)の異なる時系列データを用て、元となる時系列データを重回帰分析してみました。

jp225

【仕様説明】

f(x)=a+bX1+cX2+dX3+eX4+fX5+gX6 (a~g:定数 X1~X6:説明変数)
の場合の定数(a~g)を重回帰分析で求め近似モデルf(x)を表示させる。また、上記式に代入した式と近似モデルからの残差二乗和を精度としコメント表示させる。
重回帰分析の範囲は、各説明変数の最古の時系列データの中で最新データのものを基準としパラメータで指定した範囲で分析する。また、欠損データがある場合は、欠損データからもっとも近い過去のデータを使用し補完を行う。

例)
USDJPY(M5)を説明変数(EURUSDとEURJPY)で分析した場合は、以下のようになる。
[USDJPY]f(x)=[97.80020313]+[72.30209778]*EURUSD+[0.73920198]EURJPY 精度=0.00525783
USDJPYTT ※本来USDJPYをEURUSD及びEURJPYで説明する場合は、f(x)=EURJPY/EURUSDですよね^^;
ちなみに、一番初めの画像は、
[NIKKEI225]f(x)=[4382]+[372]*USDJPY+[0.3139]*DJ30+[0.2165]*GOLD+[-501]*JAPAN_BOND+[-223]*5Y_T-NOTES+[225]*10Y_T-NOTES   精度=6954558 (※細部省略) です。

【使い方】

1.ココから重回帰.zipをDLし解凍する。
2.gsl.dllとgslcblas.dllをMT4のterminal.exeと同じ位置に置き、!重回帰v2.mq4をindicatorsフォルダに置く。
3.説明変数となる通貨ペアを表示させ、なるべく多くのデータを取得しておく。
4.!重回帰v2.mq4をチャートにセットし、DLLの使用を許可する。
dll
5.パラメータの説明
N:            重回帰分析の範囲 デフォルトは、1000個のデータを使用
type_price: 0:終値 1:初値 2:高値 3:底値 4:(高値+低値)/2 5:(高値+低値+終値)/3
                 6:(高値+低値+2*終値)/4
zeroset:     一日ごとに重回帰分析を更新する場合は、true そうでない場合は、false
                ※ 更新時刻は、サーバー時間の00:00です。
symb1~6:  説明変数の通貨ペアを記入します。ただし、symb1から記入し使用しない箇所は、空欄(何も記入しない)にしてください。
pea

【まとめ】

MT4の仕様(バグ?)により以下のコードが通りませんでした。(昔からだけど。。。)
//////////////////////////////////////////////////////////////////////////////////
string symb[6] = {USDJPY,EURUSD,EURJPY,AUDJPY,CHFJPY,GBPJPY};
double date[6];
for(int i=0;i<6;i++) date[i] = iMA(symb[i], 0,1,0,MODE_SMA,0,0);
////////////////////////////////////////////////////////////////////////////////////
一度、string型の変数に代入し直してやるかiMA(symb[0],…の様にしてやると通りました。
また、ラッパーを使用せず作成したため行列の範囲を指定して剥きとるgsl関数が使用できませんでした。


【PR】

MT4使いの方でMT5なんて別にいらないんじゃないの?という方へ
MT4の最適化作業にイライラしていませんか?
せっかく高スペックのCPUを積んでいるのにMT4が作動しているのはシングルスレッドのみで、ほかのスレッドが空いていませんか?Core i7 が泣いていますよ!
MT5ならマルチスレッドに対応しているので、最適化作業の時間短縮が測れます。また、複数のPCを繋いだり専用クラウドサービスを使用して最適化作業をすることも可能です。
数日前に『メタトレーダー4&5』という書籍が発売になりました。著者は、豊嶋先生です。これを機にMT5の世界に足を組み入れてみませんか?
いつやるの?。。。いまでしょう! <フルw



【code】


//+------------------------------------------------------------------+
//|                                                    !重回帰v2.mq4 |
//|                        Copyright 2012, MetaQuotes Software Corp. |
//|                                        http://www.metaquotes.net |
//+------------------------------------------------------------------+
#property copyright "Copyright 2012, MetaQuotes Software Corp."
#property link      "http://www.metaquotes.net"
#property indicator_chart_window
#property show_inputs
#property indicator_buffers 1
#property indicator_color1 Snow

//DLL
#import "gsl.dll"
 //ベクトル 
   //長さ n のベクトルを生成し、生成したベクトル構造体へのハンドルを返す。
   int gsl_vector_alloc (int n);
   //v:成形されているベクトル構造体へのハンドル。すでに確保されているベクトル v を解放する。
   void gsl_vector_free (int v);
   //ベクトル v の i 番目の要素を返す。
   double gsl_vector_get (int v, int i);
   //ベクトル v の i 番目の要素に x の値を代入する。
   void gsl_vector_set (int v, int i, double x);
   
 //行列
   //大きさか. n1 行× n2 列の行列を生成し、新しい初期化された行列構造体へのポインタを返す。
   int gsl_matrix_alloc (int n1, int n2);
   //すでに確保されている行列 m を解放する。
   void gsl_matrix_free (int m);
   //行列 m の (i, j) 成分に x の値を代入する。
   void gsl_matrix_set (int m, int i, int j, double x);
   //行列 m の (i, j) 成分を返す。
   double gsl_matrix_get (int m, int i,int j); 
 
 //重回帰
    //この関数は n 変数て. p 個のパラメータを持つモデルで近似するための作業領域を確保する。
    int gsl_multifit_linear_alloc (int n, int p);
    //この関数は作業領域 w に割り当てられたメモリを解放する。
    void gsl_multifit_linear_free (int work); 
    //この関数は観測データ y と予測子変数の行列 X に対する、モデル y = Xc の最良近似パラ
    //メータ c を計算する。モデルパラメータの分散共分散行列 cov は、近似に対するデータの
    //ばらつきから計算される。近似モデルからの残差二乗和 χ2 は chisq に返される。
    int gsl_multifit_linear (int X, int y, int c,int cov, double & chisq[], int work);
#import

//+--------------------------------------------------------------------------------------+
extern int    N          = 1000;          //回帰範囲
extern int    type_price = PRICE_CLOSE;   //取得値
extern bool   zeroset    = false;         //true: 00:00に再計算 false: 再計算しない
extern string symb1      = "EURUSD";      //説明変数1
extern string symb2      = "EURJPY";      //説明変数2
extern string symb3      = "AUDUSD";      //説明変数3
extern string symb4      = "AUDJPY";      //説明変数4
extern string symb5      = "GBPUSD";      //説明変数5
extern string symb6      = "GBPJPY";      //説明変数6
//+--------------------------------------------------------------------------------------+

string symb[6];
int vcont;
double line1[];
//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int init(){
   SetIndexStyle(0,DRAW_LINE);
   SetIndexBuffer(0,line1);
   symb[0]=symb1;
   symb[1]=symb2;
   symb[2]=symb3;
   symb[3]=symb4;
   symb[4]=symb5;
   symb[5]=symb6;
   for(int i=0;i<6;i++)if(symb[i]!="")vcont++; else break;
   return(0);
}
//+------------------------------------------------------------------+
//| Custom indicator deinitialization function                       |
//+------------------------------------------------------------------+
int deinit(){
   Comment("");
   return(0);
}
//+------------------------------------------------------------------+
//| Custom indicator iteration function                              |
//+------------------------------------------------------------------+
int start(){
   string text;
   static double ix[8];
   static int barsht;
   int matricx_xn;
   double line;
   string name,sym;
   int i,q,z,j,s;
   int counted_bars=IndicatorCounted();
   if(counted_bars<1){
     //各データ数で最も少ないbar数の取得
      datetime oldbars;
      for(z=0;z<vcont;z++)oldbars = MathMax(oldbars,iTime(symb[z],0,iBars(symb[z],0)-1));
      int lastbar = iBarShift(NULL,0,oldbars)-1;
      barsht = Bars-lastbar;
     //領域の確保 
      int matricx_xt =  gsl_matrix_alloc(lastbar, vcont+1);
      int matricx_xa=  gsl_matrix_alloc(N, vcont+1);
      double date;
     //全データの取得
      SelectPrices(symb,type_price,0,lastbar,matricx_xt);
     /*確認用
      for(i=0;i<lastbar;i++)
         Print(
            gsl_matrix_get(matricx_xt,i,0),":"
            ,gsl_matrix_get(matricx_xt,i,1),":"
            ,gsl_matrix_get(matricx_xt,i,2),":"
            ,gsl_matrix_get(matricx_xt,i,3),":"
            ,gsl_matrix_get(matricx_xt,i,4),":"
            ,gsl_matrix_get(matricx_xt,i,5),":"
            ,gsl_matrix_get(matricx_xt,i,6)
         );
       */
     //重回帰分析用のデータの取得
      for(i=(lastbar-N);i>0;i--){
         if(i== lastbar-N ||(Time[i]%86400 ==0 && zeroset)){
            for(j=0;j<N;j++){
               for(q=0;q<vcont+1;q++){
                  date = gsl_matrix_get(matricx_xt,i+j,q);
                  gsl_matrix_set (matricx_xa,j,q,date);
               }
            }
            //重回帰分析 
            Multiple_regression(i,N,type_price,matricx_xa,ix);
         }
         line=0.0;
         for(q=0;q<=vcont;q++)line += ix[q]*gsl_matrix_get(matricx_xt,i,q);
         line1[i] = line;
      }
      gsl_matrix_free(matricx_xt);
      gsl_matrix_free(matricx_xa);
   }else{
      i=Bars-counted_bars;
      while(i>=0){
         if((ix[0]==0)||(Time[i]%86400 ==0 && zeroset)){
            matricx_xn = gsl_matrix_alloc(N, vcont+1);
            SelectPrices(symb,type_price,i,N,matricx_xn);
            Multiple_regression(i,N,type_price,matricx_xn,ix);
            gsl_matrix_free(matricx_xn);
         }     
         
         line1[i] = ix[0];
         for(s=0;s<vcont;s++){
            sym = symb[s];
            line1[i] += ix[s+1]*iMA(sym,0,1,0,MODE_SMA,type_price,i);
         }
         i--;
      }
   }
   
   text = "["+Symbol()+"]f(x) =["+ix[0] +"]";
   for(s=0;s<vcont;s++) text = text + " + ["+ix[s+1]+"]*"+symb[s];
   text= text +"      精度="+ix[7];
    Comment(text);
   return(0);
  }
//+------------------------------------------------------------------+
//|データの取得
//+------------------------------------------------------------------+
void SelectPrices(string& sym[],int pricetype,int index,int range,int matricx_handle){
   datetime sytime;
   int s,i,q,t;
   double price;
   //f(x) = a + b*x1+c*x2+.....g*x6
   //a
   for(int z=0;z<range;z++)gsl_matrix_set (matricx_handle, z, 0,1.0);
   //x1...x6
   for(s=0;s<vcont;s++){
      if(sym[s]=="")break;
      SelectPrice_symb(sym[s],s,pricetype,index,range,matricx_handle);
   }
}
//+------------------------------------------------------------------+
//|データの取得(通貨ペア)
//+------------------------------------------------------------------+
void SelectPrice_symb(string symbl,int symb_i,int pricetype,int index,int range,int matricx_h){
   datetime sytime;
   int q =0,t=0;
   for(int i=0;i<range;i++){
         sytime = iTime(symbl,0,q+i+index);
         for(t=0;t<100;t++)   {
            if(Time[i+index]>=sytime)break; else q++;
            sytime = iTime(symbl,0,q+i+index);
         }
         gsl_matrix_set (matricx_h, i, symb_i+1,iMA(symbl,0,1,0,MODE_SMA,pricetype,q+i+index));
         if(Time[i+index] > sytime)q--;
   }
}
//+------------------------------------------------------------------+
//|重回帰分析
//+------------------------------------------------------------------+
void  Multiple_regression(int index,int range,int type,int matricx_xd,double & ans[]){
   int matrix_cov,vector_y,vector_c,multifit_work;
   int i;
   double chisq[1];
   //f(x) = a + b*x1+c*x2+.....g*x6
   //作業領域の確保
    matrix_cov     =  gsl_matrix_alloc(vcont+1,vcont+1);
    vector_y       =  gsl_vector_alloc(range);
    vector_c       =  gsl_vector_alloc(vcont+1);
    multifit_work  =  gsl_multifit_linear_alloc (range,vcont+1);

  //y
   for(i=0;i<range;i++)gsl_vector_set (vector_y , i,iMA(NULL,0,1,0,MODE_SMA,type,i+index));
   
  //処理
   multifit_work = gsl_multifit_linear_alloc (range,vcont+1);
   gsl_multifit_linear(matricx_xd,vector_y,vector_c,matrix_cov,chisq,multifit_work);
   gsl_multifit_linear_free (multifit_work);
   
  //結果の取得 
  //ans[8] = {a,b,c,d....chisq}
   for(i=0;i<vcont+1;i++) ans[i] = gsl_vector_get(vector_c,i);
   ans[7] = chisq[0];
   
  //メモリ解放
   gsl_matrix_free(matrix_cov);
   gsl_vector_free(vector_y);
   gsl_vector_free(vector_c);
}