Julia vs Python: ビットコインオプションのモンテカルロシミュレーション


Justin Huang(@rawjustin)さんのJuliaとPythonの比較記事(原題:Julia vs. Python: Monte Carlo Simulations of Bitcoin Options)が非常に面白かったので、翻訳してみました。

Juliaは科学計算・数値計算向けの言語として高速さを売りにしている言語で、Rからの移行を進める人たちも増えているようです。

以前訳した科学計算にPythonが適しているという記事を読んでいた時に、そういえば、Juliaって言語があったな、と思って調べていました。(はじパタでもJuliaのLTしましたし、kawasaki.rb 8回で出てくるベンチマークも印象に残っています)

今回のJustinさんの指摘では、普通のプログラマーにとってPythonよりもJuliaの方が、特定の高速化を意識せずにコードを書けるため、プロトタイプに向いているのでは、という点が興味深い視点でした。

なお、元の文章とコードの権利はJustin氏に(いくつかのPythonのコードはDr. Hilpischに)帰属します。
翻訳を快諾していただいたJustin氏に感謝します。

なお、翻訳の間違い等ありましたら、chezouまでご指摘ください。

[2014/05/07 15:30追記]
Hackers Newsのコメントで vectorizeしないで4倍速くなるコードが紹介されています。

jwmerrill 氏曰く、

In general, carefully written devectorized Julia will be faster than equivalent vectorized code because it is possible to avoid allocating containers at intermediate stages of the computation. This is surprising for people coming from e.g. Matlab, R, or Numpy, because vectorized code is often faster in those environments.

ということなので、Juliaは驚くべきことに、vectorizeしない方が高速になるようです。やばい


私は最近、HackerSchoolに参加した時にStefan Karpinskiに聞いたJuliaを紹介した。
Juliaはとても高速で高性能な科学計算言語をターゲットとしており、その速度はCのネイティブコードに近づきつつある。
Python quants in NYCに参加した後、Dr. Yves J. Hilpischが金融分析におけるPythonの処理速度についてのプレゼンテーションを聞いた時に、JuliaとPython/numpy stackを比較することを決めた。

これはJuliaとPythonの対決ではなく、むしろ同じ問題を2つの言語で解いて比較をする練習である。
ipython/ijuliaのコードを添付したので、よかったら試して欲しい。

Monty Python Bitcoinとは

もしかすると、的はずれなことを言うかも知れないが、やろうとしていることを簡単に説明してみる。
基本的に、古典的なオプション価格問題のモンテカルロシミュレーションをJuliaとPythonの両方で解き、2つの言語で取られているアプローチを比較する。

モンテカルロ法

モンテカルロ・シミュレーションは、確率的なテクニックである。
ランダムな確率分布を使いサンプリングを繰り返すことで、複雑なシステムをシミュレートする。典型的には精度のために極端に大きな数をサンプリングする。

シンプルで概念的な例としては、カジノでクラップ(私がお金を捨てるのに好む方法だ 訳注:Wikipedia)を繰り返しプレイし結果を記録して、あるベットに対する特定のオッズを見つけることが挙げられる。

ビットコイン オプション

私がビットコインオプションと言う時は、ヨーロッパのコールオプション形式の理論的なビットコインのストック・オプションを意味する。なぜなら、ビットコインに関するものはすべて興味深い(し、もっとクリックを稼げる)

幾何ブラウン運動

株価はランダムウォークに従い、予測できないという金融理論がある。
この仮説に従い、モンテカルロシミュレーションによりビットコインのヨーロッパコール・オプションの評価したい。ビットコインの価格変動をランダムウォークでシミュレートする。(今回は、ランダムな連続値をシミュレートする架空の手法として幾何ブラウン運動を考える)

株価モデル

ビットコインの価格が効率的市場仮説(Efficient Market Hypothesis)とマルコフ過程に従うと仮定する。

関連するパラメータ:

S0: 開始価格。ビットコインの現在の価値として扱う

r: 期待されるリターン、または成長率の平均。リスク中立な決済(例:米国債)で参照される

sigma: ボラティリティのパーセンテージ。 このケースで私は後述するbitstampの年刊の計算に基づく

