8000万レコードをInsertする

概要

・とある物件で8000万弱のレコードを新規のDBに効率的にぶち込む方法を考える必要に
・MySQL(Amazon RDS) + Python(mysql.connector)を想定
・mysql.connector + pythonについては拙文ながら(http://nekopuni.holy.jp/?p=927)に書いております。
・コミットの位置には気をつけよう + Multiple Insert最強ねというお話。

方法その1

まずは最初にやった方法。1レコードごとにInsertしてコミットしていく方法。
今までDB関連でやったコードはレコード量も大したことなかったので以下の方法でも特に問題なかった。

具体的な環境としてはテキストデータ(csv)を読み込み、それをDBにInsertしていくというもの。
csvの中身は日付(DATE)とデータ値(VALUE)がカンマ区切りになっているものであり、エッセンスだけ抜き出すと以下の様なコードを書いていました。

1000レコード弱を処理するのに33秒ほど。8000万行だと30日かかる計算に( д) ゚ ゚
これはイカンザキ

方法その2

コミットするのに時間がかかるということを風のうわさで聞いたので、コミットに時間がかかるのなら最後にコミットすれば良いじゃないというマリーアントワネット作戦

1000レコード弱を処理するのに8秒ほど。たっだ1行のインデントを変えただけで25%ほどの速度で処理できるように。
ただこれでも8000万レコードだと7.5日かかる計算。。

方法その3

どうやら世の中にはMultiple Insertというのがあるらしい。これで1000レコードごとに束ねてInsertすればさらに速度が向上するのでは
以下の例では2レコードのMultiple Insertですが、実際には1000レコード続けて書くようにして実施。

で、出てきたエラー(!)がこちら

なんでやねん!と思ったらAmazon RDSで借りられるデフォルトのDBの容量5GBでは足りていない模様。
ということで急遽100GBまで増額して再トライ。
1000レコードあたりの速度が0.18秒まで高速化!方法その2の3%ほどの速度で処理できるように。

ということで最終的に8000万レコードをぶっ込んだ際にかかった時間が以下

当初30日かかるという無謀な計画が、無事2時間で終えることができました。
最初どんだけ効率悪いコード書いてたんだって話ですねorzorz
Multiple Insert最強ね。

参考

インサート(insert)の処理方式別のパフォーマンスを検証 : 株式会社インターオフィス
http://www.inter-office.co.jp/contents/194/

MCMCサンプリング(非農業部門雇用者数編)

概要

・前回(http://nekopuni.holy.jp/?p=983)、非農業部門雇用者数(NFP)に対してローカルレベルモデルを適用し、そのパラメータは最尤推定により推定した。
・今回は、巷で話題沸騰のMCMCを用いて、パラメータの不確実性を考慮した推定を行なう。

最尤推定の問題点

最尤法でパラメータ \( \psi \)を推定できるからいいジャマイカというのが尤もな感覚なのですが、
『フィルタリング漸化式や平滑化漸化式にMLE\( \hat{ \psi } \)を使うのは一般的な慣習ではあるが、\( \psi \)の不確実性を適切に考慮しづらいという欠点がある。』(P.150)

最尤推定のアイディアとしては、現在取得できるデータに最もフィットするパラメータを選択するというものですから、真のパラメータとは外れたところで推定されている可能性もあり、そのパラメータのバラ付き具合も含めて推定してあげたほうが紳士的なんじゃないですかということですかね(ベイズ推定)

データ取得

前回同様(というかいつもですな)、NFPのデータはセントルイス連銀先生から。

All Employees: Total nonfarm (PAYEMS) – FRED – St. Louis Fed
http://research.stlouisfed.org/fred2/series/PAYEMS/

diffを取ったものを変数nfpに格納しておきます。

パラメータ(分散)の分布を仮定

ローカルレベルモデルも再掲

$$
y_t = \mu_t + v_t \\
\mu_t = \mu_{t-1} + w_t
$$
ここで\( v_t, w_t \)は正規分布に従い
$$
v_t \sim \left( 0, V \right) \\
w_t \sim \left( 0, W \right)
$$

今回推定する未知パラメータは\( V, W\)の2つです。
ここでパラメータ(分散)がガンマ分布に従うという仮定を置きます。
正確には、分散の逆数(以下、精度と言う)がガンマ分布に従うというものでして、

$$
V_t^{-1} = \phi_y \sim G \left( \alpha_y, \beta_y \right) \\
W_t^{-1} = \phi_\mu \sim G \left( \alpha_\mu, \beta_\mu \right)
$$

とします。

実際にMCMC推定を行なう

dlmパッケージを用いることで簡単にMCMCサンプリングを行なうことが出来ます。
中身でやっていることは、ギブスサンプリングでして、

\( \pi \left( \theta_{0:T} | y_{1:T}, \psi \right) \)から\( \theta_{0:T} \)
をサンプリングし、ここでサンプリングされた\( \theta_{0:T} \)を用いて
\( \pi \left( \psi | y_{1:T}, \theta_{0:T} \right) \)から \( \psi \)
をサンプリング

というのを何回も行なうというもの(らしい)です。片方を固定して、もう片方をサンプリングってイメージですかね。
で、どうやってコードで書くかというと

となります。
引数の意味としては
・a.y: \( \phi_y \)の事前平均
・b.y: \( \phi_y \)の事前分散
・a.theta: \( \phi_\mu \)の事前平均
・b.theta: \( \phi_\mu \)の事前分散
・n.sample: サンプリングする数
・thin: 何回おきにサンプルを保存するか
・save.states: 状態\( \mu \)を保存するか
となるようです。
事前平均と事前分布でなく、ガンマ分布の事前パラメータを指定することも可能で、その場合はshape.y, rate.y, shape.theta, rate.thetaとなるようです。ガンマ分布の事前平均・分散のほうが直感的に入れやすいということなんですかね。
thin = 1とすることで1個おきのサンプルのみ最終的に保存されます。すなわち、実際には44000回サンプリングは行なわれているということですね。

で、この関数は実行すると無茶苦茶時間がかかります。(Mac book proで1時間半ほど…)

また、このdlmGibbsDIG関数ですが、公式のドキュメントにもあるように、
『The d-inverse-gamma model is a consistent univariate DLM with unknown observation variance, diagonal system variance with unknown diagonal entries.』
となっておりますので、多変量DLMや、\( W \)が対角行列以外を想定した場合には対応していないということですかね。

結果

サンプルの生データをプロットしてみますと

140321_1

140321_2

といった具合でサンプルされたようです。
まあ一口にパラメータと言ってもこれだけバラつきますよというお話ですね。

ここで、最初の2000サンプルをバーンイン期間として除外し、パラメータの移動エルゴード平均を計算すると以下。
バーンイン期間を除いた上で、500個目から移動エルゴード平均を計算しています。次第にある程度の水準で収束してそうな気配がしておりますな。
(このバーンイン期間をどの程度取ればよいのかが謎ですが。)

140321_3

140321_4

また、バーンイン期間を除く、MCMCサンプルの平均をとると

となりまして、\( V \)の推定値が22042.9、\( W \)の推定値が4117.3となりました。
カッコ内の数字はモンテカルロ標準偏差の推定値。

最尤推定との比較

ちなみに最尤推定による結果は

となりました(前回記事投稿時から、NFPのデータが更新されているので推定しなおしております)。
まあこんなもんかなあという感じ

この手法による欠点

・dlmGibbsDIG関数を一般化できない
前述したように、多変量DLMには使えないですし、分散行列も対角行列でなければ対応してない。
(単変量DLMで分散行列に非対角要素が存在するケースを考えなければならない場合はそんなに無いかもしれないですが…)

・オンライン推定できない
新たな観測値が観測された場合、再度MCMCを行う必要がある。
非常に計算コストがかかるので、カルマンフィルタのように逐次的に計算できる方法のほうが良い
この問題点に関しては逐次モンテカルロ法により解決できる(らしい)。
ということで当面の課題としては、多変量の推定ができるようになるということとオンライン推定ができるようになることということで(― ―)(― ―)(― ―)

参考文献

Amazon.co.jp: Rによるベイジアン動的線形モデル (統計ライブラリー): G.ペトリス, S.ペトローネ, P.カンパニョーリ, 和合 肇, 萩原 淳一郎: 本
http://www.amazon.co.jp/dp/4254127960

Giovanni Petris: dlm: an R package for Bayesian analysis of Dynamic Linear Models
http://cran.r-project.org/web/packages/dlm/vignettes/dlm.pdf

Package ‘dlm’
http://cran.r-project.org/web/packages/dlm/dlm.pdf

メモ: Online Algorithms in High-Frequency Trading

Online Algorithms in High-frequency Trading – ACM Queue
http://queue.acm.org/detail.cfm?id=2534976

内容としては、HFTで用いられているone-pass algorithmsというものの紹介。
HFTで使われている一般的な方法を、私は全く知らないのでこの記事の内容が正しいのかどうかはわからないですが。

NYSEでは1秒間に215,000回もクオートが更新されるらしく、それを適切にインプットしてトレードのための意思決定を行なわなければならない。かつ、それらを高速に処理する必要がある。

そのためには以下の計算が主に重要だとか。
・時系列平均のオンライン推定
・ボラティリティのオンライン推定
・単回帰の回帰推定のオンライン推定

日次データを元に計算するというタイムスケールではなく、マイクロ秒ごと(かは知らんが)に計算する必要があるので、過去データをいちいち読みに行ってはメモリが足らんし計算も遅いぞなもしとのことで、再帰的に計算できるように通常の移動平均などの計算式を少し変えておりますな。
また、メモリ上に保存する値も最小限にすると。
再帰的に計算できるようにすることで、過去データの情報は全て1つの値に集約させているというイメージですな。

HFTが具体的にどうやって実装されているのか情報が中々手にはいらないので、たとえ古い情報であったとしてもなかなか参考になるお話ということでした。

まあ最近は
高頻度取引事業が表舞台に-米バーチュがIPO申請 – Bloomberg
http://www.bloomberg.co.jp/news/123-N2AQC16JTSF201.html

高頻度トレーディングをNY州が調査、不当な優位性が焦点 – Bloomberg
http://www.bloomberg.co.jp/news/123-N2MZ8J6K50YJ01.html

といったお話もあり時々ニュースにもなるHFT業界ということで。
FOMCステートメント発表の数秒前に米債だか為替が反応するとかいうこともあるとかないとかw

参考

Online Algorithms in High-frequency Trading – ACM Queue
http://queue.acm.org/detail.cfm?id=2534976

HFT における Online one-pass algorithm – Sparsegraph DotCom
http://i.ho.lc/blog/2014/03/online-onepass-algorithm-in-hft/

見せかけの回帰のまとめ

概要

・俺様まとめというヘッジ文言
・回帰分析(クロスセクション or 時系列)において発生する見せかけの回帰の対処法の備忘録
・ただ回帰を精緻にやったからといって儲かるかというとそれは別のうわ何をするやめ

見せかけの回帰の種類

1. クロスセクション回帰における見せかけの回帰
2. 時系列回帰における見せかけの回帰
2.1 定常過程の場合
2.2 非定常過程の場合

巷でよく言われるのは、2.2の非定常過程の場合ですが、この記事ではそれ以外のものも含めてまとめておきます。また、後から見返した時に(主に自分が)分からなくならないように参考文献のページも書いておりますので、文献をお持ちの方は適宜ご参照ということで何卒

1. クロスセクション回帰における見せかけの回帰

クロスセクションの回帰でも、本当は影響していないファクター同士の関係が有意に出ることがある。
良くある話としては、子どもの身長と足の速さに正の相関が見られるという類のもの。
言わずもがな、この背後には「年齢」という共通のファクターが存在する。

◆解決策:重回帰(上記の例の場合、年齢も説明変数に加える)を行なう。
これによって、年齢の影響を除外した上での、身長と足の速さに関する回帰係数を求めることができる。
(森崎 P.101)

2. 時系列回帰における見せかけの回帰

2.1 定常過程の場合

$$
y_t = a + b x_t + \epsilon_t
$$
というモデルを考える。上記を(1)式とする。

◆解決策:
(1)\(x \)を自己回帰モデル(AR過程)にぶち込み、定常かどうかを判断。
$$
x_t = \delta_x + \phi_1 x_{t-1} + \phi_2 x_{t-2} + \cdots + \phi_q x_{t-q} + v_{xt}
$$

定常かどうかは、Rではadf.test()で検定することができる(要tseriesパッケージ)。

帰無仮説:非定常過程
対立仮説:定常過程

(2)定常と判断されたら、(1)式のOLSを行ない、推定誤差\( \hat{\epsilon_t} \)を求める。
\( x \)が非定常と判断された場合は2.2に進む。

(3)推定誤差\( \hat{\epsilon_t} \)がiid(自己相関なし)ならば(2)のOLS推定量が最良漸近正規推定量(BAN推定量)となる。
自己相関があるかどうかはBreusch-Godfrey検定(森崎P.189)

ライブラリを読み込み、回帰式を指定してあげればその回帰式の推定誤差に関する自己相関を検定してくれる様子。
帰無仮説:自己相関なし
対立仮説:帰無仮説ではない
p値が大きかったら、系列相関があるとは言えないので、OLS推定量が一致推定量になりますと。
逆にp値が小さく、系列相関があると判断される場合は(4)に進む。

Rを使って計量経済分析 – Akihiko NODA
http://at-noda.com/econwiki/index.php?(URL長いので略)

(4)推定誤差\( \hat{\epsilon_t} \)が自己相関ありの場合、OLS推定量に一致性が得られない。
(5)yも自己回帰モデルにぶちこむ。
$$
y_t = \delta_y + \theta_1 x_{t-1} + \theta_2 x_{t-2} + \cdots + \theta_q x_{t-q} + v_{yt}
$$
ここで(1)と(5)で得られた\( v_{xt} \)と\( v_{yt} \)を用いてOLSを行なう。(森崎P.197から)
$$
v_{yt} = \beta v_{xt} + w_t
$$
このとき、\( \beta \)の推定量\( \hat{\beta} \)は(1)式\( b \)のBAN推定量となる。

2.2 非定常過程の場合

$$
y_t = a + b x_t + \epsilon_t
$$
という回帰において、\( x \)が非定常と判定されているとする。

◆解決策
(1)推定誤差\( \hat{\epsilon_t} \)が定常であるかどうかを検定
Engle-Granger共和文検定(沖本P.129、森崎P.247)
ここでADF検定を用いないことに注意。Engle-Granger検定に用いるインプットがそもそも推定値であるため、通常の単位根検定とは異なる棄却点を用いる必要がある。

\( \hat{\epsilon_t} \)が定常過程となった場合、\( x \)と\( y \)は共和分の関係となるので、この場合は例外的にOLS推定が一致推定量となるのでここで終了。(森崎P.242)

\( \hat{\epsilon_t} \)が非定常過程の場合、見せかけの回帰が発生する。
(1)ラグ変数の導入、または VARモデルによる推定を行なう
(2)\( x \)と\( y \)の階差を取ってOLS(沖本P.128)

メモ

OLSの推定誤差に関して、系列相関を検定したいときはBreusch-Godfrey検定、定常かどうかを検定したいときはEngle-Granger共和文検定となるというところが細かいところ

参考文献

・現代計量経済学: 森崎 初男: 本 – http://www.amazon.co.jp/dp/4904341031
・経済・ファイナンスデータの計量時系列分析 (統計ライブラリー): 沖本 竜義: 本 – http://www.amazon.co.jp/dp/4254127928

非農業部門雇用者数を状態空間モデルで予測する

概要

・状態空間モデルの練習も兼ねて非農業部門雇用者数(の変化)をモデル化し、1期先予測を行なう。
・使うモデルはローカルレベルモデル

ローカルレベルモデルとは

学んだばかりのものをドヤ顔で紹介するわけですが、ローカルレベルモデルとは以下のように定式化されます。
$$
y_t = \mu_t + v_t \\
\mu_t = \mu_{t-1} + w_t
$$
ここで\( v_t, w_t \)は正規分布に従い
$$
v_t \sim N \left( 0, V \right) \\
w_t \sim N \left( 0, W \right)
$$
と表されます。解釈としては、\( \mu_t \)が\( t \)時点での観測不可能な状態、これがランダムウォークに従い、また実際に観測される値(今回は非農業部門雇用者数(NFP)の変化が該当する)は、その状態から観測誤差( \( v_t \) )を伴って現れるということですね。

このローカルレベルモデルは、状態空間モデルの中の動的線形モデル(DLM: Dynamic Linear Models)というクラスに分類されます。状態がダイナミック(定数でない)ということですね。Rではdlmという素晴らしいライブラリがありますので、こちらを使用させていただきます。

データ取得

例の如く、セントルイス連銀先生からデータを取得します。

All Employees: Total nonfarm (PAYEMS) – FRED – St. Louis Fed
http://research.stlouisfed.org/fred2/series/PAYEMS/

diffを取り、plotすると以下のようになります(xtsオブジェクトに変換しています)。
dataが取得したデータを格納している変数です。

140301_1

NFP自体は1939年1月からの系列ですので、diffを取ると、1939年2月からの系列となります。
リーマン・ショック期や、第二次世界大戦中は大きく減少しています。

パラメータ推定

上述したDLMには正規分布の分散という未知パラメータ\( V, W \)がありますので、まずはこれを推定します。
参考文献(後述)のP.150あたりを参考にすると

dlmMLE関数で、初期値(0.1)、パラメータの下限(未知パラメータが分散なので負の値は取らない)を与えた上で最適化してくれます。

『DLMの尤度関数には、多数の局所最大が存在する可能性がある』(P.147)
『特に、要素convergenceは、常に確認する必要がある。この値が非零であれば、最小値への収束が達成されていない知らせとなる。』(P.148)

とありますので、先ほどの結果は、最小値は達成されているが、局所最大の可能性があるということですね。実際初期値を変えて推定し直すと、異なったパラメータが推定されます。が、そこまで大きく違いもないようなので、今回はこのまま使います。

カルマンフィルタ

先ほどの推定されたパラメータをもとにカルマンフィルタを動作させます。カルマンフィルタにより各時点での状態\( \mu_t \)を推定できます。

実際に推定結果をプロットすると以下

140301_2

黒線が実際の観測値、赤線がカルマンフィルタにより推定された状態\( \mu_t \)です。観測誤差分を除いた比較的滑らかな線が引かれています。これが観測不可能な状態の推定結果ということですね。

各時点の状態を推定することにより、1期先の観測値も予測することが出来ます。

1期先予測

2014年1月時点での値が152.277なので、2月のNFPの予測も152.277となりますね(正規分布の期待値は0なので)。
うん、何か当たらない気がする笑

どの程度予測力があるのか

何やら小難しいモデルを使っていますが、そもそもこのローカルレベルモデルは予測力あるんですかというのは当然の疑問。業界でよく使ってる3ヶ月平均とかでも良いんじゃねーのと思いますわな。
ということで予測力の比較を行います。

そこで使われるのがTheilのUというもので詳しくは参考文献P.100を参照。TheilのU自体は以下のようにして計算されます。
$$
U = \sqrt{ \frac{ \sum_{t=2}^n \left( y_t – f_t \right)^2 }{ \sum_{t=2}^n \left( y_t – y_{t-1} \right)^2 } }
$$

変化のないモデル(今期の値を来期の値の予測値とするモデル)でのMSEと、比較したいモデルのMSEの比をとっています。Uの値が小さいほど、当該モデルのMSEが小さいことを意味しますので、予測力も高いということですね。

まずはDLMの場合のTheilのU

カルマンフィルタがデータに適応する猶予を与えるため、最初の10サンプルは除外した上で計算しています。で、いくつかというと

ということで、変化のないモデルに対して、67%のMSEで予測を行なっていることが分かります。予測誤差が33%減少しているイメージですかね。

移動平均の計算自体は手前味噌ですが過去作った関数を流用します。

[R]テクニカル指標の計算
http://nekopuni.holy.jp/?p=974

3ヶ月移動平均と、6ヶ月移動平均も計算しておきます。

DLMのときと同じサンプルに対する予測力を比較するため、こちらも最初の10サンプルを除外しています。
で、いくつになるかというと

ということで、3ヶ月移動平均を用いた予測、6ヶ月移動平均を用いた予測とも、DLMよりは劣る予測力ということが分かります。とはいっても3ヶ月移動平均なんかは3%ほどの予測力の違いですので、コストパフォーマンスを考えたら移動平均でも良いんじゃないかという気がしないでもない結果に・・汗

観測値、DLM、移動平均のプロット

最後に、どんな感じでそれぞれのモデルが計算しているかを見ておきます。

140301_3

黒線が観測値、赤線がDLM、緑線が3ヶ月移動平均です。やはり2013年12月の74,000が効いていますね笑

参考文献

Amazon.co.jp: Rによるベイジアン動的線形モデル (統計ライブラリー): G.ペトリス, S.ペトローネ, P.カンパニョーリ, 和合 肇, 萩原 淳一郎: 本
http://www.amazon.co.jp/dp/4254127960