Rubyでカルマンフィルタを実装し、株価に適用する

このエントリーをはてなブックマークに追加

前回はVIDYAという指標を使ってシステムの検証を行っていました。この指標は次のようにかけました。

VIDYAt = VIDYAt-1 + a (CLOSE - VIDYAt-1)

このaには短期分散と長期分散の比率が入っていて、比率が高いほど誤差が反映される仕組みになっています。ですので、終値がVIDYAを頻繁に交差するような場合は、そのたびにVIDYAが方向を変えます。それは少し困るので、もうちょっと賢いフィルタを導入することで、どうにかできないかと思いまして、カルマンフィルタを実装していくことにしました。カルマンフィルタを導入すると上記の式は、

予測値t = 予測値t-1 + a (観測値 - 予測値t-1)

みたいな感じでかけて、aは予測値と観測値のそれぞれの誤差の内分比で評価されることになります。例えば株価の場合は観測値に終値を用い、予測値には移動平均を用いて見ることなんてことが出来るかもしれません。これまでの指標と同様の式の形をしていますので、予測値の誤差が大きい場合は観測値の値が反映され、観測値の誤差が大きい場合は予測値の値が反映されることになります。

詳細な説明は他の多くのサイトをご覧ください。なんでみんなあんなに上手に説明できるんだろう。。。

目次

実践

実装

愚直にWikipediaをコーディングしていきます。カルマンフィルターの項に予測と更新の節があり、そこに更新式が説明されています。いかつい式ですが、コードにしてしまえばさっぱりしていますね。行列の計算はNumo::NArrayを利用させてもらいましたが、逆行列のためのLinalgが導入できない悲しい環境でしたので、そこだけMatrixを使っています。

class Kalman
  include Numo

  attr_accessor :ft0, :pt1, :qt0, :gt0, :xht1, :ht0, :rt0

  def predict(ut0)
    ut0 = 0 if not ut0
    @xht0 = @ft0.dot(@xht1) + ut0
    @pt0 = @ft0.dot(@pt1).dot(@ft0.transpose)+@gt0.dot(@qt0).dot(@gt0.transpose)
    [@xht0, @pt0]
  end

  def update(zt0)
    @xht1 = @xht0 + kt0.dot(et0(zt0))
    ktht = kt0.dot(@ht0)
    @pt1 = (Int32.eye(*ktht.shape) - ktht).dot(@pt0)
    [@xht1, @pt1]
  end

  def observe(zt0, ut0=nil)
    predict(ut0)
    update(zt0)
  end

  def et0(zt0)
    zt0 - @ht0.dot(@xht0)
  end

  def st0
    @rt0 + @ht0.dot(@pt0).dot(@ht0.transpose)
  end

  def kt0
    @pt0.dot(@ht0.transpose).dot(Matrix[*st0.to_a].inv.to_a)
  end
end

サンプル1

まずはWikipediaにある、設定例を参照して1次元のトロッコ問題の設定を行っていきます。設定っていうのは行列F、G、Q、H、R、Pを決めて行くことです。これらの行列は式が決まらない限り設定できませんので、Wikipediaを見ながらご確認ください。 コードにするとこのようになります。

Bundler.require
include Numo

sa = 0.1
sz = 0.1
dt = 0.01

dir = File.expand_path '../../data/kalman/', File.dirname(__FILE__)
FileUtils.mkdir_p dir
file = dir + '/sample1.jpeg'

f = DFloat.new(2,2).fill 0
f[0,true] = [1, dt]
f[1,true] = [0, 1]

g = DFloat.new(2,1).fill 0
g[0,0] = dt**2 / 2
g[1,0] = dt

q = DFloat.new(2,1).fill 0
q[0,0] = sa
q[1,0] = sa

h = DFloat.new(1,2).fill 0
h[0,true] = [1,0]

r = DFloat.new(1,1).fill 0
r[0,0] = sz


kalman = Kalman.new
kalman.ft0 = f
kalman.gt0 = g
kalman.qt0 = q
kalman.ht0 = h
kalman.rt0 = r

これで設定が完了です。続いて実際に乱数を生成してどのようにフィルタリングされるかを確認します。

obs = []
trs = []
inp = DFloat.new(1,1).fill 0
x = []
1000.times do |i|
  inp[0,0] = Random.rand
  out, p = kalman.observe inp
  obs << inp[0,0]
  trs << out[0,0]
  x << i
  puts "#{i}\t#{inp[0,0].round(2)}\t#{out[0,0].round(2)}"
end

Numo.gnuplot do
  reset
  set terminal: 'jpeg'
  set output:  file
  set grid: true
  plot [x, obs, with: :lines, title: 'observe', lc: "'salmon'"],
    [x, trs, with: :lines, title: 'true', lc: "'blue'", lt: 0]
end

obsが観測値でtrsが真の値などと呼ばれる予測値になります。Inpは0〜1のランダム値でそれを位置情報として観測しています。観測は誤差のため、かなりギザ付いた感じなんですね。その観測誤差を考慮して実際の位置を予測した結果が真の値として帰ってきます。実際にグラフを確認してみます。

profit_histgram

なんとなーくそれっぽい結果が帰ってきてますね。はじめは予測値の誤差がでかいため、大きく振れてますが次第にそれっぽい位置に落ち着いています。

サンプル2

株価の場合はどうするのかと言うことになるので、適当にモデルを作って確認してみましょう。

Averaget = Averaget-1 + ut

