日々のメモ

統計解析や機械学習、ファイナンスについての勉強メモ

Rによる機械学習 最尤推定法その1

機械学習 統計

よくある例題を使って最尤推定法をやってみます

> training
 [1]  0.02841715  0.34097049  1.37103609  0.52697056  0.88287394  0.20673168
 [7] -0.92268387 -1.41716486 -0.94624667 -0.56053292  0.21504500
> test
 [1]  0.2744391  0.2354888  1.0005627  0.7502364  0.5542613  0.1038660
 [7] -0.7504576 -0.2694111 -0.9737203 -0.5597848  0.2300776

f:id:Kenjiqun:20160622203105p:plain

グラフはtrainingのデータをプロットしています

ある関数f(x)にノイズεが含まれていますことがわかっています
ノイズεは確率変数で何らかの確率分布に従っています

y=f(x)+ε

このデータからf(x)を推定することを考えます
こういった問題は最小二乗法や最尤推定法が用いられます

今回は最尤推定法を用いてf(x)を推定してみようと思います。
最尤推定法は以下のステップからなります

⓪ノイズの確率分布を仮定する。確率分布を定めるパラメータを導入する(平均、分散など)
①関数形を仮定する。関数を定めるパラメータを導入する(例えばg(x)=axのaがパラメータ)
②観測されたデータが出現する確率Pを計算する
③Pが最大になるようなパラメータを決める

最尤推定法は

「観測されたデータが最も出現確率が高いデータのはずだ」

という仮説のもとでパラメータを決定します

これは例えば表が確率q、裏が確率1-qででるコイン(pは未知)を5回投げて
表、表、表、裏、表
と出た時にq=4/5と推定するのが最も合理的だ、というのと同じ考え方です

ちなみにコインの例で最尤推定法を行うと
観測されたデータが出現する確率PはP=(1-q)q^4となってこれを最大化するようなpは

∂P/∂q=0よりq=4/5となります



TeXかなんかで数式を入れたかったのですが、何が一番良いのか模索中なので続きはまた次回

Rによる統計解析 日経平均先物その2

統計

日経平均先物の分布について確認してみたいと思います
Black-Scholesモデルなど、株価のリターンの分布は対数正規分布と仮定されることが多いですが、実際はどのような分布になっているか調べてみます

今回は2010/6/21〜2016/6/21の日足データを使います
Macであれば楽天証券のMarketSpeedからデータをcsvで保存します
f:id:Kenjiqun:20160621103244p:plain
この中で、

> file <- file("NikkeiSakimono.csv", open="r", encoding="cp932")
> data <- read.table(file, header=TRUE, sep=",")

ですがこのままだとデータがfactor型になってしまいます
factor型のままでは解析できません

> sapply(data,class)
    日付     始値     高値     安値     終値   前日比   出来高    X5DMA
"factor" "factor" "factor" "factor" "factor" "factor" "factor" "factor"
  X25DMA   X75DMA   X5DVMA  X25DVMA  X75DVMA
"factor" "factor" "factor" "factor" "factor"

桁区切りを無くせばread.tableした時にnumeric型になってくれます
ということでExcelで桁区切りを一旦なくしてから再度読み込みました


あと日付のデータ型を日付型に変換します

> data$日付<-as.Date(data$日付)

f:id:Kenjiqun:20160621105210p:plain

株価のリターンとしては前日比を用いることにして、ヒストグラムをかきます

リターンは収益率にするべきでした
リターン=前日比/前日の終値
とします

> return <- data$前日比[-1]/data$終値[1:length(data$終値)-1]

[-1]で第1成分以外を指定し、[1:length(data$終値)-1]で最後の成分以外を指定します
リターンのヒストグラムは

> hist(return,main="Histogram of Return",xlab="return",freq=FALSE,breaks=50)

f:id:Kenjiqun:20160621215630p:plain

freq=FALSEとすることで縦軸を頻度ではなく確率密度にしています
breaksは分割数を指定しています

さて、パッと見正規分布に見えますがどうでしょう

> mu=mean(data$前日比)
> sigma=sd(data$前日比)
> curve(dnorm(x,mu,sigma),add=T)

dnorm(x,mu,sigma)で平均mu分散sigmaの正規分布を作ってヒストグラムに重ねてみます

