GARCHによるボラティリティをチャートに表示して眺めてみたかったので作りました。
推定するのはモデルにARMAや平均値を含めないピュアGARCH(p,q)です。
RにはfGarchとtseriesのインストールが必要です。
mt4Rは開発者7bit氏のオリジナルではなく、miccllyさんによるMT4 build600以降対応版で。
build670で作成しましたが722でも動作することを確認しています。
MT4からRを呼んで関数を記述したソースファイルを読み込ませる仕様です。
バックグラウンドのRの動作状況はDebugViewで確認することができます。
先にインジケーターのパラメーターの説明をしてしまいます。
Sizeで指定したサンプルサイズでのモデル推定を、PastBars分遡った足からローリングします。
デフォルトではどちらも1000なので、この場合は階差分1本を含めると最低2000本の足が必要です。
描画されるのは推定されたモデルを用いて予測した1期先の値（条件付き標準偏差）、
いわゆるひとつのボラティリティです。
pとqは次数で、もちろん増やすこともできますがあまり良い結果は出ません。
最適化に失敗する確率が高くなり、係数の信頼性も低下します。
ティック毎に計算するのはしんどいので足が増えたときにのみ新しく計算をします。
それでもできるだけ最新の価格を使いたかったので、現在の足の始値をデフォルトとしましたが、
PrevCloseをTRUEにすると前足の終値を使います。
これについてはちょっとした発見というか謎がありましたので後述。
fGarchをFALSEにするとtseriesで計算します。
この違いについても後ほど。
Testモードはこのインジケーターを参照するEAのバックテストをMT4で行う時のためのものです。
通常はFALSEのままで。
設定にもよりますがモデル推定には時間が掛かるため、
Testモードによる同期実行ではRが計算している間MT4はほぼフリーズ状態で、
ティック受信すらままならなくなるので注意が必要です。
※実際にテストするにはこれだけではたぶん不十分で、
EAの方にも変数や配列にきちんと値が格納されたか確認する記述が必要でしょう。
RcheckはRの生存確認用で、死亡が確認された場合、Rの再起動を試みます。
DebugViewで監視しているときはティック毎に返事が返って来るのがうるさいし、
R側のソースにバグがあるような場合、無限に再起動を繰り返してMT4が固まったりするので、
私はFALSEがデフォルトです。
Infoは最新の計算結果に関する情報を返すか、返すならどこに返すか、です。
計算に掛かった時間（正確には計算開始から終了後の次のティックまでの時間）や、
モデル推定に失敗した数などが出力されます。
計算終了前に足が増えた場合の遅延分の補正はできていると思いますが、
タイムフレームの切り替え、スリープや回線障害からの復帰時等は処理があやしいというか、
そのような場面でのMT4の挙動に合わせ切れていないので、インジケーターを再起動した方がよいでしょう。
インジケーターの基本的な構造は上述の7bit氏のものや、
faiさんのBTLMを参考にさせて頂いておりますが、
本プログラムに関する実装ミスやヒヨコードな部分の責任はもちろんすべて私にあります。
インジケーターの起動時、プロパティ画面でDLLの使用にチェックが必須です。
// MT4R_GARCH_Predict
#property indicator_separate_window
#property indicator_minimum 0
#property indicator_buffers 1
#property indicator_color1 Red
#property indicator_width1 2
#include <mt4R.mqh>
extern int Size = 1000;
extern int p = 1;
extern int q = 1;
extern int PastBars = 1000;
extern bool PrevClose = FALSE;
extern bool fGarch = TRUE;
extern bool Test = FALSE;
extern bool Rcheck = FALSE;
extern int Info = 1; // 0:off, 1:print, 2:comment
int R;
double SIGMAt[];
int init(){
SetIndexLabel(0, "SIGMAt");
SetIndexBuffer(0, SIGMAt);
SetIndexStyle(0, DRAW_LINE);
// 使用するRterm.exeの場所を指定
R = RInit("C:/R/R-3.1.1/bin/x64/Rterm.exe --no-save", 2);
// ソースファイルの名前と場所
RExecute(R, "source('c:/R/src_GARCHpred.R')");
RAssignInteger(R, "size", Size);
RAssignInteger(R, "p", p);
RAssignInteger(R, "q", q);
if(Info == 1) Print("Init: R(" + R + "), Bars=", Bars);
return(0);
}
int deinit(){
RDeinit(R);
Comment("");
return(0);
}
int start(){
int SampleSize;
int Loop;
int i;
static int LastCalcBar = 0;
static datetime CalcStartTime = 0;
double NewPrice[];
double Result[];
// Rの非同期計算が終了しているかどうかの確認
if(RIsBusy(R)) return(0);
// Rの生存確認と再起動の試行
if(Rcheck){
bool Ronline = false;
Ronline = RGetBool(R,"online");
if(!Ronline){
Alert(Symbol() + Period() + " R(" + R + ") Offline, Trying Restart Indicator.");
deinit();
init();
CalcStartTime = 0;
LastCalcBar = 0;
return(0);
}
}
// 計算結果の描画
if(CalcStartTime != 0){
bool Rstatus = RGetBool(R, "status");
if(Rstatus){
// Rの計算結果の取得
int ResultSize = RGetInteger(R, "length(result)");
ArrayResize(Result, ResultSize);
RGetVector(R, "rev(result)", Result, ResultSize);
int Fault = RGetInteger(R, "fault");
// 計算実行中に増えた足の確認
int Delay = Bars - LastCalcBar;
// 計算結果に関する情報出力
string InfoStr = StringConcatenate("Delay=" + Delay +
", Calculated: " + TimeToStr(TimeLocal(), TIME_MINUTES) +
" in " + (TimeLocal() - CalcStartTime) +
" sec, CalcPeriod: " + TimeToStr(Time[Size - 1 + Delay], TIME_DATE|TIME_MINUTES) +
" - " + TimeToStr(Time[Delay]) +
", Length=" + ArraySize(Result) +
", Fault=" + Fault);
if(Info == 1) Print(InfoStr);
else if(Info == 2) Comment(InfoStr);
// 描画配列への格納
for(i = 0; i < ResultSize; i++){
SIGMAt[i + Delay] = Result[i];
}
// 次回計算の準備
CalcStartTime = 0;
RExecute(R, "last.pred <- tail(result, 1)");
RAssignInteger(R, "fault", 0);
RAssignBool(R, "status", FALSE);
}
}
// タイムフレームの切り替え等でラインが消えた場合の対処(うまく行くとは限らない)
if(SIGMAt[0] == EMPTY_VALUE) LastCalcBar = 0;
// 足が増えていなければ戻る
if(LastCalcBar == Bars) return(0);
// Rに渡す変数や価格配列の準備と実行
// 初回計算時は必要な足すべて、2回目以降は増えた足のみ渡す
if(LastCalcBar == 0){
SampleSize = PastBars + Size;
Loop = PastBars;
}else{
SampleSize = Bars - LastCalcBar;
Loop = SampleSize;
}
RAssignInteger(R, "loop", Loop);
ArrayResize(NewPrice, SampleSize);
for(i = 0; i < SampleSize; i++){
if(PrevClose) NewPrice[i] = Close[i + 1];
else NewPrice[i] = Open[i];
}
RAssignVector(R, "new.price", NewPrice, SampleSize);
RExecute(R, "price <- c(price, rev(new.price))");
CalcStartTime = TimeLocal();
LastCalcBar = Bars;
// GARCHモデル推定関数の実行
// 非同期実行時(RExecuteAsync)、RIsBusy()で計算終了が確認されるまでRに次の命令を送ってはいけない
if(Test){
if(fGarch) RExecute(R, "result <- fg.GARCH.pr(price, size, p, q, loop, last.pred)");
else RExecute(R, "result <- ts.GARCH.pr(price, size, p, q, loop, last.pred)");
}
else if(fGarch) RExecuteAsync(R, "result <- fg.GARCH.pr(price, size, p, q, loop, last.pred)");
else RExecuteAsync(R, "result <- ts.GARCH.pr(price, size, p, q, loop, last.pred)");
return(0);
}
次にRの環境構築と関数のソースファイルです。
ファイル名とその保存場所はインジケータープログラム内の記述と一致させなければなりません。
予測を行う箇所でGARCH(1,1)の計算式をコメントアウトしているのはpredictがうまく行かない場合に備えて。
GARCH(1,1)なら必要になることはおそらくないと思いますが。
# src_GARCHpred
library(tseries)
library(fGarch)
# グローバル変数の準備
online <- TRUE
status <- FALSE
price <- NULL
last.pred <- 0
fault <- 0
#-------------------------------------------------------------------------------+
# fGarchによる推定関数
fg.GARCH.pr <- function(price, size, p, q, loop, last.pred){
# 計算開始の宣言
status <<- FALSE
form <- paste('~garch(', p, ', ', q, ')', sep = '')
# form <- paste('~arma(', m, ', ', n, ') + garch(', p, ', ', q, ')', sep = '')
# 原系列を対数化してから階差を取る、100倍に深い意味なし
dl.price <- diff(tail(log(price), size + loop)) * 100
pred <- rep(NA, loop)
# ローリング
for(i in seq_along(pred)){
temp <- dl.price[i:(i + size - 1)]
last.temp <- tail(temp, 1)
# モデルの推定
g.fit <- tryCatch(garchFit(formula(form), temp, include.mean = F, trace = F),
error = function(err) FALSE, warning = function(warn) FALSE)
# 推定に成功したら予測実行、いずれかに失敗したら値はNA
if(!is.logical(g.fit)){
pred[i] <- tryCatch(as.numeric(predict(g.fit, n.ahead=1, doplot=F)[3]),
error = function(err) NA, warning = function(warn) NA)
# GARCH(1,1)の予測式
# pred[i] <- sqrt(sum(g.fit@fit$matcoef[,1] * c(1, last.temp ^ 2, tail(g.fit@sigma.t, 1) ^ 2)))
}
}
# 推定または予測に失敗した要素の総数をMT4に返すために取得
fault <<- length(which(is.na(pred)))
# NAの要素には前の要素の値を代入
if(is.na(pred[1])) pred[1] <- last.pred
na.result <- !is.na(pred)
pred <- pred[cummax((seq_along(pred)) * na.result)]
# 計算終了の宣言
status <<- TRUE
return(pred)
}
#-------------------------------------------------------------------------------+
# tseriesによる推定関数。やってることはfGarchの場合と一緒
ts.GARCH.pr <- function(price, size, p, q, loop, last.pred){
status <<- FALSE
dl.price <- diff(tail(log(price), size + loop)) * 100
pred <- rep(NA, loop)
for(i in seq_along(pred)){
temp <- dl.price[i:(i + size - 1)]
last.temp <- tail(temp, 1)
g.fit <- tryCatch(garch(temp, order = c(p, q), control = garch.control(trace = F)),
error = function(err) FALSE, warning = function(warn) FALSE)
if(!is.logical(g.fit)){
pred[i] <- tryCatch(as.numeric(predict(g.fit, temp, genuine = T)[length(temp) + 1, 1]),
error = function(err) NA, warning = function(warn) NA)
# pred[i] <- sqrt(sum(coef(g.fit) * c(1, last.temp ^ 2, tail(fitted(g.fit)[,1] ,1) ^ 2)))
}
}
fault <<- length(which(is.na(pred)))
if(is.na(pred[1])) pred[1] <- last.pred
na.result <- !is.na(pred)
pred <- pred[cummax((seq_along(pred)) * na.result)]
status <<- TRUE
return(pred)
}
さて、走らせてみればわかりますが、tseriesに比べてfGarchは計算に10倍くらい時間がかかります。
私の安物PCではサンプル数1000の1000回回しに1分以上待たされます。
2回目以降は増えた足の分しか回さないのですぐですが。
待っている間、不安になるのでDebugViewやタスクマネージャーで様子を窺います。
これがtseriesなら数秒。
もうtseries一択で行きたいところですが、それはそれでちょっと問題がありまして。
デフォルト設定でドル円日足の2007年8月頃から今日までの画像をご覧ください。
緑がtseries、その上に描画したベージュがfGarchです。
ほとんどの期間でほぼぴったり重なっていますが、2012年9月から今年の4月あたりはどうしたことでしょう？
拡大してみるとこんなです。
最尤法最適化の結果で、しかもローリングした1期ずつの描画なので、
どちらが正しいと言える根拠も私は持ち合わせておりませんが、
普通の印象として不自然なのはtseriesの方ではないでしょうか。
いろいろと試してみて分かったのは、サンプルサイズを小さくしていった場合に、
tseriesは前期の推定結果との整合性に欠けるような係数が突然選ばれることがあるということです。
ほぼ1年分の260サンプルだとこのありさま。
でもこれ実はエラーの回数はfGarchの方が多いという。
tseries 12/1000に対し、fGarch 61/1000。
ちなみにサンプル数1000の時はどちらもエラーゼロです。
tseriesの最適化の最大試行回数がデフォルトでは100回なのでこれは少な過ぎ？と考え、
fGarchと同じ10000回に増やしてみましたが、結果はほとんど変わりませんでした。
まったくの素人考えですが、tseriesはfGarchに比べて探索範囲が狭いのかも？
とはいえ初期値を与えるような引数を探してみましたがこれは見当たらず。
さすがに実装を潜るまでの情熱は持てませんでした。
まぁ、速さを追求すれば失うものがあっても仕方ないよね・・・ということで個人的にはtseriesは却下、
fGarchの方で遊んでいたところ、また別の疑問が。
次の画像はどちらもデフォルト設定のfGarchで、ベージュが当日始値で計算したもの、水色が前日終値です。
基本24時間営業の為替市場では週明けの窓以外は前日終値と当日始値はほぼ同じ水準のはずです。
週明けだってそうそう頻繁に大きな窓は開けません。
それでもこれだけの違いが出てくるんですね。
ではtseries(オレンジ)とfGarch(水色)、どちらも前日終値でさっきの問題の期間を見てみると・・・。
あれ？今度はほぼ揃っちゃった？
なんかGARCHって私が思っていたよりセンシティブなものなのでしょうか。
1000本くらい食わせとけばまぁまぁ安定するだろうというのは甘かったのかもしれません。
MT4のチャートだとわかりづらい部分もあるし、ローリングするのと、しないで同じ期間を見るのとではまったく別の話なので、
興味のある方はデータを用意してR上で遊んでみてください。
ローリングしていない結果のプロットですらちょっとしたことでかなり違う印象の絵になることもあり、
何をどう見るかにもよりますが、注意が必要だなと思いました。
あとfGarchのgarchFitはARMA+GARCHの同時推定を行ったり、仮定する分布を標準正規分布以外から選ぶこともできます。
計算所要時間が増えたり、予測も含め失敗する確率が上がったりするし、何よりオーバーフィッティング沼ですけど。
モデル推定ができているのにpredictで失敗するのは関数内部でARMA部分の予測をarima関数でやっているせいかも。
仮定している分布が違うのに係数と新値だけ渡して無理矢理やらせるから？
諸々、実装や記述に誤りがありましたらごめんなさい。
お気付きの点、どうぞご指摘ください。
＜おまけ画像＞ ヒストリカルボラティリティ(20日、青)と一緒に。