z: ブラウン運動、つまり平均0の正規分布からランダムなサンプル

rsigma は年度ベースで計測している。

これらのパラメータが与えられた時、ビットコインの株価をモデル化するために確率微分法を使う。

dS_t = r S_t\,dt + \sigma S_t\,dZ_t

このモデルに従い、 確立過程 St を幾何ブラウン運動として定義できる。
いくつかの微積分とStに自然対数の導関数の後、Sをtで離散化した式に到達する:

S_t =  S_{t-\Delta t}e^{\left((r-\frac{\sigma^2}{2})\Delta t+\sigma \sqrt{\Delta t}z_t )\right)}

ビットコインのボラティリティを算出する

ビットコインに対する一般的な論点はボラティリティである。
今日ビットコインを買うことのリスクは、明日には半分の価値になっているかもしれないことだ。
rに対する案を出すためには、ビットコインのボラティリティを計算する必要がある。
Bitstampで交換されたビットコインの(分単位でリサンプリングした)価格変動の対数に対する標準偏差を利用した。2013/03/01以降の年率を算出した。

import pandas as pd
import numpy as np

# read in bitstamp USD trade data with columns time, price, volume
bitstamp = pd.read_csv("bitstampUSD.csv" , names=['time', 'price', 'volume'], index_col=0)
bitstamp.index = pd.to_datetime((bitstamp.index.values*1e9).astype(int))

# resample the csv data to minutely data for volume and price
volume = bitstamp.volume.resample('1min', how='sum')
value = bitstamp.prod(axis=1).resample('1min', how='sum')

# obtain the volume weighted average price
vwap = value / volume

# pad data with repeating values
# then grab a years worth of data starting from 01 March 2013
bitstamp_pad = vwap.fillna(method='pad', limit=10)
bitstamp_year = bitstamp_pad.ix['2013-03-01':'2014-03-01']

# calculate the return and volatility
returns = np.log(bitstamp_year / bitstamp_year.shift(1))
returns.std()*np.sqrt(returns.size) * 100

## output: 260.70492740919627

残念ながら、Juliaのツールには時系列データに対するものが現状では全くない。
DataFrames.jlと呼ばれる素晴らしいパッケージがあり、データフレームの処理と変換を行ってくれる。 例:

using DataFrames

bitstamp = readtable("bitstampUSD.csv", colnames = ["time", "price", "amount"])

残念ながら、これが生のJuliaコードで時間周期に対するfillnaresampleのような関数を書かずにできる最大限のことである。

また、DataFrames.jlはpandasのパッケージよりCSVの読み込みは遅いし、インデキシングもサポートされていない。

Pythonでのモンテカルロシミュレーション


これ以降のほとんどのPythonコードはDr. Hilpischのもので彼にクレジットがあることを注意して欲しい。
(ソース)

ここで、sigmaと他のパラメータをPythonで書かれた単純なモンテカルロシミュレーションを実行することで得られる。
ここでのコードでは基本的に、上述した数式で表されるビットコインの価格が取りうる経路の数Iを生成している。

#
# Simulating Geometric Brownian Motion with Python
#
import math
from random import gauss

# Parameters
S0 = 600; # current bitcoin price
r = 0.02; # risk neutral payoff, assumed 2% for this exercise, in reality probably less.
sigma = 2; # extremely high sigma due to spike in bitcoin prices late last year
T = 1.0; # 1 Time cycle
M = 100; # 100 steps
dt = T / M # dt

# Simulating I paths with M time steps
def genS_py(I):
    # initialize array to hold all of our paths
    S = []
    # for each path i to I
    for i in range(I):
        path = []
        # for each step t to M + 1
        for t in range(M + 1):
            if t == 0:
                # append S0, our starting value to the front of the path
                path.append(S0)
            else:
                # take a random normally distributed number z (mean = 0, std = 1)
                # and append it to our current path
                z = gauss(0.0, 1.0)
                St = path[t - 1] * math.exp((r - 0.5 * sigma ** 2) * dt + sigma * math.sqrt(dt) * z)
                path.append(St)
        S.append(path)
    return S

これは、基本的なPython的な経路を作る方法である。10万経路をシミュレーションして、時間を計測するとこうなる:

I = 100000
%time S = genS_py(I)

# CPU times: user 40.8 s, sys: 887 ms, total: 41.7 s
# Wall time: 42.1 s

くだらないミスがないか確認するために、matplotlibで簡単なグラフを描く。

bitcoin_plot

もし、ビットコインの価値が凄い大きさになっていたり0になっていたら、間違いなく調べ直すべきだ。

Juliaでのモンテカルロシミュレーション

今度は、Juliaでやってみよう

#
# Simulating Geometric Brownian Motion with Julia
#

# Parameters
const S0 = 600; # current bitcoin price
const r = 0.02; # risk neutral payoff, assumed 2% for this exercise, in reality probably less.
const sigma = 2; # extremely high sigma due to spike in bitcoin prices late last year
const T = 1.0; # 1 Time cycle
const M = 100; # 100 steps
const dt = T / M # dt

# Simulating I paths with M time steps
function genS_jl(i::Int64)
    S = {}
    for i in 1:i
        path = Float32[]
        for t in 1:(M + 1)
            if t == 1
                push!(path, S0)
            else
                z = randn()
                st = path[t - 1] * exp((r - 0.5 * ^(sigma,2)) * dt + sigma * sqrt(dt) * z)
                push!(path, st)
            end
        end
        push!(S,path)
    end
    return S
end

Juliaの文法はいくつかの注意すべき点を除いてほとんど同じだ:

  • arrayのindexは1始まりである。Pythonは0始まりなのに対して
  • randn()は平均0、標準偏差1の正規分布からランダムな数を生成する組み込み関数である。
  • Juliaの ^() vs Pythonの **push!() vs .append()などの関数
  • グローバル変数に対する定数宣言(10x倍速くなる)
I = 100000
@elapsed S = genS_jl(I)

# 2.243841008

しかしながら、Juliaで書かれた全く同じコードが20倍も速いのだ! これは、Juliaが勝ったことを意味するのだろうか? このアルゴリズムの一実装に対する速度の面においては、Yesだ。しかし、実際にはそうではない。 これは、後述するようにもっとも効率的なモンテカルロシミュレーションのモデル化手法ではない。

ビットコインのコールオプションの評価

さて、演習の一部としてビットコインの価格に対して、理論的なヨーロッパのコール・オプションの評価をしたい。
満期になると、コール・オプションの価値は次式で表すことができる。

max((S_t - K), 0)

ここで、Kはオプションの行使価格とする

オプションの行使価格を K = $1000 とすると、どちらの言語でも単純に計算できる

python

K = 1000.0
C0 = math.exp(-r * T) * sum([max(path[-1] - K, 0) for path in S]) / I
round(C0, 3)

# 361.203

julia

K = 1000.0
C0 = exp(-r * T) * sum([max(path[end] - K, 0) for path in S]) / I
round(C0, 3)

# 374.496


注意:これらの2つの計算における値は、試行回数を増やすと361に収束する。

ここでは、arrayの最後の要素にアクセスする方法([end] vs [-1])以外、言語による違いはほとんど無い。

ベクトル化した数式の処理

しかしながら、ベクトル計算をすることでより高速に処理をすることができる。
ベクトル化のテクニックを使うことで、乗算を1つの要素を何回もループで回すのではなくベクトルに対し一度で処理ができる。
これは、コンパイラやハードウェアによる最適化の恩恵を受けやすく、かなり処理時間が速くなる。

PythonではNumPyというパッケージを使ってベクトル化できる。