f:id:Kenjiqun:20160621220143p:plain

これを見ると日経平均先物は正規分布と比べて端の方に分布が偏っているように見えます
実際、分散は0.0145くらいなので6σ程度まで分布が存在することになります

6σといえば1億分の1以下なので正規分布だとまあほぼほぼ出現しないはずですが、実際にはその辺まで現れます
暴騰暴落が起こりやすいということでその辺に気をつけないと死んでしまいそうです

リターンの分布がべき分布になるとか、フラクタル性を持つとかいう話もありますがそれはまた今度にしようと思います

Rによる統計解析 日経平均先物その1

統計

Rで統計解析を行ってみます(参考文献 : Rによる統計解析)

とは言ってもHello World!から行きましょう

まずはRのインストール(Macでやってます)
The Comprehensive R Archive Network

インストールしたら

$ R

R version 3.3.0 (2016-05-03) -- "Supposedly Educational"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin13.4.0 (64-bit)

R は、自由なソフトウェアであり、「完全に無保証」です。
一定の条件に従えば、自由にこれを再配布することができます。
配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。

R は多くの貢献者による共同プロジェクトです。
詳しくは 'contributors()' と入力してください。
また、R や R のパッケージを出版物で引用する際の形式については
'citation()' と入力してください。

'demo()' と入力すればデモをみることができます。
'help()' とすればオンラインヘルプが出ます。
'help.start()' で HTML ブラウザによるヘルプがみられます。
'q()' と入力すれば R を終了します。

>

で起動できる(簡単)

では早速

> cat('Hello World!\n')
Hello World!

無事Hello World!できました



具体的に何かのテーマを決めてやった方がわかりやすそうなので
日経平均先物のデータを使ってみたいと思います

データはこちらから
http://k-db.com/futures/

今回は2016年09月限の6/17の5分足データを扱ってみます
f:id:Kenjiqun:20160620153629p:plain

こいつをcsvでダウンロードしてきます。
このデータに関してですが、とりあえず6/17に取引された値段が5分ごとに記録されたものです。
Rにcsvファイルを読み込ませるにはread.tableを用います

> ( data <- read.table("futures_F101-1609_5m_20160617.csv", header=TRUE, sep=",") )

ここで
headerは1行目の"日付 時刻 始値 高値 安値 終値 出来高 売買代金"をヘッダとして扱うことを意味します
またカンマ区切りなので sep="," としてやります
Rはデータフレームというリスト構造を持っており、read.tableによってデータフレームdataにcsv内のデータを構造を保ったまま保存することができます
ここ大事

Macでやると

>>  ( data <- read.table("futures_F101-1609_5m_20160617.csv", header=TRUE,sep=",") )
 make.names(col.names, unique = TRUE) でエラー:
   '<93><fa><95>t' に不正なマルチバイト文字があります

とエラーになってしまいます
Rが想定している文字コードがMacではutf-8なんですがhttp://k-db.com/futures/のファイルの文字コードがcp932であることに起因します
そのため、encodeしてやる必要があります

> filename <- file("futures_F101-1609_5m_20160617.csv", open="r", encoding="cp932")
> ( data <- read.table(filename, header=TRUE, sep=",") )
          日付  時刻  始値  高値  安値  終値 出来高 売買代金
1   2016-06-18 03:00 15580 15580 15580 15580     82  1277560
...

これで準備ができました。
日付や高値などは変数として扱うことができます
データフレーム名$変数名
データフレーム名 [, 列番号]
データフレーム名 [列番号]
によって変数にアクセスすることができます
例えば

> data$安値
  [1] 15580    NA 15560 15560 15550 15540 15540 15540 15540 15540 15530 15500
 [13] 15500 15500 15490 15490 15510 15510 15510 15500 15500 15500 15500 15480
 [25] 15470 15450 15460 15450 15440 15430 15410 15420 15440 15430 15420 15440
 [37] 15440 15470 15460 15470 15460 15440 15440 15460 15490 15510 15500 15490
 [49] 15480 15490 15500 15490 15480 15480 15480 15510 15500 15500 15500 15500
 [61] 15510 15520 15520 15510 15510 15520 15520 15530 15540 15540 15530 15530
 [73] 15530 15520 15510 15510 15520 15510 15500 15510 15500 15500 15520 15530
 [85] 15530 15530 15530 15530 15530 15540 15550    NA 15540 15530 15520 15530
 [97] 15530 15530 15550 15560 15560 15560 15550 15540 15530 15540 15540 15530