Close`t = Averaget + vt

ut = a(Closet-1 - Averaget-1)

Averageは予測でCloseは観測、制御ベクトルutは1期前のAverageとCloseの差です。vtはノイズです。これでF,Hは要素が1のみの1行1列の行列になることがわかります。G,R,Qは適当に設定して、下のコード のようになりますので、それを走らせて確認してみます。

Bundler.require
include Numo
include Kabu

sa = 1
sz = 50
alpha = 0.1
count = 200

dir = File.expand_path '../../data/kalman/', File.dirname(__FILE__)
FileUtils.mkdir_p dir
file = dir + '/sample2.jpeg'

f = DFloat.new(1,1).fill 0
f[0,0] = 1

g = DFloat.new(1,1).fill 0
g[0,0] = 1

q = DFloat.new(1,1).fill 0
q[0,0] = sa

h = DFloat.new(1,1).fill 0
h[0,0] = 1

r = DFloat.new(1,1).fill 0
r[0,0] = sz


com = Company.find_by_code 'I201'
closes = Soks.parse(com.soks.order(:date).limit(count+10), :close)

kalman = Kalman.new
kalman.ft0 = f
kalman.gt0 = g
kalman.qt0 = q
kalman.ht0 = h
kalman.rt0 = r

kalman.pt1 = DFloat.new(1,1).fill closes[0..10].dev(11)[-1]
kalman.xht1 = DFloat.new(1,1).fill closes[10]

obs = []
trs = []
inp = DFloat.new(1,1).fill 0
ut = DFloat.new(1,1).fill 0
x = []
closes[11..-1].each_with_index do |c,i|
  inp[0,0] = c
  ut[0,0] = alpha * (c - closes[i-1])
  out, p = kalman.observe inp, ut
  obs << inp[0,0]
  trs << out[0,0]
  x << i
  puts "#{i}\t#{inp[0,0].round(2)}\t#{out[0,0].round(2)}"
end

Numo.gnuplot do
  reset
  set terminal: 'jpeg'
  set output:  file
  set grid: true
  plot [x, obs, with: :lines, title: 'observe', lc: "'salmon'"],
    [x, trs, with: :lines, title: 'true', lc: "'blue'", lt: 0]
end

profit_histgram

移動平均とのクロスに対して、鈍感になっていることは確認できますね。

サンプル3

最後に重みを予測値にして、AR(1)の重みを推定させてみます。AR(1)は

Xt = w1 Xt-1 + w0 + u

と書け、wiが重み、uはノイズです。この場合は行列Hに過去の株価を設定し、重みwを予測します。株価の前日比をrtとして

xt = xt-1

rt = w0 + w1 rt-1 + ut

xt = [w0, w1]T

こんな感じに書けるので、Fが2行2列の単位行列で、Hは1行2列の、要素に1と前日の前日比が入った行列であれば良いことがわかります。Q,R,Gは式を満たすように適当に設定します。 コードはこのようになります。

Bundler.require
include Numo
include Kabu

sz = 0.01
sa = 1
p_ = 1
count = 200

dir = File.expand_path '../../data/kalman/', File.dirname(__FILE__)
FileUtils.mkdir_p dir
file = dir + '/sample3.jpeg'

f = DFloat.new(p_+1,p_+1).eye

g = DFloat.new(1,1).fill 0

q = DFloat.new(1,1).fill 0

h = DFloat.new(1,p_+1).fill 0

p = DFloat.new(p_+1,p_+1).eye

r = DFloat.new(1,1).fill 0
r[0,0] = sz

x = DFloat.new(p_+1,1).fill(0)

com = Company.find_by_code 'I201'
closes = Soks.parse(com.soks.order(:date).limit(count+10), :close)
logs = closes.log

kalman = Kalman.new
kalman.ft0 = f
kalman.gt0 = g
kalman.ht0 = h
kalman.qt0 = q
kalman.rt0 = r
kalman.pt1 = p
kalman.xht1 = x

obs = Soks.new
trs = Soks.new
inp = DFloat.new(1,1).fill 0
ut = DFloat.new(1,1).fill 0
x = []
logs.each_cons(p_+1).to_a.each_with_index do |c,i|

  kalman.ht0[0,true] = [1] + c[0..-2]
  inp[0,0] = c[-1]

  out, pt = kalman.observe inp
  obs << inp[0,0]
  trs << out[0]
  x << i
  puts "#{i}\t#{inp[0,0].round(2)}\t#{out[true,0].to_a.join("\t")}"
end

Numo.gnuplot do
  reset
  set terminal: 'jpeg'
  set output:  file
  set grid: true
  set y2tics: true
  set ytics: :nomirror
  plot [x, obs[40..-1].cumu, with: :lines, title: 'observe', lc: "'salmon'"],
    [x, trs[40..-1], with: :lines, axes: :x1y2, title: 'true', lc: "'blue'", lt: 0]
end

私が興味があるのは、w0のドリフト項ですので、それと前日比の累積値をプロットしてみます。

profit_histgram

うーん。これはなんか微妙ですねー。

まとめ

とりあえず作って動かしてみましたよ。ただ、ぼんやりとした評価しかできないあたりに、作っては見たもののってところの限界を感じます。バックテストに落としこんで結果が出れば早いのですが、今のところ、で?って感じにしかなりません。

この辺をかじりだすと時系列統計学なんてあたりに手を出さなきゃならんのですかね。でも、それはシステムにどう落とし込めばいいのかが分からないので結局行き詰まりになりそうな気もするところです。

Recent Entries
Categories
    Tags
    Archives
    Search