import numpy as np
def genS_np(I):
    S = S0 * np.exp(np.cumsum((r - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * np.random.standard_normal((M + 1, I)), axis=0))
    S[0] = S0
    return S

%time S = genS_np(I)

# CPU times: user 1.07 s, sys: 147 ms, total: 1.22 s
# Wall time: 1.38 s

おお、ベクトル化した計算で桁違いに処理が速くなった! 先ほどと似た方法で、オプションの理論値を$1000の行使価格として計算ができる。

C0 = math.exp(-r * T) * np.sum(np.maximum(S[-1] - K, 0)) / I
round(C0, 3)

# 367.457

まったく同じベクトル化がJuliaでもできる。

function genS_jlvec(I)
    S = S0 * exp(cumsum((r - 0.5 * ^(sigma, 2)) * dt + sigma * sqrt(dt) * randn((M + 1, I)), 1))
    S[1,1:end] = S0
    return S
end

@elapsed S = genS_jlvec(I)

# 1.430265529

JuliaのメリットとしてはNumPyのようなライブラリーが不要だということだ。 我々は素のJuliaを使っている。 驚くべきことに、Juliaのナイーブなloopでの実装は、ベクトル化された実装と同程度に速い! これは、アルゴリズムの単純なプロトタイピングに非常に向いていることを示す。 ある種の計算でナイーブな実装を駆逐することで猛烈な高速化を期待する、Pythonが提案できない価値のある論点である。

最後に、Juliaでオプションの値を計算する文法はほとんど同じである:

K = 1000.0
C0 = exp(-r * T) * sum(max(S[end, 1:end] - K, 0)) / I
round(C0, 3)

# 359.783

結論

Python Loops: 42.1 s
Julia Loops: 2.24 s
Python Vectorized: 1.38 s
Julia Vectorized: 1.43 s

speeds

Juliaは、科学/数値計算にとってとても大きなポテンシャルを持っているように思える。
ナイーブなループを使った、ビットコインの価格のモンテカルロの経路シミュレーションはPythonよりもかなり速く、ベクトル化して最適化された計算と同じオーダーの速度だった。
つまり、今回の分析での重要な知見だと考えている:
Juliaは計算的に困難な問題を効率的に解くにあたって、ベクトル化された処理を必要としない。

Juliaはベクトルや行列の計算にNumPyのような外部ライブラリを必要としない。
他の化学計算言語が険しい学習曲線であるのとは異なり、Juliaは普通のプログラマーにとって極めて読みやすくなる。
処理速度は前述のとおり明確に速い。

残念ながら、DataFrames.jlなどJuliaのデータ分析のためのツールは、Pythonのpandasと比較して性能面、機能面ともに成熟してるとも安定しているとも言いがたい。
Juliaのコミュニティは若く成長している。
だから、様々な人々が言語に対して貢献することで、私はこれらの点で多くの改善や成長がなされることを期待している。しかし、それまでは私はPythonとpandasの行く末を見ていく。

すべてを考慮した上で、メインストリームとして採用されることを約束はできないが、私はJuliaを試すことをおすすめする。
私としてはJuliaのコミュニティに貢献したいと強く思っている。

Files

julia code
python code
bitcoin code

Thanks

Dr. Yves J. Hilpisch (inspiration)
Neeraj Wahi (options pricing)
James J. Porter (Julia)

カテゴリー: program, translation | タグ: , , | 2件のコメント

kawasaki.rb #011 を開催しました #kwskrb


4/23にkawasaki.rb #011を開催しました。
次回で1年になるのかと思うと、感慨深いです。
なお、kawasaki.rbは毎月第4水曜の19:30-開催しています。(最終週だと思われていたのでご説明)

togetter

パーフェクトRuby読書会

2-9 例外 〜 2-12 組込みの変数/定数までを読みました。

組込み変数の、$-Fはワンライナーの-Fと互換だよねー、といった話や、RailsとMySQLでprecisionを使おうとしたら予約語でハマったという話が出ました。

次回は、ついに3章に突入です。

AngularJS+SinatraでTurbolinks: @snowcrush さん

@snowcrush さんのblogは、AngularJS + Sinatraで構成されているそうですが、
google botsはJSで動的に生成されたページにindexを貼ってくれないので、
なんとかして拾えるようにしましたという話。
“これってTurbolinksじゃん!Railsで良いじゃん!”という結論には、笑ってしまいました。

カテゴリー: Ruby, 勉強会 | タグ: , , | コメントをどうぞ

TokyuRuby会議07でLTしてきました #tqrk07


無限プレモルが出るという危険な楽しいイベント、TokyuRuby会議07rubyistokeitimerにドラの音を付けたよ!という話をしました。

http://rubyistokei.herokuapp.com/timer?limit=1&gong=on
こんな感じでご家庭でもドラの音を楽しみながらLTできます!

本当はデモをするつもりだったのですが、いい感じに出来上がっていたので肝心の音を鳴らすことができませんでした。せっかく、iPhone, iPad, Androidで動作確認をして会場で準備までしてたのに!

あと、同僚に教えてもらった大事な情報として、フリーの音素材を探すのにはfreesoundがいいよ!ということをお知らせしそびれたのが悔やまれます。

大抵、Webで音を鳴らす場合再圧縮したり切ったり貼ったりすることが多いので、加工をできる音素材が必要なのですが、
自作のものでない場合、ライセンス的にどうなのか不安になると思います。

そんな時、CC0ライセンスの音素材を探せるfreesoundはめっちゃ助かりました。

ちなみに、持ち込み飯はキッシュを作りました(レシピはそのうち投稿します)。ラップにつつまずホールで持っていけばよかった。

切る前はこんな感じでした

### 青梗菜とベーコンのキッシュ

ツナと長ネギの味噌キッシュ

なお、飯王の写真を見せた所、妻から次回飯王を目指せとの指令が下りました。

カテゴリー: 勉強会 | タグ: , , | コメントをどうぞ

kawawaki.rb #010 を開催しました #kwskrb


kawasaki.rbもいよいよ10回目を開催してきました。
まさかの1日二回LTという展開でしたが、無事に話すことが出来て良かったです。

togetterはこちら

パーフェクトRuby読書会

2-7-7 Regexp, 2-7-8 %記法の括弧, 2-7-9 手続きオブジェクト(Proc), 2-8 様々な代入式まで読みました。

はじめて翻訳記事を書いたら300ブクマ超えた話 (chezou)

前々回の話を翻訳した時に思ったことを書きました。

Jekyll ドキュメントの日本語翻訳レポジトリを作った話(@kk_Atakaさん)

@kk_Atakaさんによる、Jekyllのドキュメントの翻訳の話。
ご本人のブログの記事をもとにお話していただきました。
Jekyllドキュメントの日本語翻訳リポジトリ「jekyllrb.com.jp」を作成しました – kk_Atakaの日記

特に、本家のリポジトリが英語しか想定しておらず、「言語ごとの翻訳を本体に入れるつもりはない」と言われているため、
本家との差分をどう管理していけば良いのか、という話を議論しました。
懇親会で、@merborneさんのTogglateを使って本体の更新を反映させやすくできると良いですねーと盛り上がりました。

Railsのhelperでnamed capture使ったらハマった話 (chezou)

当日人気のない話でしたが、ActiveSupport::SafeBufferはRegexp.last_matchがそのままだと取れないので気をつけよう、という話でした。

次回は、4/23(水)に開催予定です。
申し込みは下記からおねがいします。

Kawasaki.rb #011 – Kawasaki.rb | Doorkeeper

カテゴリー: Ruby, 勉強会 | タグ: , , | コメントをどうぞ

4ステップでDoorkeeperにコメント欄を追加する方法


画像

Doorkeeperはコミュニティなどのイベントを行うのに便利ですよね。
そんなDoorkeeperですが、他のzusaarなどのサイトと比較してコメント欄がデフォルトでついていないので、
Disqusと連携してコメント欄を追加する方法を書きます。

僕がやった時は日本語で情報が見つからなかったのですが、
@yumu19ですらはまっていたので、結構挫折している人がいるのかも…、と思いブログに書くことにします。

1. Disqusのアカウントをとる

http://disqus.com/

アカウントの取り方はメモを取っていないのですが、多分普通にメール認証すれば良いはず…。

2. “Add Disqus to Your Site”をクリック

ログインした後のtop右下の”Add Disqus to Your Site”をクリックする

画像

3. 必要な情報を埋めていく

“Site name”はdescriptionなので適当に決めればOK。

“Choose your unique Disqus URL”で入力する情報が、後のDoorkeeperで入力する情報。
“Category”も適宜入力すればOK

“Finish registration”を押せば、完了。

その後に出てくるscriptの生成画面はdoorkeeperには関係ないので、気にしなくて良い。

画像

4. Doorkeeperのコミュニティ管理画面の”連携機能”を開く

“Disqus shortname”に、さっき作ったDisqus URLの.disqus.comの前の文字列を入れる。
今回の例だと”kawasakirbtest”を入れて保存する

画像

 

これで、冒頭のようなコメント欄がイベントページの下部にできます!

カテゴリー: 雑談 | タグ: , | コメントをどうぞ

kawasaki.rb #009を開催しました #kwskrb


去る2/26(水)にkawasaki.rb #009を開催しました。

togetterのまとめはこちらです。
kawasaki.rb #009 まとめ #kwskrb – Togetterまとめ

パーフェクトRuby読書会

2-7-4 配列, 2-7-5 ハッシュ, 2-7-6 範囲(Range)をやりました。

Hashでは、自分で作ったクラスをHashのkeyに用いる場合はどうすればいいのか、という話がでました。
自作クラスにhashメソッドとeql?(other)メソッドを定義すれば、
自作クラスでも同値性判定が行われるという確認をしました。
(やっててよかった、レシピブック&yokohama.rb…)

コードで書くと、こういう感じですね。

class Foo
  attr_reader :a, :b

  def initialize(a, b)
    @a, @b = a, b
  end

  def hash
    @a%10 + @b%7
  end

  def eql?(obj)
    @a == obj.a && @b == obj.b
  end
end

f1 = Foo.new(1, 2)
f2 = Foo.new(1, 2)

f1.eql?(f2) # => true

h = {}
h[f1] = "foo"
h[f2] # => "foo"

Rangeでは、(1..-1)のような負の数までのRangeってなんであるんだろう?という話をしました。
使用するケースとしては、Arrayにアクセスするときary[1..-1]というのを渡すよねー、という話になりました。

LT

BestGems.org (ぺけみさおさん)

LTはぺけみさおさんによる ランキングから見るRubyGems -BestGems.orgのご紹介- と、@aflc_jp さんによるpandasの紹介(資料公開されたら追記します)が発表されました。

ぺけみさおさんの、BestGems.orgはRubyGems.orgを毎日クローリングして、
DL数の変化を見ているというシンプルながら価値のありそうなアプローチでランキングを作られていました。
実際にランキング作って分析するのって、楽しいですよね。

ランキングを分析しているとnaughtというgemの人気が上がっているとのことでした。

pandas紹介するよ(aflc_jpさん)

aflc_jpさんのpandasの紹介は、Pythonのデータ分析によく使われるライブラリである
pandasの基本と使い方のexampleを紹介してもらいました。

資料
デモ

説明を聞くと、pandasはRでよく使われているDataFrameを扱えるようにしているところが特徴のようで、
SQL likeにgroupbyして集計をしやすくなっているという話でした。
例としてはオリンピックのメダルのデータを使って、日本が歴代メダル数を稼いだ種目などを表にしていました。
Excelに出力することもできるので、分析はpandasでreportingはexcelで、という使い分けができるようです。

改めて実感したのは、IPython notebookは便利だなーということ。
IJuliaで便利さは実感していたのですが、Rubyの実行もブラウザベースでやれるようになるそうなので、
kawasaki.rbのパーフェクトRuby読み会でもIPythonベースにしようかなー。

次回は、3/26(水)に開催予定です

カテゴリー: 勉強会 | タグ: , , , | コメントをどうぞ

「はじめてのパターン認識」読書会に参加しました #はじパタ


先日開催された第12回「はじめてのパターン認識」読書会@yamakatuさんに誘われるままに参加してきました。

部分空間法の話の続きとカーネル主成分分析分析のお話が中心でした。
この本、タイトルはやさしそうな雰囲気を出していながら、普通に数式をきちんと書いてくれている本で、わりときちんと読まないと読み落とす感じがたまりません。
でも、こうやってわりと網羅的に書いてくれる本って一昔前にはなかったので、自分の中の情報を整理する意味でも良いかなーと思いました。

また、「Juliaを使った機械学習」という釣り気味のタイトルでLTしてきました。

実際には、Gadfly使うとおしゃれなグラフが書けますよーというレベルでしかなかったのですが、D3.jsなんかとも親和性高い描画エンジンなのでよければ使ってみると楽しいですよ。

カテゴリー: program | タグ: , | コメントをどうぞ