[109] 15530 15530 15510 15510 15520 15510 15490 15480 15470 15480 15460 15480
[121] 15480 15490 15540 15540 15530 15530 15540 15530    NA 15570 15540 15560
[133] 15550 15550 15570 15570 15550 15540 15550 15560 15550 15530 15540 15560
[145] 15570 15580 15580 15580 15580 15570 15560 15550 15560 15560 15580 15580
[157] 15610 15620 15610 15630 15640 15680 15680 15680 15670 15640 15640 15640
[169] 15640 15640 15630 15630 15640 15630 15620 15620 15620 15590 15600 15590
[181] 15560 15640 15630 15610 15610 15640 15650 15680 15700 15700 15690 15680
[193] 15660 15620 15610 15630 15630 15630 15620 15590 15620 15630 15610

いくつか NA となっていますが、これはデータの欠損を示しています
日経平均先物では取引が行われない時間帯があります
www.rakuten-sec.co.jp
実際、取引が行われない時間帯のデータが欠損していることがわかります

今は問題ないですが、実際に統計値を求めるときにはNAの扱いに注意が必要になります

次に読み込んだデータをプロットしてみようと思います
日付と時刻データが分離しているため、単純にプロットすると分かりづらくなりそうです(例えば横軸を時刻にすると日付の順がおかしくなるかもしれない)
そこで日付と時刻をくっつけてみます
Rでは日付型(年月日)とPOSIXct,POSIXlt,POSIXt(日付時分秒)があります
POSIXctは1970年元旦からの経過秒で、POSIXltは可読性が高くなるように経過秒を日付時分秒の文字列リストにしたものです
POSIXtはその両方を含んでいます
今回はpaste関数で日付と時刻を結合し、as.POSIXlt関数で日付時分秒の文字列リストに変換します

> data$時刻 <- as.POSIXlt(paste(data$日付,data$時刻))

これで日付データもいらないのでデータフレームから削除します

> data <- data[, colnames(data) != "日付"]

colnamesが"日付"でない部分集合をdataに再代入しています

> data
                   時刻  始値  高値  安値  終値 出来高 売買代金
1   2016-06-18 03:00:00 15580 15580 15580 15580     82  1277560
...

こんな感じになります
ここまできたら終値をプロットできます

> plot(data$時刻, data$終値, type="l", xlab = "time", ylab = "closing price")

f:id:Kenjiqun:20160620173134p:plain

簡素なグラフが出来上がりました

最後にsummary関数を利用した簡単なデータ解析をしてみたいと思います

summary(データフレーム名[列番号の範囲])

としてやります
実際、日付や時刻の平均や最大値は無意味なので始値から売買代金までのデータをまとめてみます

> summary(data[2:7])
      始値            高値            安値            終値
 Min.   :15420   Min.   :15440   Min.   :15410   Min.   :15420
 1st Qu.:15510   1st Qu.:15520   1st Qu.:15500   1st Qu.:15510
 Median :15540   Median :15550   Median :15530   Median :15540
 Mean   :15554   Mean   :15564   Mean   :15544   Mean   :15554
 3rd Qu.:15590   3rd Qu.:15610   3rd Qu.:15580   3rd Qu.:15590
 Max.   :15720   Max.   :15740   Max.   :15700   Max.   :15730
 NA's   :3       NA's   :3       NA's   :3       NA's   :3
     出来高          売買代金
 Min.   :   0.0   Min.   :       0
 1st Qu.:  88.5   1st Qu.: 1372845
 Median : 166.0   Median : 2579150
 Mean   : 356.0   Mean   : 5549347
 3rd Qu.: 403.5   3rd Qu.: 6292590
 Max.   :5224.0   Max.   :81704200

Min.は最小値、1st Qu.は第1四分位数、Medianは中央値、Meanは平均値、3rd Qu.は第3四分位数、Max.は最大値となっています
また、NA'sはNAの出現回数を表しています
summary関数を用いることで簡単な統計量を求めることができました