2015年6月3日水曜日

MT4インジケータ 移動平均より滑らかに追従よく! ウェーブレット変換 その2

2015.6.10 このプログラムですが逆変換にバグがあります。正しく動作しません。

Gewinn2のOANDAデモ口座での運用を、みんなのMT4で公開始めました。ブログパーツにはっ付けてあります。
FXCMジャパンでのリアル運用も近いうちに公開予定です。

さて、昨日に引き続き、ウェーブレット変換です。
MT4インジケータ ウェーブレット変換を実装する。その1

引き続き日銀のレポートを参考に実装していきます。
http://www.imes.boj.or.jp/research/papers/japanese/kk23-1-1.pdf

離散ウェーブレット変換ですが、可逆性を保持しています。つまり変換結果から元データが再現可能です。
この特性を生かしたのがウェーブレット逆変換です。

ウェーブレット変換は、ローパスフィルタとハイパスフィルタで信号成分を分離する仕組みです。ハイパスで分離された成分を除去した上で、逆変換を行うと高周波成分が失われた値となります。これを利用してある特定の閾値を持って高周波成分を除去するとノイズフィルタが実装可能です。

ハイパス成分の標本標準偏差未満の値は除去する形で実装したのが、このインジケータとなります。
waveletSNFilter.PNG

EURUSD1時間足のチャートに7EMAとウェーブレットノイズフィルタの結果を表示した結果です。
クローズ値に対する値で、黄色が7EMA、赤がウェーブレットになっています。

水色の枠の中を見ていると、7EMAが暴れているのに対してウェーブレットは暴れていません。その上で全体的に7EMAよりも強く追従していることが判ります。

ウェーブレット変換によるノイズ除去効果が表れています。ちなみに極端に動作したときに逆にウェーブレットの方が暴れている箇所があります。これはウェーブレットフィルタ特性となります。ウェーブレットフィルタを変更することで軽減できます。

ウェーブレット変換の基本的な活用方法を示しましたが、変換結果を持って、類似波形を検索したり、エネルギーを算出して大きな影響を与えている波長を検出したりと活用方法がたくさんあるようです。色々研究してEAに活用したいと思います。

あ、ちなみにこのプログラム結構重たいので、チャートの要素数をそこそこ減らして実行してください。
//------------------------------------------------------------------
// ウェーブレット逆変換ノイズフィルタ

#property copyright "Copyright 2015,  Daisuke"
#property link      "http://mt4program.blogspot.jp/"
#property version   "1.00"
#property strict

#property indicator_chart_window

//バッファーを指定する。
#property indicator_buffers 5

//プロット数を指定する。
#property indicator_plots   1

#property indicator_label1  "reverseWavelet"
#property indicator_type1   DRAW_LINE
#property indicator_color1  clrRed
#property indicator_style1  STYLE_SOLID
#property indicator_width1  1

#property indicator_type2   DRAW_NONE
#property indicator_type3   DRAW_NONE
#property indicator_type4   DRAW_NONE
#property indicator_type5   DRAW_NONE

//入力パラメータ ノイズ除去レベル
input double Noise = 1;

//バッファ
double waveletBuffer[];

//ディテールバッファ
double detail1[];
double detail2[];
double detail3[];
double detail4[];

// 解析レベル
const int level = 2;

//必要データ数
int dataCount = 0;

//------------------------------------------------------------------
//初期化
int OnInit()
{
dataCount = (int)MathPow(2, level);

SetIndexBuffer(0, waveletBuffer);
SetIndexBuffer(1, detail1);
SetIndexBuffer(2, detail2);
SetIndexBuffer(3, detail3);
SetIndexBuffer(4, detail4);

SetIndexDrawBegin(0, dataCount);

string short_name;
short_name = "WV("+ IntegerToString(level)  +")";
IndicatorShortName(short_name);

SetIndexDrawBegin(0, dataCount);

return(INIT_SUCCEEDED);
}

//------------------------------------------------------------------
//計算イベント
int OnCalculate(
const int rates_total,          //各レート要素数
const int prev_calculated,      //計算済み要素数
const datetime &time[],         //要素ごとの時間配列
const double &open[],           //オープン価格配列
const double &high[],           //高値配列
const double &low[],            //安値配列
const double &close[],          //クローズ価格配列
const long &tick_volume[],      //ティック数(要素の更新回数)
const long &volume[],           //実ボリューム(?)
const int &spread[])            //スプレット
{
// フィルタ
double wavebletFilter[];
double scalingFilter[];

ArrayResize(wavebletFilter, 2);
ArrayResize(scalingFilter, 2);

// ハールフィルタ
wavebletFilter[0] =-0.70710678;
wavebletFilter[1] =0.70710678;
scalingFilter[0] = 0.70710678;
scalingFilter[1] = 0.70710678;

//元となる値を計算する。
for( int i = (rates_total - prev_calculated - 1); i>=0 ; i-- )
{
if( i >= (rates_total - 256 - 1) )
{
continue;
}

// 最終的には一個の値に復元される。
double scalingValues[];
ArrayResize(scalingValues, 1);

CalulateWavelet(close, wavebletFilter, scalingFilter, dataCount , level, i, i, scalingValues);
waveletBuffer[i] = scalingValues[0];
}

return(rates_total - 1);
}

