Sageで再現してみよう：集合知プログラミング第１０章株価
今回は、集合知プログラミング の第１０章にでてきます株価の因子分析を取り上げます。
あらかじめ読み込んでおくライブラリ
NMFの計算を行うnnmf.pyを読み込んでいます。
また、銘柄コードと銘柄名をプリントするために、_plist関数を定義します。
import urllib2
attach(DATA + 'nnmf.py')
attach(DATA + 'RUtil.py')
attach(DATA+'clusters.py')
|
|
# 日本語を含むリストを一行でプリントする
import json
def _plist(obj):
if isinstance(obj, list):
sys.stdout.write("[")
for i in range(len(obj)):
if i != 0: sys.stdout.write(", ")
orig = json.dumps(obj[i], indent=4)
sys.stdout.write(eval("u'''%s'''" % orig).encode('utf-8'))
sys.stdout.write("]\n")
|
|
銘柄一覧の作成
１０章では、Yahoo!ファイナンスを利用していましたが、日本の市場には対応していないみたいです。 そこで、フリーで株情報を公開化している株価データ ダウンロードサイト を利用することにします。
個別銘柄 株価時系列データから上位50個の銘柄を取り出して、銘柄コードリストtickersと銘柄辞書namesを作成します。
stocks = ["7751-T キヤノン", "8306-T 三菱UFJ", "7201-T 日産自", "3632-T グリー", "2432-T DENA",
"6753-T シャープ", "7203-T トヨタ", "9984-T ソフトバンク", "6502-T 東芝", "6857-T アドバンテ",
"8316-T 三井住友", "6301-T コマツ", "6752-T パナソニック", "6758-T ソニー", "5411-T JFEHD",
"8604-T 野村HD", "6954-T ファナック", "7267-T ホンダ", "6501-T 日立", "6762-T TDK",
"8411-T みずほ", "4503-T アステラス薬", "7974-O 任天堂", "4063-T 信越化", "8031-T 三井物",
"5401-T 新日鉄", "8058-T 三菱商", "9202-T ANA", "2914-T JT", "9432-T NTT",
"8802-T 菱地所", "8801-T 三井不", "9503-T 関西電", "4689-T ヤフー", "9437-T NTTドコモ",
"8035-T 東エレク", "8766-T 東京海上", "8830-T 住友不", "7731-T ニコン", "9983-T ファーストリテイ",
"4502-T 武田", "9433-T KDDI", "8591-T オリックス", "7733-T オリンパス", "6594-O 日電産",
"5214-T 日電硝", "3382-T 7&I-HD", "4452-T 花王", "7752-T リコー", "9831-T ヤマダ電"]
tickers = [s.split(' ')[0] for s in stocks]
names = {}
for s in stocks:
nt = s.split(' ')
names[nt[0]] = nt[1]
|
|
株価データの取得
urlの部分を修正して、株価データのcsvファイルをダウンロードし、日付、出来高を取り出し、 pricesとdatesにセットします。
shortest=250
prices = {}
dates = None
|
|
for t in tickers:
# urlをオープン
url = 'http://k-db.com/site/jikeiretsu.aspx?c=%s&year=0&download=csv'%t
rows = urllib2.urlopen(url).readlines()
# 各行の取引量のフィールドを抽出
prices[t] = [float(r.split(',')[5]) for r in rows[2:] if r.strip() != '']
if len(prices[t]) < shortest: shortest = len(prices[t])
# 日付をセット
if not dates:
dates = [r.split(',')[0] for r in rows[2:] if r.strip() != '']
|
|
因子分析
データが揃ったので、NMFを使って因子分析をします。
因子分析によって、銘柄の特徴ベクトルと日付の特徴ベクトルが求まります。
# 解析用データをセットアップ
l1 = [[ prices[tickers[i]][j] for i in range(len(tickers)) ] for j in range(shortest)]
|
|
w, h = factorize(matrix(l1), pc=10, iter=200)
6.05977578675e+18 3.20754906265e+17 2.49372897957e+17 1.93040816834e+17 1.51952713087e+17 1.28544850401e+17 1.14018781526e+17 1.03932334015e+17 9.61949585547e+16 8.97094612896e+16 8.38712056951e+16 7.86487234366e+16 7.43046936711e+16 7.06981681954e+16 6.76117710223e+16 6.49549803065e+16 6.26296608829e+16 6.06215150688e+16 5.89213120721e+16 5.74971625714e+16 6.05977578675e+18 3.20754906265e+17 2.49372897957e+17 1.93040816834e+17 1.51952713087e+17 1.28544850401e+17 1.14018781526e+17 1.03932334015e+17 9.61949585547e+16 8.97094612896e+16 8.38712056951e+16 7.86487234366e+16 7.43046936711e+16 7.06981681954e+16 6.76117710223e+16 6.49549803065e+16 6.26296608829e+16 6.06215150688e+16 5.89213120721e+16 5.74971625714e+16 |
# print h
# print w
|
|
結果の表示
特徴ベクトルの値大きな順にソートし、上位６個の銘柄と日付の上位３個の銘柄を表示します。
特徴ベクトルの上位銘柄の内、トップの値が他の値と比べて突出している場合には、 トップの企業のスキャンダルや特別なイベントによって出来高が大きく変動したことを表し、 逆に値が揃っている場合には、株価の出来高に関するグループの候補になると考えられます。
Feature 8を例にとると、2011-11-10にオリンパスの粉飾決算で「オリンパス」が管理銘柄になり、 その翌日の2011-11-11に大量の売りがでたことに起因していると思われます。
株価の情報については、 日々の日経平均株価を記録するブログを参考にさせて頂きました。 このブログは、日々の株価の情報をきちんと記録してあるので、今回のような結果の検証にはとても助かりました。
# Loop over all the features
for i in range(shape(h)[0]):
print "Feature %d" %i
# Get the top stocks for this feature
ol=[(h[i,j],names[tickers[j]]) for j in range(shape(h)[1])]
ol.sort()
ol.reverse()
for j in range(6):
_plist(list(ol[j]))
print
# Show the top dates for this feature
porder=[(w[d,i],d) for d in range(shortest)]
porder.sort()
porder.reverse()
print [(p[0],dates[p[1]]) for p in porder[0:3]]
print
Feature 0 [32620663.395954657, "日立"] [9487259.9985786844, "みずほ"] [5704258.1124020452, "三菱UFJ"] [5664853.7776350277, "日産自"] [5376607.0328148548, "東芝"] [3310128.5223182668, "新日鉄"] [(2.5192355540360323, '2011-08-04'), (1.6023115788708648, '2012-02-03'), (1.5610543992850998, '2012-04-19')] Feature 1 [30490097.004402213, "野村HD"] [15525867.294821035, "三菱UFJ"] [9054191.8785016369, "新日鉄"] [6282340.5948974723, "東芝"] [5674870.6713392464, "パナソニック"] [3528139.1283912286, "シャープ"] [(2.3935164256115891, '2012-03-02'), (2.3285530583857188, '2011-11-08'), (2.1979519659808759, '2011-11-09')] Feature 2 [47684429.391790703, "みずほ"] [26341818.955080986, "三菱UFJ"] [2467344.1345816916, "日立"] [2030294.1606793641, "パナソニック"] [1967484.947207479, "新日鉄"] [1890776.8176482248, "野村HD"] [(2.600235535777176, '2012-03-09'), (2.4691405420226942, '2011-08-09'), (2.4275550869614042, '2011-10-27')] Feature 3 [28045026.598636635, "シャープ"] [16269443.135748301, "みずほ"] [7624325.7699214518, "三菱UFJ"] [5790908.8732109191, "パナソニック"] [3537168.0596196144, "日立"] [3260164.1475366801, "ソニー"] [(3.9224885229722863, '2012-03-29'), (3.06718640263529, '2012-03-14'), (2.8986486073852111, '2012-02-03')] Feature 4 [21822077.011053085, "東芝"] [18896045.585721053, "三菱UFJ"] [12631138.268594489, "日産自"] [8486027.9185644463, "みずほ"] [5151400.829650715, "新日鉄"] [4072078.5143724168, "グリー"] [(2.086093388601614, '2011-08-24'), (2.0583185990393442, '2011-08-25'), (1.9529996931095019, '2012-07-25')] Feature 5 [18314236.598199803, "オリンパス"] [5893444.9026885722, "東芝"] [4294919.9610484838, "日立"] [3969107.318491233, "野村HD"] [3253220.7821490159, "新日鉄"] [2318741.2907765107, "三井物"] [(4.6085283937260808, '2011-11-11'), (4.0564409840004565, '2011-10-18'), (3.8584932529132638, '2011-10-27')] Feature 6 [22447161.82133168, "ANA"] [9385166.5665901937, "三菱UFJ"] [6597214.6294900449, "みずほ"] [4099716.4186763056, "シャープ"] [3809638.9915798255, "パナソニック"] [3776402.5976034272, "新日鉄"] [(4.4852823808006139, '2012-07-03'), (3.4359807031676879, '2012-07-04'), (3.330666172522645, '2012-07-26')] Feature 7 [41117978.470087767, "みずほ"] [13358547.630855544, "東芝"] [3153311.0035421075, "日立"] [2096301.6571017692, "グリー"] [2067527.1205180509, "日産自"] [1628227.7816972502, "トヨタ"] [(3.7874413297685976, '2011-08-05'), (3.0447265440916613, '2011-08-09'), (2.5914566809561776, '2011-09-06')] Feature 8 [23723279.898510188, "みずほ"] [17132214.508064192, "新日鉄"] [9968688.5243205298, "三菱UFJ"] [3791076.7238644976, "日立"] [2445764.2275681193, "日産自"] [1792965.3679678829, "三井物"] [(2.9054428270557868, '2011-07-27'), (2.6921681041609165, '2011-10-28'), (2.298575164748399, '2012-02-08')] Feature 9 [34339148.518358395, "みずほ"] [28142120.391353182, "野村HD"] [4048380.7955130935, "三菱UFJ"] [1950625.7868780512, "三井住友"] [1413048.3704661257, "日立"] [1340900.773408686, "三菱商"] [(4.0108105767297788, '2011-11-08'), (3.5348559499404941, '2011-11-09'), (2.3390343512254157, '2012-03-22')] Feature 0 [32620663.395954657, "日立"] [9487259.9985786844, "みずほ"] [5704258.1124020452, "三菱UFJ"] [5664853.7776350277, "日産自"] [5376607.0328148548, "東芝"] [3310128.5223182668, "新日鉄"] [(2.5192355540360323, '2011-08-04'), (1.6023115788708648, '2012-02-03'), (1.5610543992850998, '2012-04-19')] Feature 1 [30490097.004402213, "野村HD"] [15525867.294821035, "三菱UFJ"] [9054191.8785016369, "新日鉄"] [6282340.5948974723, "東芝"] [5674870.6713392464, "パナソニック"] [3528139.1283912286, "シャープ"] [(2.3935164256115891, '2012-03-02'), (2.3285530583857188, '2011-11-08'), (2.1979519659808759, '2011-11-09')] Feature 2 [47684429.391790703, "みずほ"] [26341818.955080986, "三菱UFJ"] [2467344.1345816916, "日立"] [2030294.1606793641, "パナソニック"] [1967484.947207479, "新日鉄"] [1890776.8176482248, "野村HD"] [(2.600235535777176, '2012-03-09'), (2.4691405420226942, '2011-08-09'), (2.4275550869614042, '2011-10-27')] Feature 3 [28045026.598636635, "シャープ"] [16269443.135748301, "みずほ"] [7624325.7699214518, "三菱UFJ"] [5790908.8732109191, "パナソニック"] [3537168.0596196144, "日立"] [3260164.1475366801, "ソニー"] [(3.9224885229722863, '2012-03-29'), (3.06718640263529, '2012-03-14'), (2.8986486073852111, '2012-02-03')] Feature 4 [21822077.011053085, "東芝"] [18896045.585721053, "三菱UFJ"] [12631138.268594489, "日産自"] [8486027.9185644463, "みずほ"] [5151400.829650715, "新日鉄"] [4072078.5143724168, "グリー"] [(2.086093388601614, '2011-08-24'), (2.0583185990393442, '2011-08-25'), (1.9529996931095019, '2012-07-25')] Feature 5 [18314236.598199803, "オリンパス"] [5893444.9026885722, "東芝"] [4294919.9610484838, "日立"] [3969107.318491233, "野村HD"] [3253220.7821490159, "新日鉄"] [2318741.2907765107, "三井物"] [(4.6085283937260808, '2011-11-11'), (4.0564409840004565, '2011-10-18'), (3.8584932529132638, '2011-10-27')] Feature 6 [22447161.82133168, "ANA"] [9385166.5665901937, "三菱UFJ"] [6597214.6294900449, "みずほ"] [4099716.4186763056, "シャープ"] [3809638.9915798255, "パナソニック"] [3776402.5976034272, "新日鉄"] [(4.4852823808006139, '2012-07-03'), (3.4359807031676879, '2012-07-04'), (3.330666172522645, '2012-07-26')] Feature 7 [41117978.470087767, "みずほ"] [13358547.630855544, "東芝"] [3153311.0035421075, "日立"] [2096301.6571017692, "グリー"] [2067527.1205180509, "日産自"] [1628227.7816972502, "トヨタ"] [(3.7874413297685976, '2011-08-05'), (3.0447265440916613, '2011-08-09'), (2.5914566809561776, '2011-09-06')] Feature 8 [23723279.898510188, "みずほ"] [17132214.508064192, "新日鉄"] [9968688.5243205298, "三菱UFJ"] [3791076.7238644976, "日立"] [2445764.2275681193, "日産自"] [1792965.3679678829, "三井物"] [(2.9054428270557868, '2011-07-27'), (2.6921681041609165, '2011-10-28'), (2.298575164748399, '2012-02-08')] Feature 9 [34339148.518358395, "みずほ"] [28142120.391353182, "野村HD"] [4048380.7955130935, "三菱UFJ"] [1950625.7868780512, "三井住友"] [1413048.3704661257, "日立"] [1340900.773408686, "三菱商"] [(4.0108105767297788, '2011-11-08'), (3.5348559499404941, '2011-11-09'), (2.3390343512254157, '2012-03-22')] |
株価チャートの表示
ここで、オリンパスの株価チャートがどのように変換しているか、RのRFinanceYJパッケージを使って表示してみましょう。
RFinanceYJ、quantmodがインストールされていない場合には、コメントを外して実行してください。 RFinanceYJのインストールには、libxml2-devが必要となります。
#r('install.packages("RFinanceYJ")')
#r('install.packages("quantmod")')
r('library(RFinanceYJ)')
r('library(quantmod)')
r('library(xts)')
[1] "quantmod" "TTR" "Defaults" "RFinanceYJ" "xts" "zoo" "XML" [8] "stats" "graphics" "grDevices" "utils" "datasets" "methods" "base" [1] "quantmod" "TTR" "Defaults" "RFinanceYJ" "xts" "zoo" "XML" [8] "stats" "graphics" "grDevices" "utils" "datasets" "methods" "base" |
graph = preGraph('olympus.pdf')
r('olympus <- quoteStockTsData("7733.t", "2011-11-01")')
r('names(olympus)<-c("Date","Open","High","Low","Close","Volume")')
r('olympus<-read.zoo(olympus,tz="")')
r('candleChart(olympus)')
postGraph(graph)
|
|
特徴的な銘柄の比較
NMF分析で抽出した特徴的な銘柄の株価と出来高をプロットしてみます。 ここでは、野村HD、三菱UFJ、グリー、DENAを例とします。
# Feature 0の確認
graph = preGraph('Nomura.pdf')
r('Nomura <- quoteStockTsData("8604.t", "2011-11-01")')
r('names(Nomura)<-c("Date","Open","High","Low","Close","Volume")')
r('Nomura<-read.zoo(Nomura,tz="")')
r('candleChart(Nomura)')
offGraph()
graph = preGraph('UFJ.pdf')
r('UFJ <- quoteStockTsData("8306.t", "2011-11-01")')
r('names(UFJ)<-c("Date","Open","High","Low","Close","Volume")')
r('UFJ<-read.zoo(UFJ,tz="")')
r('candleChart(UFJ)')
offGraph()
graph = preGraph('DENA.pdf')
r('DENA <- quoteStockTsData("2432.t", "2011-11-01")')
r('names(DENA)<-c("Date","Open","High","Low","Close","Volume")')
r('DENA<-read.zoo(DENA,tz="")')
r('candleChart(DENA)')
offGraph()
graph = preGraph('Gree.pdf')
r('Gree <- quoteStockTsData("3632.t", "2011-11-01")')
r('names(Gree)<-c("Date","Open","High","Low","Close","Volume")')
r('Gree<-read.zoo(Gree,tz="")')
r('candleChart(Gree)')
offGraph()
|
|
html.table([[getGraph('Nomura.pdf', 0.5), getGraph('UFJ.pdf', 0.5)],
[getGraph('DENA.pdf', 0.5), getGraph('Gree.pdf', 0.5)]])
|
|
銘柄の特徴ベクトルによるクラスタリング
最後に、銘柄の特徴ベクトル（計算では$h^T$）を使って株の銘柄の クラスタリングを行います。
グリーとDENA、三井物と三菱商、パナソニックとソニー等きれいに分類できています。
株の出来高をNMFで因子分析するだけで、銘柄のクラスタリングや関連づけができることに 正直驚きました。
# デンドログラムで日本語を表示するための修正
from PIL import Image,ImageDraw,ImageFont
font = ImageFont.truetype('/usr/local/share/Fonts/MS Mincho.ttf', 14)
|
|
# データのセットアップ
data = h.transpose()
(ii, jj) = shape(data)
v = [[data[i, j] for j in range(jj)] for i in range(ii)]
namelist = [names[t] for t in tickers]
# 株価のクラスタリング
clust=hcluster(v)
drawdendrogram(clust,namelist,png=DATA + 'stockclust.png')
showPNG('stockclust.png', fac=1)
|
|
|
|