2012/12/30

ATC2012が終わりました。

10月から開催されていたATC2012が終了しました。ということで、今回は、その後報告です。

 

結果報告

損益:▲5678.91$  PF0.44 で総合295位

詳しくは以下を参照してください。

kekka 

残念ながら今回も惨敗しました^^;

日本人の方では、MrAwajiさんが48位と健闘され、10名中4名の方がプラスで大会を終了されています。

詳しくは、ここを参照ください。

 

検証

せっかく大会に出たので、大会のトレードとMT5のストラージテスターがどれだけ異なるか確認してみました。

まずは、損益グラフの確認です。

seiseki

やはりストラージテスターとのズレがあるようです。

約定時間を調べてみましたが、分単位まで完全に一致していましたので、約定金額を調べてみると以下の様な結果となりました。(データ区間は、pips単位で各数値は、パーセント形式に加工してあります。)

sas

sa

 

ストラージテスターとATC2012の約定差が発生しているのが図ります。これが、数秒差(分単位以下)の時間のズレか滑ったのが原因かは、解りません。しかし、前回までの大会と比べると飛躍的にストラージテスターとの整合性が取れていると思います。(私のEAがまともにできていたってことかな(・・?)

 

たられば

ココからは、単なる負け惜しみですのでご了承ください。

今回のEAは、単純なBreakOut系のEAがロット調整でどこまでいけるかが課題だったので、その検証です。

まずロット調整なしの損益グラフです。

kote

PF1.09で辛うじてプラスとなっています。

ロット調整を最適化してみると

tes

 

120個のスイッチの内26個を変更した結果です。今回の正解率は、78%で惨敗したわけですw

 

まとめ

たらればの内容は、大会用のお遊びネタですよ。間違ってもリアルのトレードで使おうなどとは考えない方が得策です。どちらかといえば、悪徳EA販売業者が使う手口のような気がします。

本年もいろいろありましたが、何とかすべて乗り越えることが出来ました。ゆっくり新年を迎えたいと思います。

ATCに参加されたみなさんお疲れ様でした。また、来年もお目にかかれることを楽しみにしております。

では^^;

bighope

2012/10/30

続・相関係数

2012/10/30追記
前回の内容にて、相関係数は、計算コストが高すぎてリアルタイムで随時処理していくには不向きなものと記述しましたが、kartzさんからコメントを頂き、それが誤りであったことがわかりましたので記述します。
内容は、私が説明するよりkartzさんのコメントを見て頂いた方が良いと思いますのでそちらをご覧ください。
※正直完全に理解できたか疑問が残るというのが本音だったりしますが。。(汗

コード化

kartzさんのコメントをコード化すると以下の様になりました。
//+------------------------------------------------------------------+
//|                                                   相関係数_i.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_separate_window
#property indicator_buffers 1
#property indicator_color1 Yellow

extern int pd = 100;//期間

double Sx[]; //S(x,a,b) := Σ[i=a~b; x[i]]: xの総和
double Sy[]; //S(y,a,b) := Σ[i=a~b; y[i]]: yの総和
             //M(v,a,b) := S(v,a,b) / (b-a+1): 相加平均
double Fxy[];//F(x,y,a,b) := Σ[i=a~b; (x[i] - M(x,a,b)) * (y[i] - M(y,a,b))]
double Fxx[];//F(x,x,a,b) := Σ[i=a~b; (x[i] - M(x,a,b)) * (x[i] - M(y,a,b))]
double Fyy[];//F(y,y,a,b) := Σ[i=a~b; (y[i] - M(x,a,b)) * (y[i] - M(y,a,b))]
double p[];//相加係数
//+------------------------------------------------------------------+
// | Custom indicator initialization function                        |
//+------------------------------------------------------------------+
int init() {
//---- indicators
   IndicatorBuffers(6);
//---- drawing settings
   SetIndexStyle(0,DRAW_LINE);
   SetIndexStyle(1,DRAW_NONE);
   SetIndexStyle(2,DRAW_NONE);
   SetIndexStyle(3,DRAW_NONE);
   SetIndexStyle(4,DRAW_NONE);
   SetIndexStyle(5,DRAW_NONE);
   SetIndexBuffer(0,p);
   SetIndexBuffer(1,Sx);
   SetIndexBuffer(2,Sy);
   SetIndexBuffer(3,Fxy);
   SetIndexBuffer(4,Fxx);
   SetIndexBuffer(5,Fyy);
//----
   return(0);
  }

//+------------------------------------------------------------------+
//|初期値計算                                                        |
//+------------------------------------------------------------------+
 void Base(int index            //開始位置 
             ,int range           //適用範囲
             ,double &S_x[]     //Sx[];
             ,double &S_y[]     //Sy[];
             ,double &F_xy[]   //Fxy[];
             ,double &F_xx[]   //Fx[];
             ,double &F_yy[]   //Fy[];
                                   ){
   
   int i;
   double x,y,avx,avy,sigxy,sigx,sigy;
   
   //x.yの相加平均
   for(i=range;i>0;i--){
      //xの加算
      avx += High[index+i-1] - Low[index+i-1];
      //yの加算
      avy += Volume[index+i-1];
   }
   avx /=range;
   avy /=range;
            
   for(i=range;i>0;i--){
      x = High[index+i-1] - Low[index+i-1];
      y = Volume[index+i-1];
      sigxy += (x-avx)*(y-avy);
      sigx   += (x-avx)*(x-avx);
      sigy   += (y-avy)*(y-avy);
   }
   S_x[index]  = avx*range;
   S_y[index]  = avy*range;
   F_xy[index] = sigxy;
   F_xx[index] = sigx;
   F_yy[index] = sigy;
}
 
//+------------------------------------------------------------------+
//|相関係数                                                          |
//+------------------------------------------------------------------+
int start(){
   double x,y,xn,yn;
   int    counted_bars=IndicatorCounted();
   int limit=Bars-counted_bars -1;//2012/10/30変更
   //初期設定
   if(counted_bars<=0){
   limit -=pd;//2012/10/30追記
   Base(limit+1,pd,Sx,Sy,Fxy,Fxx,Fyy);
      p[limit+1] = Fxy[limit+1]/MathSqrt(Fxx[limit+1]*Fyy[limit+1]);
   }
   for(int i=limit; i>=0; i--){
      //追加項
      x      = High[i] - Low[i];
      y      = Volume[i];
   //削除項(2012/10/30変更) 
      xn    = High[i+pd] - Low[i+pd];//変更前 xn= High[i+pd+1] - Low[i+pd+1];
      yn    = Volume[i+pd];          //変更前  yn= Volume[i+pd+1];
      //総和の更新
      Sx[i]     = Sx[i+1] + x -xn;
      Sy[i]     = Sy[i+1] + y - yn;
      //F項の更新
       Fxy[i]   = Fxy[i+1] + (Sx[i+1]*Sy[i+1] - Sx[i]*Sy[i])/pd - xn*yn +x*y;
      Fxx[i]   = Fxx[i+1] + (Sx[i+1]*Sx[i+1] - Sx[i]*Sx[i])/pd - xn*xn +x*x;
      Fyy[i]   = Fyy[i+1] + (Sy[i+1]*Sy[i+1] - Sy[i]*Sy[i])/pd - yn*yn + y*y;
      //相関係数
       p[i] = Fxy[i]/MathSqrt(Fxx[i]*Fyy[i]);
   }
return(0);
}
  

表示・確認

従来の計算方法で求めた結果を赤線で表示させ、今回の計算結果を黄線で重ね合わせてみました。
(変更前)
photo


(変更後)※わかりやすくするために従来の計算方法を途中で終了させています。

 

まとめ

この方法なら、随時処理しても十分に使用できます。また、一つ勉強になりました。kartzさんありがとうございました。
//2012/10/30追記
kartzさんから指摘して頂いた、【削除項】の変更の結果、ほぼ一致したラインを作ることが出来ました。

2012/10/27

忘却型相関係数

2012/10/29修正しました。

相関係数は、非常に計算コストの高い関数で、私の様な非力なPCを使っている者にとっては、リアルタイムで随時処理していくには不向きなものだと思っています。ということで、近似関数を用いて相関係数の計算コスト削減を行ってみました。

変換

まずは、相関係数の公式の確認です。対象となるデータをxとyとし期間をnとすると以下のようになります。
soukans
この公式の分母と分子にそれぞれ1/nを乗算すると、以下のようになります。
soukann
ここで、各項を見ると相加平均であるのがわかります。そこで、見慣れた記号SMAに置き換えて見ると
soukanm
となります。そこで、SMAをEMAで置き換えると(近似させると)、
soukane
となり、EMAの公式を当てはめると以下のようになります。
susiki
これで忘却型相関係数の完成です。

確認作業

本来の相関係数と忘却型相関係数をコードにすると以下のようになります。
(※以下のコードは、値幅とティックボリュームの相関関係を求める場合です。)
//+------------------------------------------------------------------+
//|相関係数                                                               |
//+------------------------------------------------------------------+
int start(){
   double x,y,avx,avy,sigxy,sigx,sigy;
   int i,q;
   for(i=limit; i>=0; i--){
      // x:(High[i] - Low[i]) y:Volume[i];
      for(q=pd;q>0;q--){
      //xの加算
         avx += High[i+q-1] - Low[i+q-1];
      //yの加算
         avy += Volume[i+q-1];
      }
      //x.yの相加平均
      avx /=pd;
     avy /=pd;
      
      //変数の初期化
      sigxy =0;
     sigx   =0;
     sigy   =0;
      
      for(q=pd;q>0;q--){
         x = High[i+q-1] - Low[i+q-1];
         y = Volume[i+q-1];
         sigxy += (x-avx)*(y-avy);
         sigx   += (x-avx)*(x-avx);
         sigy   += (y-avy)*(y-avy);
      }
      cn[i] = sigxy/(MathSqrt(sigx)*MathSqrt(sigy));
      
   }
return(0);
}

//+------------------------------------------------------------------+
//|忘却型相関係数                                                          |
//+------------------------------------------------------------------+
int start(){
   double x,y;
   int    counted_bars=IndicatorCounted();
   if(counted_bars>0) counted_bars--;
   int limit=Bars-counted_bars-1;
   
   for(int i=limit; i>=0; i--){
      x      = High[i] - Low[i];
      y      = Volume[i];
      avx[i]   = (1-alfa)*avx[i+1] + alfa*x;
      avy[i]   = (1-alfa)*avy[i+1] + alfa*y;
      sigxy[i] = (1-alfa)*sigxy[i+1] + alfa*(x-avx[i])*(y-avy[i]);  
      sigx[i]  = (1-alfa)*sigx[i+1] + alfa*(x-avx[i])*(x-avx[i]);
      sigy[i]  = (1-alfa)*sigy[i+1] + alfa*(y-avy[i])*(y-avy[i]);
      //忘却型相関係数
      cn[i] = sigxy[i]/(MathSqrt(sigx[i])*MathSqrt(sigy[i]));
   }
return(0);
}
  
チャート上で確認すると,このようになります。
soukan
※チャート上のMAは、それそれの適用期間に合わせたSMAとEMAです。
※上記のサンプルをココに置いておきますので確認してみてください。

まとめ

本来、EMAは、忘却性能と逐次性能が有り、逐次型とか指数平均型とかのネーミングが良いのでしょうが、忘却という響きが好きで、忘却型としました。また、今回は、わかりやすいようにバッファを多用しましたが、1つ前の値しか使用しないので変数(static)で十分対応できます。
ChangeFinder、SDARモデルなどで検索すると、この手法の発展系を見ることが出来ます。

2012/10/13

Notepad++ de MQL5

Notepad++ v6.2が登場したことでハイライト機能が充実しました。と言うことで今回は、MQL5の実装をしてみました。
main


初めに

自動補完用の設定ファイル(MQL5.xml)とハイライト用の設定ファイル(MQL5hl.xml)をココからDLしてくだしい。


自動補完の設定

MQL5.xmlを前回と同様の場所[Notepad++]‐[plugins]‐[\APIs]にコピペします。
※今回の自動補完用ファイルには、関数の引数表示は、設定していません。


ハイライトの設定

Notepad++のバージョンアップに伴い[ユーザー定義ダイアログ]がなくなりました。
設定場所は、【言語】‐【Define-your language…】です。
後は、前回同様です。


ヘルプ(chm)の設定

プラグイン【LanguageHelp】をココからDLします。※plugin manager からは、DLできませんでした。
解凍後、[Notepad++]‐[plugins]にLanguageHelpA.dllとLanguageHelpU.dllをコピペします。
Notepad++を再起動すると、注意表示されますが無視ししてください。
【プラグイン】‐【Language Help】‐【option】と進み【Add】をクリックします。設定は、以下のとおりです。
help
検索したい文字を指定しCTRL+F1で実行できます。


コンパイラーの設定

MetaTrader5は、MQL5.DLLでコンパイルを行いますが、コマンドプロンプトでコンパイル作業が行えるようにMQL5.EXEが用意されています。引数などは、MetaEditor Help → Creating Programs → Compilingに記載されていますので確認してください。あとは、MQL4の時と同じです。
comp


その他

MQL5のエディタで作成されるファイルのフォーマットは、UNICODEです。MQL5を新規に作成する場合は、UT-8のフォーマットにした方が良いとは思いますが、SJISのフォーマットでもコンパイルできました。(※作動確認はしていません)
また、MQL5のMetaEditorで行えるデバッグ機能をNotepad++で行う方法がわかりませんでした。もし可能ならコメント頂けると幸いです。

2012/10/04

Notepad++ de MQL4

秀丸エディタをMQL4エディタ(MetaEditor)にする方法』と言う記事を読んで、面白そうだったのでNotepad++でやってみました。実施範囲は、自動補完・ハイライト・コンパイル作業です。
open


初めに

Notepad++とは
notepad++の導入方法などは、コチラの記事を参考してください。
※日本語入力及び表記に多少問題があるようです。気になるかたは、Notepad++ EUC-JP 対応版の方がよいかもしれません。ただし、今回、確認したのは、Notepad++v6.18(UNICODE)版です。
インストールと初期設定が修了したらココからMQL4_bag.zipをDLし解凍してください。


自動補完の設定

DLしたファイル(mql4.xml)を[Notepad++]‐[plugins]‐[\APIs]にコピペします。
標準インストールした場合はココです。[C:\Program Files\Notepad++\plugins\APIs]
Notepad++を機動させ[設定(T)]‐[環境設定]‐[自動保存/入力補完]と進み以下のように設定します。
para
[入力補完フォーム]
┣入力の毎に補完を行うにチェック
┣関数名補完にチェック
┗入力の直後関数パラメーターヒントを表示にチェック
※お好みによって◎文字目から補完するか設定。
以上で完成です。


テスト:拡張子が.mq4のファイルを作成し、[言語(L)]からMQL4を指定すれば確認できます。
また、引数を伴う関数を指定し(かっこ)をつけるとパラメーターのヒントが表示されます。
sort※関数の引数で()で囲まれた部分がデフォルトです。


ハイライトの設定

Notepad++を機動させ[表示(V)]‐[ユーザー定義ダイアログ]‐[Import...]と進みDLしたファイル(mql4hr.mxl)を指定します。
yuza
yuzaoo
お好きな色を設定して、Notepad++を再起動してください。
※演算子に*と/を追加するとコメントが色抜けします。(原因不明です。)
再起動後MQL4ファイルを開いてもハイライト機能が現れなかったら、[言語(L)]のmql4をクリックしてください。


コンパイラーの設定

Notepad++を機動させ[プラグイン]-[Plugin Manager]-[Show Plugin Manager]と進み[NppExec]をインストールします。
NPPE
再起動後[プラグイン]‐[NppExec]-[Execute(F6)]を表示させ画面の様に入力します。
NPPE3

// save current file
NPP_SAVE 
//compile
C:/Program Files\MetaTrader 4/metalang.exe  "$(FULL_CURRENT_PATH)" 

※C:/Program Files\MetaTrader 4/metalang.exe は、コンパイル先のmetalang.exe のpathです。
入力後[OK]で作動確認をし、問題がなければ名前を付けて[save]します。

次にショートカットの割り当てをします。
[プラグイン]‐[NppExec]‐[Advanced Option]-[Addociated script]にて先ほど保存したファイルを指定し[Add/Modify]をクリックし[OK]で終了します。
NPPE5

[設定(T)]‐[ショートカット管理]‐[Plugin Commands]からファイル名を探しショートカットの設定をします。
NPPE6
NPPE7

これで設定完了です。


まとめ

私自身ほとんどnotepad++を使用していないので何らかのトラブルが生じる可能性がありますのでご了承ください。notepad++のおすすめ機能は、折り畳みと[ctrl]+スクロールでの文字拡大縮小機能です。目が疲れてきたときに助かります。やはり、日本語表記及び記入に問題があるようで、日本語記入後のEnterで補完ログが発生したりします。
また、本来ならHelpの表示もしたかったのですが解りませんでした。

2012/09/29

ATC2012開催前夜

今年で3回目の出場になるATC2012が来週から始まります。
ということで、今回は、3回の大会出場のコンセプト(自分なりの。。)をまとめます。

ATC2010

初出場となったこの年は、MT5をMT4の様に扱う方法にこだわりました。MT5とMT4で大きく異なる点は、ポジション管理です。MT4では、オーダーごとの情報になっているのに対し、MT5では、通貨ごとにポジションがまとめられます。そこで仮想ポジション(オーダー単位)の履歴を把握しMT4の様に作動させることにしました。結果は、思う様に作動しませんでした。
参考:豊嶋先生のブログに仮想ポジションについての記事があります。非常に参考になります。

ATC2011

2回目の出場となったこの年は、複数ペアの取引に挑戦してみました。同一のシステムを複数ペアに取り組むものです。また、メモリの使用量削減にも挑戦しました。が結果は、散々なものとなりました。

ATC2012

3年目となる今回は、資金管理(ロット調整)にこだわって作成してみました。システムの内容は、以下のとおりです。
  • 対象通貨:EURUSD(単一ペア)
  • 対象時間軸:1H
  • システム概要:ブレイクアウト
  • 特徴:必ずBUYorSELLのポジションを持ち続ける。
資金管理なしの場合の成績(0.01lot固定)
  • 期間:2012.01.01~2012.08.01
  • 取引数:326
  • PF:1.11
  • 総損益:631.54

このような、EAに資金調整(ロット調整)を行った結果以下の様になった。
  • 期間:2012.01.01~2012.08.01
  • 取引数:466(ロット調整しかしていないのに取引数が上がっている@@;
  • PF:1.62
  • 総損益:118282.73

※資金調整(ロット調整)の方法は、利益が上がった時間帯には、多くのロットを賭け、利益の上がらない時間帯は、少ないロットで、と言う感じです。本来なら、仮想ポジションを設けて、取引をしないという選択肢もありますが、今回の大会から0.01lotが採用されたのでポジションを持ち続ける事にしました。(単なる手抜き・・

まとめ

今年も何とか、ATC2012に出場することが出来ました。あとは、変なバグや資金が吹っ飛ぶことなく大会終了時まで楽しませてもらえることを祈るばかりです。


おまけ

MQL5リファレンスの機械翻訳(今回は、埋め込み系)を作ってみましたが、9割程度完成したところで心が折れました。。orz 不完全なものですがココに置いておきます。※CHMファイルを見る場合は、ブロックの解除が必要です。

2012/09/08

ヒストリカルデータからエクスポートしたCSVファイルのカスタマイズ

今回は、MT4からヒストリカルデータをエクスポートしたCSVファイルをExcelで加工しやすいようにするVBスクリプトを作ってみました。

【問題点】

MT4からエクスポートしたCSVファイルは、日付のフォーマットがExcelに対応していないため、日付が文字列として認識されます。その為、Excelのありがたい関数郡(例:月別計算or曜日計算or時間の抽出)がそのまま使えません。
histmae

【対応】
そこで、MT4の時間フォーマットをExcelで認識できるフォーマットに変更するスクリプトを作りました。
MT4の時間フォーマット【年.月.日,時間】をExcelフォーマット【年/月/日 時間】に変更します。
スクリプトのDLは、こちら【Time_format_conversion.vbs
Option Explicit
Dim datename,ob,reader,writer,text
Dim re,maches,name,ctext
'ファイル名を取得
If WScript.Arguments.Count = 0 Then
    WScript.Echo "ファイルをドラッグ&ドロップしてください。"
    WScript.Quit
Else
 datename = WScript.Arguments(0)
End If

'オブジェクトの設定
Set ob = CreateObject("Scripting.FileSystemObject")
Set re = New RegExp

'ファイルの読み込み
Set reader = ob.OpenTextFile(datename,1,false)

'新しいファイル名の作成
re.Pattern = "\.csv"
name = re.Replace(datename, "c.csv")

'書き込み用のファイル作成
Set writer = ob.OpenTextFile(name,2,true)

'正規表現のパターン設定
re.Pattern = "(\d{4})\.(\d{2})\.(\d{2})\,"

'書き換え
Do Until reader.AtEndOfStream
    text = reader.ReadLine
    ctext = re.Replace(text,"$1/$2/$3 ")
    writer.WriteLine ctext
Loop
'ファイルを閉じる
reader.Close()
writer.Close()

'完了メッセージ
WScript.Echo name &"を出力しました。"

【使い方】
スクリプトをDL後、ヒストリカルデータからエクスポートしたCSVファイルをスクリプトファイルの上にドラッグ&ドロップしてください。エクスポートしたCSVファイルのあるフォルダにCSVファイル名c.csvというファイルが作成されます。作成されたファイルを開くと以下の様に表示され、Execlで日付が認識されるようになります。
histato
さて、ここまで来たら完成ですが少し遊んでみます。例えば、月別に見た時間別のボラリティの計算をしてみます。
各列にラベルを付けて、月関数【MONTH(シリアル数値)】、時間関数【HOUR(シリアル値)】を追加し、高値と安値の差をHi-Lw(pips)に入れたものです。
hiskakou

これをピボットテーブルで表すと以下の様になり月間別に見た時間ボラリティの平均が出来ます。
pibo
※横軸:月 縦軸:時間です。
グラフ化すると以下の様になります。
burafu
サマータイムの影響が出ているかな(・・?

【まとめ】
VBスクリプトの勉強もかねて作成してみました。



それと。。。今年もATC2012に出場できそうです^^;

2012/05/31

GSLの乱数と確率密度関数を使ってみる。

GSLの乱数と確率密度関数は、条件を付ければラッパーを使わなくてもそのまま使えそうなので試してみました。
※条件とは、GSLには、多くの乱数器がありますが、初期設定の乱数器を使用するという条件です。

ヘッダーの作成

代表的なものだけピックアップし、ヘッダーファイルを作成しました。

//+------------------------------------------------------------------+
//|                                                     gsl_rand.mq4 |
//|                        Copyright 2012, MetaQuotes Software Corp. |
//|                                        http://www.metaquotes.net |
//+------------------------------------------------------------------+
#property copyright "Copyright 2012, MetaQuotes Software Corp."
#property link      "http://www.metaquotes.net"
#property show_inputs

#import "gsl.dll"
//乱数の基本設定
   //環境設定の取得
   int gsl_rng_env_setup ();
   //型T の乱数発生器のインスタンスを生成して、そのハンドルを返す。
   //T:gsl_rng_env_setupで取得したハンドル。
   int gsl_rng_alloc (int T);
   //乱数発生器を初期化
   //r:gsl_rng_allocで取得したハンドル
   //s:初期化種
   void gsl_rng_set (int r, datetime s);
   //乱数発生器のインスタンス r に割り当てられているメモリを解放する。
   void gsl_rng_free (int r);

//以下 r:gsl_rng_allocで取得したハンドル
//乱数

   //0<=x<1の乱数を返す。
   double gsl_rng_uniform (int r);
   //0<x<1の乱数を返す。
   double gsl_rng_uniform_pos (int r);

//分布に従う乱数と確率密度関数
  
   //正規分布乱数
   double gsl_ran_gaussian (int r, double sigma);
   //正規分布の確率密度関数
   double gsl_ran_gaussian_pdf (double x, double sigma);
   
   //平均値mu の指数分布にしたがう乱数を返す。
   double gsl_ran_exponential (int r, double mu);
   //x における平均値mu の指数分布の確率密度関数
   double gsl_ran_exponential_pdf (double x, double mu);
   
   //幅a のラプラス分布にしたがう乱数を返す。
   double gsl_ran_laplace (int r, double a);
   //x における幅a のラプラス分布の確率密度関数
   double gsl_ran_laplace_pdf (double x, double a);
   
   //底がa、指数がb のべき指数分布にしたがう乱数を返す。
   double gsl_ran_exppow (int r, double a, double b);
   //x における底がa、指数がb のべき指数分布の確率密度関数
   double gsl_ran_exppow_pdf (double x, double a, double b);
   
   //係数a のコーシー分布にしたがう乱数を返す。
   double gsl_ran_cauchy (int r, double a);
   //x における係数a のコーシー分布の確率密度関数
   double gsl_ran_cauchy_pdf (double x, double a);
   
   //ガンマ分布にしたがう乱数を返す。
   double gsl_ran_gamma (int r, double a, double b);
   //引数がa とb のときの、x におけるガンマ分布の確率密度関数
   double gsl_ran_gamma_pdf (double x, double a, double b);
   
   //a からb の間に一様に分布する乱数を返す。
   double gsl_ran_flat (int r, double a, double b);
   //a からb の間の一様分布の確率密度関数
   double gsl_ran_flat_pdf (double x, double a, double b);
   
   //対数正規分布にしたがう乱数を返す。
   double gsl_ran_lognormal (int r, double zeta, double sigma);
   //パラメータ値がzeta とsigma で与えられたときの、対数正規分布の確率密度関数
   double gsl_ran_lognormal_pdf (double x, double zeta, double sigma);
   
   //自由度n のカイ二乗分布にしたがう乱数を返す。
   double gsl_ran_chisq (int r, double nu);
   //自由度n のカイ二乗分布の確率密度関数
   double gsl_ran_chisq_pdf (double x, double nu);
   
   //自由度nu1、nu2 のF 分布にしたがう乱数を返す。
   double gsl_ran_fdist (int r, double nu);
   //自由度nu1、nu2 のF 分布の確率密度関数
   double gsl_ran_fdist_pdf (double x, double nu);
   
   //t 分布にしたがう乱数を返す。
   double gsl_ran_tdist (int r, double nu);
   //自由度nu のt 分布の確率密度関数
   double gsl_ran_tdist_pdf (double x, double nu );
   
   //ベータ分布にしたがう乱数を返す。
   double gsl_ran_beta (int r, double a, double b );
   //パラメータ値がa とb で与えられたとき、ベータ分布の確率密度関数
   double gsl_ran_beta_pdf (double x, double a, double b);
   
   //ロジスティック分布にしたがう乱数を返す。
   double gsl_ran_logistic (int r, double a);
   //係数a が与えられたとき<鴻Wスティック分布の確率密度関数
   double gsl_ran_logistic_pdf (double x, double a);
#import


使い方

先回のコードに移植してみました。※変更ヶ所は、赤字部分です。
gsl.dllとgslcblas.dllをterminal.exeと同じ位置に置いています。

//+------------------------------------------------------------------+
//|                                                SMCFilter_gsl.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_separate_window
//#property indicator_chart_window
#include <gsl_rand.mqh>

#property indicator_buffers 2
#property indicator_color1  Yellow
#property indicator_width1  1
#property indicator_color2  Red
#property indicator_width2  1

extern int    Particles_vol = 1000;//粒子の数
extern int    Range         =  100;//試行範囲
extern int    sigmaperiod   =  100;//偏差取得期間
extern double Set           = 3000.0;//初期化範囲(pips)


double Linebuf[];
double Particles[][2];
double ans[];
int r;
//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int init()
  {
//---- indicators
   SetIndexStyle(0, DRAW_LINE);
   SetIndexBuffer(0,Linebuf);
   SetIndexEmptyValue( 0,0);
   
   SetIndexStyle(1, DRAW_LINE);
   SetIndexBuffer(1,ans);
   SetIndexEmptyValue( 1,0);
   
   ArrayResize(Particles,Particles_vol);
   
   //乱数の初期化 
   int T=gsl_rng_env_setup();
   r = gsl_rng_alloc(T);
   gsl_rng_set (r,TimeLocal());
//----
   return(0);
  }
//+------------------------------------------------------------------+
//| Custom indicator deinitialization function                       |
//+------------------------------------------------------------------+
int deinit()
  {
//----
   gsl_rng_free (r);
//----
   return(0);
  }
//+------------------------------------------------------------------+
//|粒子の初期化                                                      |
//+------------------------------------------------------------------+  
void setparticles(double & pobj[][],int range,double set){
   for(int i=0;i<Particles_vol;i++){
      Particles[i][0] = data(range)+(((set*2)/Particles_vol)*i-set)*Point;
      Particles[i][1] = 0; 
    }
 }
//+------------------------------------------------------------------+
//|乱数の成形  0<=random<1                                           |
//+------------------------------------------------------------------+
double random(){
   return(gsl_rng_uniform (r));
}
//+------------------------------------------------------------------+
//|正規分布を計算 mu:平均 sigma:偏差                                |
//+------------------------------------------------------------------+
double normal(double x,double sigma){
   return(gsl_ran_gaussian_pdf(x,sigma));
}
//+------------------------------------------------------------------+
//|正規分布となる乱数の成形 mu:平均 sigma:偏差                       |
//+------------------------------------------------------------------+
double normalrand(double sigma){     
     return(gsl_ran_gaussian (r,sigma));
}
//+------------------------------------------------------------------+
//|リサンプリング q:粒子数                                           |
//+------------------------------------------------------------------+
void resample(int q,double &pobj[][] ){
   
   double wt[];
   double rd;
   int i,j;
   
   //粒子のコピー
   double copyp[][2];
   ArrayResize(copyp,q);
   ArrayCopy(copyp,pobj,0,0,q*2);
   
   //累計重み
   ArrayResize(wt,q);
   wt[0] = pobj[0][2];
   for(i=1;i<q;i++) wt[i]=wt[i-1]+pobj[i][2];
   
   //リサンプリング(ルーレット法)
   for(i=0;i<q;i++){
      rd = random()/Particles_vol;
      for(j=0;j<q;j++){
         if(rd > wt[j]){
            continue;
         }else{
            pobj[i][0] = copyp[j][0];
            pobj[i][1] = 0; 
            break;
         }
      }   
   }
}
//+------------------------------------------------------------------+
//|尤度計算および重み付け                                            |
//+------------------------------------------------------------------+
void wlikelihood(double price,int p,double sigma,double &pobj[][]){
   double alfa[];
   double wsum;
   double g = 1;    //Gの偏微分の値
   int i;
   ArrayResize(alfa,p);
   for(i=0;i<p;i++){
      alfa[i] = normal((price-pobj[i][0]),sigma)*g;//(尤度)
      wsum += alfa[i];
   }
   if(wsum<1)wsum=1;
   //正規化
   for(i=0;i<p;i++)pobj[i][1] = alfa[i]/wsum;
}
//+------------------------------------------------------------------+
//|対象データ                                                        |
//+------------------------------------------------------------------+
double data(int index){
   return(Close[index]);
}
//+------------------------------------------------------------------+
//|偏差の取得                                                        |
//+------------------------------------------------------------------+
double setsigma(int range,int index){
   int      i;
   double sum,av;
   for(i=0;i<range;i++)av += data(index+i);
   av /= range; 
   for(i=0;i<range;i++)sum += MathPow((data(index+i)-av),2);
   return(MathSqrt(sum/range)); 
}
//+------------------------------------------------------------------+
//| Custom indicator iteration function                              |
//+------------------------------------------------------------------+
int start(){
   double v;//乱数
   double sum,sigma,limit;
   int i,j;
   static datetime now = 0;
   int    counted_bars=IndicatorCounted();
   if(counted_bars<1)i=Range; else i=0;
   
   if(now==0)setparticles(Particles,Range,Set);
   while(i>=0){
      //各bar一回のみ実行
      if(now!=Time[i]){
         now=Time[i];
         sigma=setsigma(sigmaperiod,i);
         for(j=0;j<Particles_vol;j++){
            Particles[j][0] = Particles[j][0] + normalrand(sigma);      
         }

         //尤度推定          
         wlikelihood(data(i+1),Particles_vol,sigma,Particles);
         //リサンプリング
         resample(Particles_vol,Particles);   
            
         //加重付平均
         sum=0.0;              
         for(j=0;j<Particles_vol;j++)sum += (Particles[j][0])*Particles[j][1];
      
         Linebuf[i] = sum;
      }
      ans[i] = data(i);//答え
   i--;
   }
   return(0);
  }
//+------------------------------------------------------------------+


まとめ

とりあえず乱数と密度関数で時間を取られることはなくなりそうです。

2012/05/27

モンテカルロフィルタのサンプルを作ってみた。

TwitterのTLに粒子フィルタと言う言葉が出てきたので調べてみました。呼び方がいろいろある様で、モンテカルロフィルタとか、パーティクルフィルタなどと呼ばれているようです。今回は、粒子をベクトル化させない(変化量を与えていない)のでモンテカルロフィルタとしました。

どんなもの?

解りやすくまとめてあるサイトが数多くありますので、参照してください。
手順は、以下の様になります。
  1. モデルの決定
  2. 粒子の初期化(散布)
  3. 予測
  4. 尤度の計算、重み付
  5. リサンプリング
  6. 重み付平均
  7. 3に戻る

作成

今回作成したプログラムの基本設計は、以下のとおりです。また、直前の値からほとんど動かない予測値を想定しています。

【モデル】
システムモデル x[t]=x[-t1]+v (※vは標準正規分布のノイズ)
観測モデル y[t]=x[t]+w (※w=0と設定)
粒子設計 Particles[粒子数][構造] (構造[0:位置情報,1:重み情報])
尤度の計算 標準正規分布を使用
正規分布の偏差 直近データより取得
リターン値 直近の終値を基準にした予測値(黄線)と答えとなる終値(赤線)

【パラメータ】
int    Particles_vol =1000; 粒子数量を設定するパラメータ
int    Range = 100; 試行範囲(粒子数量×試行回数で計算量が増えるので注意)
int    sigmaperiod =100; 標準正規分布のパラメター(偏差)を取得するための計算範囲
double Set = 3000.0; 粒子の初期設定の一様分布の範囲(pipsで記載) -Set<= Particles[][0] <=Set

【コード】

//+------------------------------------------------------------------+
//|                                                    SMCFilter.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_separate_window
//#property indicator_chart_window
#define PI  3.1415926
#property indicator_buffers 2
#property indicator_color1  Yellow
#property indicator_width1  1
#property indicator_color2  Red
#property indicator_width2  1

extern int    Particles_vol = 1000;//粒子の数
extern int    Range         =  100;//試行範囲
extern int    sigmaperiod   =  100;//偏差取得期間
extern double Set           = 3000.0;//初期化範囲(pips)

double Linebuf[];
double ans[];
double Particles[][2]; 
//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int init()
  {
//---- indicators
   SetIndexStyle(0, DRAW_LINE);
   SetIndexBuffer(0,Linebuf);
   SetIndexEmptyValue( 0,0);   
   SetIndexStyle(1, DRAW_LINE);
   SetIndexBuffer(1,ans);
   SetIndexEmptyValue( 1,0);
   //乱数の初期化 
   MathSrand(TimeLocal());
   ArrayResize(Particles,Particles_vol);
//----
   return(0);
  }

//+------------------------------------------------------------------+
//|粒子の初期化                                                      |
//+------------------------------------------------------------------+  
void setparticles(double & pobj[][],int range,double set){
   for(int i=0;i<Particles_vol;i++){
      Particles[i][0] = data(range)+(((set*2)/Particles_vol)*i-set)*Point;
      Particles[i][1] = 0; 
    }
 }
//+------------------------------------------------------------------+
//|乱数の成形  0<=random<1                                           |
//+------------------------------------------------------------------+
double random(){
   double rd = MathRand();
   rd/= 32768;
   return(rd);
}
//+------------------------------------------------------------------+
//|正規分布を計算 mu:平均 sigma:偏差                                |
//+------------------------------------------------------------------+
double normal(double x,double sigma){
   double mu = 0;
   return(MathExp(-MathPow((x-mu),2)/(2*sigma*sigma))/(MathSqrt(2*PI)*sigma));
}
//+------------------------------------------------------------------+
//|正規分布となる乱数の成形 mu:平均 sigma:偏差                       |
//+------------------------------------------------------------------+
double normalrand(double sigma){
     double mu = 0;
     double t,u,r1,r2;
     sigma = MathSqrt(sigma);
     t=MathSqrt(-2.0*MathLog(1-random()));
     u=2*PI*random();
     r1=t*MathCos(u);
     return(r1*sigma+mu);
}
//+------------------------------------------------------------------+
//|リサンプリング q:粒子数                                           |
//+------------------------------------------------------------------+
void resample(int q,double &pobj[][] ){
   double wt[];
   double rd;
   int i,j;
   
   //粒子のコピー
   double copyp[][2];
   ArrayResize(copyp,q);
   ArrayCopy(copyp,pobj,0,0,q*2);
   
   //累計重み
   ArrayResize(wt,q);
   wt[0] = pobj[0][2];
   for(i=1;i<q;i++) wt[i]=wt[i-1]+pobj[i][2];
   
   //リサンプリング(ルーレット法)
   for(i=0;i<q;i++){
      rd = random()/Particles_vol;
      for(j=0;j<q;j++){
         if(rd > wt[j]){
            continue;
         }else{
            pobj[i][0] = copyp[j][0];
            pobj[i][1] = 0; 
            break;
         }
      }   
   }
}
//+------------------------------------------------------------------+
//|尤度計算および重み付け                                            |
//+------------------------------------------------------------------+
void wlikelihood(double price,int p,double sigma,double &pobj[][]){
   double alfa[];
   double wsum;
   double g = 1;    //Gの偏微分の値
   int i;
   ArrayResize(alfa,p);
   for(i=0;i<p;i++){
      alfa[i] = normal((price-pobj[i][0]),sigma)*g;//(尤度)
      wsum += alfa[i];
   }
   if(wsum<1)wsum=1;
   //正規化
   for(i=0;i<p;i++)pobj[i][1] = alfa[i]/wsum;
}
//+------------------------------------------------------------------+
//|対象データ                                                        |
//+------------------------------------------------------------------+
double data(int index){
   return(Close[index]);
}
//+------------------------------------------------------------------+
//|偏差の取得                                                        |
//+------------------------------------------------------------------+
double setsigma(int range,int index){
   int      i;
   double sum,av;
   for(i=0;i<range;i++)av += data(index+i);
   av /= range; 
   for(i=0;i<range;i++)sum += MathPow((data(index+i)-av),2);
   return(MathSqrt(sum/range)); 
}
//+------------------------------------------------------------------+
//| Custom indicator iteration function                              |
//+------------------------------------------------------------------+
int start(){
   double v;//乱数
   double sum,sigma,limit;
   int i,j;
   static datetime now = 0;
   int    counted_bars=IndicatorCounted();
   if(counted_bars<1)i=Range; else i=0;
   
   if(now==0)setparticles(Particles,Range,Set);
   //limit =Set*Point;
   while(i>=0){
      //各bar一回のみ実行
      if(now!=Time[i]){
         now=Time[i];
         sigma=setsigma(sigmaperiod,i);
         for(j=0;j<Particles_vol;j++){
            Particles[j][0] = Particles[j][0] + normalrand(sigma); 
          }
     
         //尤度推定          
         wlikelihood(data(i+1),Particles_vol,sigma,Particles);
         //リサンプリング
         resample(Particles_vol,Particles);   
      
      
         //加重付平均
         sum=0.0;              
         for(j=0;j<Particles_vol;j++)sum += (Particles[j][0])*Particles[j][1];
      
         Linebuf[i] = sum;
      }
      ans[i] = data(i);//答え
   i--;
   }
   return(0);
  }
//+------------------------------------------------------------------+

テスト

稼働させると以下の様になります。
5598
今回のテストでは、直前の終値から次の終値は、ほとんど変化しないという条件で作成したため、このような結果となりました。

まとめ

今回は、線形・ガウス(正規分布)のモデルを使用しましたが、この手法は、非線形・非ガウスへ対応ができる手法の様です。(参考)(・・? 参考までに。。