BTCFX 仮想通貨 レンジ相場の統計モデリング
仮想通貨のbotを作りたいです。
トレンドフォロー系の論文を見つけて行けそうと思い、実装してみました。
が、あまりうまく行かず。
ネットをあさっていると、「レンジ相場がほとんど」という記事を見つけ、
レンジ相場でアルファを見つけたいと思いました。
ということで、レンジ相場を統計モデリング。
結論としてはうまく行かなかったんですが、メモとして残します。
統計モデル
レンジ相場を見ていると、三角波に近いと思いました。
ということで、ベースは三角波から出発します。
具体的には、以下のような関数です。
偶数倍音は0となって無駄な計算になるため、予め奇数倍音だけ足し合わせています。
現時点でのパラメータは周期と、初期位相の2つです。
は何倍音まで考慮するかで、ハイパーパラメータ的な感じです。
のときは以下のような感じになります。
なので、点で1周期。
なので、半周ずれる。
三角波のフーリエ級数展開の近似が良いのではもっと小さくても良いとは思う。
続いてこの波形の倍率を変えたり、縦軸に並行移動したりして、変形します。
ここで、は観測値の長さとしておきます。
なんだかややこしいのは私の変数の置き方が良くないからです。
のときは以下のような感じです。
としてあります。
いままでのパラメータを図示すると以下のような感じです。
最後に、誤差は正規分布を仮定しておきます。
つまり、観測値をとすると、
こんな感じに。
パラメータ推定
上記のモデルの6つのパラメータを推定します。
正規分布を仮定しているので、最尤推定するには最小二乗法をすれば良いですね。
だけど、非線形で複雑なので直感的に解析的には求められないと思われます。
さらに、周期と初期位相があるので局所解が多数出てくると思われます。
つまり、初期値にかなり依存します。
試行錯誤した結果、以下の流れでいい感じのパラメータが得られました。
- 線形回帰して、その係数からをおおよそ推定する。
- 観測値のレンジから、をおおよそ推定する。
- 周期と初期位相の初期値のグリッドを作成し、上記のパラメータと合わせて二乗誤差を最適化する。
曖昧な書き方ですが、こんな感じです。
ちょっとずつ最適化していくことで、初期値問題を回避していい感じのパラメータが推定することができました。
結果
冒頭のチャートに三角波を当てはめた結果、以下のようになりました。
なんかパット見はそれっぽいですよね(それっぽくなる部分を選びました)。
赤の実線は、上記統計モデルを当てはめたもの。
赤の破線は、上記統計モデルの上側を結んだ直線と、下側を結んだ直線。
ちなみに、その後はこんな感じです。
うーん、使い物になりません。
まとめ
レンジでアルファを作るなら、ここまで正確なモデルは必要ないですね。
今回のモデルが微妙だった理由としては、周期が固定されている点ですかね。
そのせいで、いい感じの支持線を推定することはできなかったと思います。
上側と下側の支持線に反発されるのがレンジ相場。
その間の周期が固定というのは今思うとセンスがない。
上側で連続で反発されたあと下側で反発されたり。
つまり、支持線の推定とかってことになりそうです。
こんな感じの支持線が推定できると嬉しい気がしますね。
Rコード
雑ですが投下しておきます。
なんかミスってるような気もするけど気にしない。
library(data.table) library(magrittr) library(ggplot2) library(dplyr) library(tidyr) library(scales) tri=function(n_point,period,phi,dim){ tri0=function(n_point,period,phi,dim){ return(sin((2*dim-1)*pi/2)*sin((2*dim-1)*2*pi*((1:n_point)/period-phi))/(2*dim-1)^2) } return(apply(sapply(1:dim,tri0,n_point=n_point,period=period,phi=phi),1,sum)*8/pi^2) } sim_tri=function(n_point,period,phi,a_ini,b_ini,a_fin,b_fin,dim){ return(seq(a_ini,a_fin,length=n_point)*tri(n_point=n_point,period=period,phi=phi,dim=dim)+seq(b_ini,b_fin,length=n_point)) } opt_tri=function(data,dim){ x=data[,2] tmp=predict(lm(x~t,data.frame(x,t=1:length(x))),data.frame(t=c(1,length(x)))) b_ini_i=tmp[1] b_fin_i=tmp[2] tmp=diff(range(x))/2 a_ini_i=tmp a_fin_i=tmp para_mat=as.matrix(expand.grid(data.frame(period_i=seq(0.2,1,by=0.2)*length(x),phi_i=tail(seq(0,1,length=6),-1)))) op=function(para_ini){ res=optim(c(para_ini[1],para_ini[2]),function(para,data,a_ini,b_ini,a_fin,b_fin){ period=para[1] phi=para[2] return(sum((data-sim_tri(n_point=length(data),period=period,phi=phi,a_ini=a_ini,b_ini=b_ini,a_fin=a_fin,b_fin=b_fin,dim=dim))^2)) },data=x,a_ini=a_ini_i,b_ini=b_ini_i,a_fin=a_fin_i,b_fin=b_fin_i) return(res$value) } para_mat %>% apply(1,op) %>% data.frame(para_mat,value=.) %>% arrange(value,period_i,phi_i) ->res0 res0=optim(c(res0[1,1],res0[1,2]),function(para,data,a_ini,b_ini,a_fin,b_fin){ period=para[1] phi=para[2] return(sum((data-sim_tri(n_point=length(data),period=period,phi=phi,a_ini=a_ini,b_ini=b_ini,a_fin=a_fin,b_fin=b_fin,dim=dim))^2)) },data=x,a_ini=a_ini_i,b_ini=b_ini_i,a_fin=a_fin_i,b_fin=b_fin_i) period_i=res0$par[1] phi_i=res0$par[2] res=optim(c(period_i,phi_i,a_ini_i,b_ini_i,a_fin_i,b_fin_i),function(para,data){ period=para[1] phi=para[2] a_ini=para[3] b_ini=para[4] a_fin=para[5] b_fin=para[6] return(sum((data-sim_tri(n_point=length(data),period=period,phi=phi,a_ini=a_ini,b_ini=b_ini,a_fin=a_fin,b_fin=b_fin,dim=dim))^2)) },data=x) return(list(data=data,res=res,dim=dim)) } plot_tri=function(res){ data=res$data dim=res$dim res=res$res data=data.frame(data,fit=sim_tri(n_point=nrow(data),period=res$par[1],phi=res$par[2],a_ini=res$par[3],b_ini=res$par[4],a_fin=res$par[5],b_fin=res$par[6],dim=dim)) data=data.frame(data, upper=res$par[4]+res$par[3]+((res$par[6]+res$par[5])-(res$par[4]+res$par[3]))/nrow(data)*(1:nrow(data)), lower=res$par[4]-res$par[3]+((res$par[6]-res$par[5])-(res$par[4]-res$par[3]))/nrow(data)*(1:nrow(data))) ggplot(data,aes(x=time))+ geom_point(aes(y=close))+ geom_line(aes(y=close))+ geom_line(aes(y=fit,color="fit"))+ geom_line(aes(y=upper,color="fit",type=2),linetype="dashed")+ geom_line(aes(y=lower,color="fit",type=2),linetype="dashed")+ scale_x_datetime(breaks = date_breaks("week"), labels = date_format("%Y-%m-%d"))+ theme(axis.text.x = element_text(angle = 90, hjust = 1),legend.position="none") } N=500 close=sim_tri(n_point=N,period=100,phi=0.5,a_ini=2000,a_fin=500,b_ini=410000,b_fin=408000,dim=10)+rnorm(N,sd=300) time=seq.Date(as.Date("2018-01-01"),by="days",length=N) %>% as.POSIXct() plot(time,close) data=data.frame(time,close) res=opt_tri(data,dim=10)