RigelのR言語メモであーる(R言語だけとは言っていない)

RigelのR言語メモであーる(R言語だけとは言っていない)

興味あることや趣味、やったことについて書くよ

BTCFX 仮想通貨 レンジ相場の統計モデリング その1

仮想通貨のbotを作りたいです。

トレンドフォロー系の論文を見つけて行けそうと思い、実装してみました。

が、あまりうまく行かず。

ネットをあさっていると、「レンジ相場がほとんど」という記事を見つけ、

レンジ相場でアルファを見つけたいと思いました。

ということで、レンジ相場を統計モデリング

結論としてはうまく行かなかったんですが、メモとして残します。

レンジ相場

f:id:strawberry_kyon:20190323121046p:plain

こういったものを想定します。

これはbitflyerのBTCJPY、の日足終値です。

2018-02-07から2018-05-17まで、100本分。

統計モデル

レンジ相場を見ていると、三角波に近いと思いました。

ということで、ベースは三角波から出発します。

具体的には、以下のような関数です。

{\begin{eqnarray}
f_{tri}(t | T, \phi)=\frac{8}{\pi^2}\sum_{d=1}^D \frac{\sin((2d-1)\frac{\pi}{2})}{(2d-1)^2} \sin((2d-1)2\pi (\frac{t}{T} - \phi) )
\end{eqnarray}}

三角波フーリエ級数展開した形です。

偶数倍音は0となって無駄な計算になるため、予め奇数倍音だけ足し合わせています。

現時点でのパラメータは周期{ \displaystyle T}と、初期位相{ \displaystyle \phi}の2つです。

{ \displaystyle D}は何倍音まで考慮するかで、ハイパーパラメータ的な感じです。

{ \displaystyle T=100, \phi=0.5, D=10}のときは以下のような感じになります。

f:id:strawberry_kyon:20190323123528p:plain

{ \displaystyle T=100}なので、{ \displaystyle 100}点で1周期。

{ \displaystyle  \phi=0.5}なので、半周ずれる。

三角波フーリエ級数展開の近似が良いので{ \displaystyle D}はもっと小さくても良いとは思う。

続いてこの波形の倍率を変えたり、縦軸に並行移動したりして、変形します。

{\begin{eqnarray}
f_{rng}(t | T, \phi, a_{ini}, a_{fin}, b_{ini}, b_{fin})=\frac{(N-t)a_{ini}+ta_{fin}}{N}f_{tri}(t | T, \phi) + \frac{(N-t)b_{ini}+tb_{fin}}{N}
\end{eqnarray}}

ここで、{ \displaystyle  N}は観測値の長さとしておきます。

なんだかややこしいのは私の変数の置き方が良くないからです。

{ \displaystyle T=100, \phi=0.5, D=10, }

{ \displaystyle a_{ini}=2000, a_{fin}=500, b_{ini}=410000, b_{fin}=408000}のときは以下のような感じです。

f:id:strawberry_kyon:20190323125902p:plain

{ \displaystyle  N=500}としてあります。

いままでのパラメータを図示すると以下のような感じです。

f:id:strawberry_kyon:20190323131035p:plain

最後に、誤差は正規分布を仮定しておきます。

つまり、観測値を{ \displaystyle  y_t}とすると、

{\begin{eqnarray}
y_t = f_{rng}(t | T, \phi, a_{ini}, a_{fin}, b_{ini}, b_{fin}) + \epsilon_t  \\
 \epsilon_t 〜 iidN(0,\sigma^2)
\end{eqnarray}}

こんな感じに。

パラメータ推定

上記のモデルの6つのパラメータを推定します。

正規分布を仮定しているので、最尤推定するには最小二乗法をすれば良いですね。

だけど、非線形で複雑なので直感的に解析的には求められないと思われます。

さらに、周期と初期位相があるので局所解が多数出てくると思われます。

つまり、初期値にかなり依存します。

試行錯誤した結果、以下の流れでいい感じのパラメータが得られました。

  1. 線形回帰して、その係数から{ \displaystyle  b_{ini}, b_{fin}}をおおよそ推定する。
  2. 観測値のレンジから、{ \displaystyle  a_{ini}, a_{fin}}をおおよそ推定する。
  3. 周期と初期位相の初期値のグリッドを作成し、上記のパラメータと合わせて二乗誤差を最適化する。

曖昧な書き方ですが、こんな感じです。

ちょっとずつ最適化していくことで、初期値問題を回避していい感じのパラメータが推定することができました。

結果

冒頭のチャートに三角波を当てはめた結果、以下のようになりました。

f:id:strawberry_kyon:20190323202957p:plain

なんかパット見はそれっぽいですよね(それっぽくなる部分を選びました)。

赤の実線は、上記統計モデルを当てはめたもの。

赤の破線は、上記統計モデルの上側を結んだ直線と、下側を結んだ直線。

ちなみに、その後はこんな感じです。

f:id:strawberry_kyon:20190323204038p:plain

うーん、使い物になりません。

まとめ

レンジでアルファを作るなら、ここまで正確なモデルは必要ないですね。

今回のモデルが微妙だった理由としては、周期が固定されている点ですかね。

そのせいで、いい感じの支持線を推定することはできなかったと思います。

上側と下側の支持線に反発されるのがレンジ相場。

その間の周期が固定というのは今思うとセンスがない。

上側で連続で反発されたあと下側で反発されたり。

つまり、支持線の推定とかってことになりそうです。

f:id:strawberry_kyon:20190323204943p:plain

こんな感じの支持線が推定できると嬉しい気がしますね。

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)