//------------------------------------------------------------------
//ウェーブレット値を計算する。
void CalulateWavelet(
const double &values[],          // 元値
const double &wavebletFilter[],  // ウェーブレットフィルタ
const double &scalingFilter[],   // スケーリングフィルタ
int count,                       // ウェーブレット計算サイズ
int targetLevel,                 // ターゲットレベル
int shift,                       // 起点インデックス
int bufferIndex,                 // 現在算出中のバッファインデックス
double &scalingValues[],         // 復元スケール値
int currentLevel = 1             // 現在のレベル
)
{
// count/2 の配列を用意する。
int total = count / 2;

double wavelet[];
double scaling[];
ArrayResize(wavelet, total);
ArrayResize(scaling, total);

int size = ArraySize(wavebletFilter);

for( int i = 0 ; i < total; i++)
{
wavelet[i] = 0;
scaling[i] = 0;
for( int k = 0 ; k < size; k++ )
{
wavelet[i] += wavebletFilter[k] * values[ i * 2 + k + shift];
scaling[i] += scalingFilter[k] * values[ i * 2 + k + shift];
}
}

// グローバルバッファに確保
switch( currentLevel )
{
case 1:
detail1[bufferIndex] = wavelet[0];
break;
case 2:
detail2[bufferIndex] = wavelet[0];
break;
case 3:
detail3[bufferIndex] = wavelet[0];
break;
case 4:
detail4[bufferIndex] = wavelet[0];
break;
}

double reverseScaling[];

if( currentLevel < targetLevel )
{
ArrayResize(reverseScaling, total / 2);
CalulateWavelet(scaling, wavebletFilter, scalingFilter, total, targetLevel, 0, bufferIndex, reverseScaling, ++currentLevel);
}
else
{
ArrayResize(reverseScaling, total);
//スケーリング値はそのまま使用する。
for( int i = 0 ; i < total; i++ )
{
reverseScaling[i] = scaling[i];
}
}

// 標準偏差以下の変動は切り捨てる。
double std = 0;

// 標準偏差用
switch( currentLevel )
{
case 1:
std = StdDev( detail1, 100, bufferIndex);
break;
case 2:
std = StdDev( detail2, 100, bufferIndex);
break;
case 3:
std = StdDev( detail3, 100, bufferIndex);
break;
case 4:
std = StdDev( detail4, 100, bufferIndex);
break;
}
std = std * Noise;

// 逆ウェーブレット変換を行う。
int scalingSize = ArraySize(scalingValues);
int reverseScalingSize = ArraySize(reverseScaling);
for( int i = 0 ; i < scalingSize; i++)
{
scalingValues[i] = 0;
for( int k = 0 ; k < size; k++ )
{
// アップサンプリングの為、偶数インデックスは飛ばす
int index = i + k;
if( index % 2 ==1 ) continue;
index = index / 2;

if( index >= reverseScalingSize ) continue;

scalingValues[i] += scalingFilter[k] * reverseScaling[index];

if( MathAbs( wavelet[index] ) > std )
{
//標準偏差以上のwavelet値だけ反映する。
scalingValues[i] += wavebletFilter[k] * wavelet[index];
}
}
}

}

//------------------------------------------------------------------
//標準偏差を求める。
// return 標準偏差
double StdDev(
const double &values[],      //元となる配列
int count,                   //計算対象数
int shift                    //シフト
)
{
if( (count + shift) > ArraySize(values) ) return 0;

double avg = 0 ;
double mu = 0 ;
for( int i = shift; i < shift + count; i++ )
{
avg += values[i];
}
avg = avg / count;

for( int i = shift; i < shift + count; i++ )
{
mu += MathPow(values[i] - avg, 2);
}

// 標本標準偏差
return MathSqrt(mu / (count - 1));
}


標準偏差は自力で求めています。iStdDevOnArrayは移動平均をとった上で標準偏差を求めるため、単純な移動平均をとるために自力で実装しました。
あと、速度向上の為、グローバルバッファに値を突っ込むなど、少々コードが汚いです・・・。

2015.6.8 プログラムを少し修正しました。ウェーブレット逆変換がうまくいっていない予感がしたので・・・。