[[FrontPage]]
#contents
2011/07/27からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www15191ue.sakura.ne.jp:8000/home/pub/20/
また、Sageのサーバを公開しているサイト（http://www.sagenb.org/, http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを
アップロードし、実行したり、変更していろいろ動きを試すことができます。
* Sageで再現してみよう：集合知プログラミング第１０章株価 [#m8f6a96b]
今回は、
[[集合知プログラミング>http://www.amazon.co.jp/dp/4873113644/]]
の第１０章にでてきます株価の因子分析を取り上げます。
*** あらかじめ読み込んでおくライブラリ [#t07157db]
NMFの計算を行うnnmf.pyを読み込んでいます。
また、銘柄コードと銘柄名をプリントするために、_plist関数を定義します。
sageへの入力:
#pre{{
import urllib2
attach(DATA + 'nnmf.py')
}}
sageへの入力:
#pre{{
# 日本語を含むリストを一行でプリントする
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")
}}
** 銘柄一覧の作成 [#p04a14f9]
１０章では、Yahoo!ファイナンスを利用していましたが、日本の市場には対応していないみたいです。
そこで、フリーで株情報を公開化している
[[株価データ ダウンロードサイト>http://k-db.com/]]
を利用することにします。
個別銘柄 株価時系列データから上位50個の銘柄を取り出して、銘柄コードリストtickersと銘柄辞書namesを作成します。
sageへの入力:
#pre{{
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]
}}
** 株価データの取得 [#i66117e5]
urlの部分を修正して、株価データのcsvファイルをダウンロードし、日付、出来高を取り出し、
pricesとdatesにセットします。
sageへの入力:
#pre{{
shortest=250
prices = {}
dates = None
}}
sageへの入力:
#pre{{
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() != '']
}}
** 因子分析 [#m3637072]
データが揃ったので、NMFを使って因子分析をします。
因子分析によって、銘柄の特徴ベクトルと日付の特徴ベクトルが求まります。
sageへの入力:
#pre{{
# 解析用データをセットアップ
l1 = [[ prices[tickers[i]][j] for i in range(len(tickers)) ] for j in range(shortest)]
}}
sageへの入力:
#pre{{
w, h = factorize(matrix(l1), pc=10, iter=100)
}}
sageの出力:
#pre{{
6.04500268713e+18
2.80667909388e+17
2.21789192007e+17
1.6802439138e+17
1.34101277419e+17
1.17386562415e+17
1.07283238074e+17
9.89422053839e+16
9.17226416359e+16
8.54569035958e+16
}}
sageへの入力:
#pre{{
# print h
# print w
}}
*** 結果の表示 [#h59e49c1]
特徴ベクトルの値大きな順にソートし、上位６個の銘柄と日付の上位３個の銘柄を表示します。
特徴ベクトルの上位銘柄の内、トップの値が他の値と比べて突出している場合には、
トップの企業のスキャンダルや特別なイベントによって出来高が大きく変動したことを表し、
逆に値が揃っている場合には、株価の出来高に関するグループの候補になると考えられます。
Feature 8を例にとると、2011-11-10にオリンパスの粉飾決算で「オリンパス」が管理銘柄になり、
その翌日の2011-11-11に大量の売りがでたことに起因していると思われます。
株価の情報については、
[[日々の日経平均株価を記録するブログ>http://hajimetekabu.at.webry.info/201111/article_7.html]]
日々の日経平均株価を記録するブログ</a>を参考にさせて頂きました。
このブログは、日々の株価の情報をきちんと記録してあるので、今回のような結果の検証にはとても助かりました。
sageへの入力:
#pre{{
# 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
}}
sageの出力:
#pre{{
Feature 0
[32195680.011513639, "野村HD"]
[3086179.4975982769, "三菱UFJ"]
[1930154.0709872572, "シャープ"]
[1524676.4505012867, "DENA"]
[1441506.620904193, "グリー"]
[1403828.3563206003, "みずほ"]
[(5.566644794683608, '2011-11-08'), (4.3412876554898618, '2011-11-09'), (2.7223269465277533, '2012-03-22')]
Feature 1
[39921451.144022748, "みずほ"]
[12339911.938650882, "日立"]
[2356448.2864209949, "日産自"]
[1595364.5536408355, "三井物"]
[1574978.9224391158, "三菱商"]
[1354693.8981276902, "新日鉄"]
[(4.3604560180911491, '2011-08-05'), (2.6085299776570623, '2011-10-27'), (2.6028653721834192, '2011-08-04')]
Feature 2
[35987745.989190437, "みずほ"]
[12519388.346624229, "三菱UFJ"]
[9412745.1724054832, "野村HD"]
[6480996.3462522496, "新日鉄"]
[5315060.9037663853, "東芝"]
[3242649.8347622748, "日立"]
[(2.8311815125392772, '2012-02-29'), (2.3629085832450984, '2012-01-20'), (2.2801930512296371, '2012-03-09')]
Feature 3
[36219908.096974887, "みずほ"]
[10593280.474430121, "野村HD"]
[9768064.1926213428, "新日鉄"]
[7475104.368119875, "東芝"]
[4195633.6495521842, "日産自"]
[3089235.5608142023, "三菱UFJ"]
[(2.3764517495483597, '2012-02-24'), (1.8743453892710862, '2012-02-17'), (1.8222639913130423, '2012-03-09')]
Feature 4
[25729092.895349674, "東芝"]
[22357497.833596095, "みずほ"]
[10084822.004651317, "三菱UFJ"]
[9453710.0034469981, "日立"]
[6572206.9072962748, "日産自"]
[6304374.0185764916, "野村HD"]
[(2.6277483882989374, '2011-09-06'), (2.290737335365054, '2011-08-09'), (1.9473346150181128, '2012-07-25')]
Feature 5
[13940864.544791378, "日立"]
[10684520.210738182, "三菱UFJ"]
[10286645.924499465, "新日鉄"]
[7827175.5684815394, "東芝"]
[6787514.2662564144, "日産自"]
[3289340.8039039453, "パナソニック"]
[(3.0882025657374186, '2011-08-04'), (3.0102643685283308, '2011-07-28'), (2.9149104140323363, '2012-05-18')]
Feature 6
[26695633.969163321, "シャープ"]
[8891764.6984640788, "三菱UFJ"]
[7694520.3495533559, "パナソニック"]
[7381981.3161944989, "みずほ"]
[6042713.222232963, "東芝"]
[4567403.6055205511, "ANA"]
[(3.7791941438265502, '2012-03-29'), (3.0578822320743528, '2012-07-24'), (2.9772394275265657, '2012-03-14')]
Feature 7
[30016063.184179123, "みずほ"]
[27681098.353298899, "三菱UFJ"]
[8667570.7121750489, "新日鉄"]
[5119809.7323025586, "野村HD"]
[3940730.7801396595, "日産自"]
[3134288.4979333193, "パナソニック"]
[(2.9815860916648962, '2012-02-15'), (2.1606569948701679, '2012-02-08'), (2.1340813374355472, '2012-02-24')]
Feature 8
[19117288.382176194, "オリンパス"]
[10367649.103981823, "みずほ"]
[7050057.8352865539, "三菱UFJ"]
[6380378.640513096, "東芝"]
[4699999.5834081303, "新日鉄"]
[2328053.7742331196, "三井物"]
[(4.3108477138906007, '2011-11-11'), (3.8831084004940113, '2011-10-27'), (3.4838826490980135, '2011-10-18')]
Feature 9
[23462047.317239508, "みずほ"]
[19277882.581798293, "三菱UFJ"]
[15789711.512183005, "ANA"]
[6873154.6707015941, "野村HD"]
[3049868.8044577157, "日立"]
[2349840.746904158, "シャープ"]
[(4.6316689004078482, '2012-07-03'), (2.9613342286902049, '2012-07-04'), (2.1278412401367111, '2012-07-26')]
}}
** 株価チャートの表示 [#ya6d6ca6]
ここで、オリンパスの株価チャートがどのように変換しているか、RのRFinanceYJパッケージを使って表示してみましょう。
RFinanceYJ、quantmodがインストールされていない場合には、コメントを外して実行してください。
RFinanceYJのインストールには、libxml2-devが必要となります。
sageへの入力:
#pre{{
#r('install.packages("RFinanceYJ")')
#r('install.packages("quantmod")')
r('library(RFinanceYJ)')
r('library(quantmod)')
r('library(xts)')
}}
sageの出力:
#pre{{
[1] "quantmod" "TTR" "Defaults" "RFinanceYJ" "xts" "zoo" "XML"
[8] "stats" "graphics" "grDevices" "utils" "datasets" "methods" "base"
}}
sageへの入力:
#pre{{
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)
}}
&ref(out1.png);
*** 特徴的な銘柄の比較 [#m615c5e1]
NMF分析で抽出した特徴的な銘柄の株価と出来高をプロットしてみます。 ここでは、野村HD、三菱UFJ、グリー、DENAを例とします。
sageへの入力:
#pre{{
# 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()
}}
sageへの入力:
#pre{{
html.table([[getGraph('Nomura.pdf', 0.5), getGraph('UFJ.pdf', 0.5)],
[getGraph('DENA.pdf', 0.5), getGraph('Gree.pdf', 0.5)]])
}}
&ref(out2.png);
** 銘柄の特徴ベクトルによるクラスタリング [#va35a82a]
最後に、銘柄の特徴ベクトル（計算では\(h^T\)）を使って株の銘柄の クラスタリングを行います。
グリーとDENA、三井物と三菱商、パナソニックとソニー等きれいに分類できています。
株の出来高をNMFで因子分析するだけで、銘柄のクラスタリングや関連づけができることに 正直驚きました。
sageへの入力:
#pre{{
# デンドログラムで日本語を表示するための修正
from PIL import Image,ImageDraw,ImageFont
font = ImageFont.truetype('/usr/local/share/Fonts/MS Mincho.ttf', 14)
}}
sageへの入力:
#pre{{
# データのセットアップ
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)
}}
&ref(out3.png);
** コメント [#c223c74e]
#vote(おもしろかった[3],そうでもない[0],わかりずらい[0])
#vote(おもしろかった[4],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
