Hatena::ブログ(Diary)

My Life as a Mock Quant このページをアンテナに追加 RSSフィード Twitter

2016-09-23

captionerパッケージで図・表に対する参照(レファレンス)を取得する

| 09:42 | captionerパッケージで図・表に対する参照(レファレンス)を取得するを含むブックマーク

http://datascienceplus.com/r-markdown-how-to-number-and-reference-tables/にあった話。

captionerパッケージを用いると、R Markdownの中での図表の参照を簡単に書くことができるというお話。

以下、簡単な例だが、こんな様に書けば図表に対する参照もできますよと。


---
title: ""
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Example of caption
Define caption.
```{r}
table_nums <- captioner::captioner(prefix = "Tab.")
tab.1_cap <- table_nums(name = "tab_1", caption = "German Bundesliga: Final Table 2015/16, Position 7-12")
```

Use it.

`r table_nums('tab_1')`
```{r fig.cap = tab.1_cap}
knitr::kable(mtcars[1:10, c(1,2,7,9)], align = c('c', 'l', 'c', 'c'), row.names = TRUE)
```

Define customized reference function.
```{r}
f.ref <- function(x) {
  stringr::str_extract(table_nums(x), "[^:]*")
}
```

Original version.
`r table_nums('tab_1')` says that ....

Customized one.
`r f.ref('tab_1')` says that ....

2016-08-28

リストに名前を付けるのは最後にした方がいいっぽい(へんなのR)

| 13:14 | リストに名前を付けるのは最後にした方がいいっぽい(へんなのR)を含むブックマーク

へんなのRというハッシュタグがある。

そこに投げつけるための話。

適当な関数の用意

3の時だけNULLであとはirisデータの頭を返す関数を用意する。

これはNULLを返す場合があるってのが本質で、あとは適当でよい。

#3の時だけNULLほかは適当
f <- function(i){
  if(i==3){
    NULL
  } else{
    head(iris,i)
  }
}

以下、この関数を適当にループさせ、その結果をリストに格納する形で使ってみる。

うまくいくケース

まずは何も考えずに使ってみたケース。

想像どおりの動き方&結果だ。

> x <- vector("list", 4)
> for(i in 1:4){
+   x[[i]] <- f(i)
+ }
> x
[[1]]
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa

[[2]]
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa

[[3]]
NULL

[[4]]
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa

うまくいかないケース

気を良くして、上のコードに一行、names(x)という処理を追加すると・・・名前がずれる。

> x <- vector("list", 4)
> names(x) <- paste0("name", 1:4)
> for(i in 1:4){
+   x[[i]] <- f(i)
+ }
> x
$name1
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa

$name2
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa

$name4
NULL

[[4]]
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa

原因と対策

にも書いてあるように、NULLを代入することは要素の削除と同義になっちゃうんでNGという話。

回避するには、この例だと最後に名前を割り当てるようにするか、あるいは汎用的な解決法としてはlapply使えと。

変なのー。

2016-08-27

Rcpp内で日付(Date)を取り扱う

| 18:24 | Rcpp内で日付(Date)を取り扱うを含むブックマーク

Rcpp内でDate型を取り扱いたいときにはDateVectorクラスを使用する。

ただ、

によると、Rcpp内で日付を扱うためのDateVectorクラスはそのうちdeprecatedになるらしいが、まぁ、まだあるので使う。

このベクトルの各要素はDateクラスのオブジェクトになっているので、そこは適当にこのドキュメントを見てこいつらを使ってやる。

例1

ためしに「引数として渡したDateのベクトルの日付(8月17日なら17)数分だけ日付を進めた日付ベクトル」をlistとして返すコードを書いてやると、

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List sample1(DateVector x) {
  List result(x.size());
  for(int i = 0; i < x.size(); ++i){
    DateVector date(x[i].getDay());
    for(int j = 0; j < x[i].getDay(); j++){
      date[j] = x[i] + j;
    }
    result[i] = date;
  }
  return result;
}

になる。

実際、動かしてみるとちゃんと動く。

> library("Rcpp")
> sourceCpp("hoge.cpp")
> sample1(x)
[[1]]
[1] "2016-08-28"

[[2]]
[1] "2016-08-29"

[[3]]
[1] "2016-08-30"
> sourceCpp("hoge.cpp")
> x <- Sys.Date() + 1:3
> x
[1] "2016-08-28" "2016-08-29" "2016-08-30"
> sample1(x)
[[1]]
 [1] "2016-08-28" "2016-08-29" "2016-08-30" "2016-08-31" "2016-09-01" "2016-09-02" "2016-09-03" "2016-09-04" "2016-09-05"
[10] "2016-09-06" "2016-09-07" "2016-09-08" "2016-09-09" "2016-09-10" "2016-09-11" "2016-09-12" "2016-09-13" "2016-09-14"
[19] "2016-09-15" "2016-09-16" "2016-09-17" "2016-09-18" "2016-09-19" "2016-09-20" "2016-09-21" "2016-09-22" "2016-09-23"
[28] "2016-09-24"

[[2]]
 [1] "2016-08-29" "2016-08-30" "2016-08-31" "2016-09-01" "2016-09-02" "2016-09-03" "2016-09-04" "2016-09-05" "2016-09-06"
[10] "2016-09-07" "2016-09-08" "2016-09-09" "2016-09-10" "2016-09-11" "2016-09-12" "2016-09-13" "2016-09-14" "2016-09-15"
[19] "2016-09-16" "2016-09-17" "2016-09-18" "2016-09-19" "2016-09-20" "2016-09-21" "2016-09-22" "2016-09-23" "2016-09-24"
[28] "2016-09-25" "2016-09-26"

[[3]]
 [1] "2016-08-30" "2016-08-31" "2016-09-01" "2016-09-02" "2016-09-03" "2016-09-04" "2016-09-05" "2016-09-06" "2016-09-07"
[10] "2016-09-08" "2016-09-09" "2016-09-10" "2016-09-11" "2016-09-12" "2016-09-13" "2016-09-14" "2016-09-15" "2016-09-16"
[19] "2016-09-17" "2016-09-18" "2016-09-19" "2016-09-20" "2016-09-21" "2016-09-22" "2016-09-23" "2016-09-24" "2016-09-25"
[28] "2016-09-26" "2016-09-27" "2016-09-28"

例2

次に、整数値としていったん受け取って、それをDateに直す方法。

これはRcppの中で単にRのas.Dateをimportして呼び出せるようにすればOK。

なんでこんな面倒なことをしないといけないのかというと、Dateクラスのオブジェクトには、数値”ベクトル”に対する加算演算子が定義されていないからだ。

#include <Rcpp.h>
using namespace Rcpp;
Function asDate("as.Date");
// [[Rcpp::export]]
List sample2(IntegerVector x, IntegerVector day) {
  List result(x.size());
  for(int i = 0; i < x.size(); ++i){
    result[i] = asDate(x[i] + seq_len(day[i]), "1970-01-01");
  }
  return result; 
}

これもちゃんと動く。

> sourceCpp("hoge.cpp")
> sample2(x, 1:3)
[[1]]
[1] "2016-08-29"

[[2]]
[1] "2016-08-30" "2016-08-31"

[[3]]
[1] "2016-08-31" "2016-09-01" "2016-09-02"

2016-08-17

三連ドット(..., dot-dot-dot, ellipsis)の取り扱い

| 16:22 | 三連ドット(..., dot-dot-dot, ellipsis)の取り扱いを含むブックマーク

これもr-wakalangに投げ込んで教えてもらった話なので、まずは簡単にまとめる。

基本的な使い方

適当な...を持つ関数を定義する。

f1 <- function(x, ...)
{
  dots <- list(...)
  print(dots)
}

これに対して以下の実行結果からわかるように、引数にマッチしなかったもの(ここではx以外)がリストのdots変数として関数内で使えるようになっているのがわかる。

要するに、多言語でいうところの可変長変数のようなもんだ。

list(...)という書き方の原点がどこにあるのかは不明だが、下記参考資料にあるR Language Definitionにも載ってるし、まぁこれはこういうもんかと思っておく。

> f1("a", 3)
[[1]]
[1] 3

> f1("b", a=3, b=7)
$a
[1] 3

$b
[1] 7

...のマッチに関する注意

普通のプログラミング言語だと可変長引数の後にはなんも書けない(引数を追加できない)ような気がするが、Rではできてしまう。

f2 <- function(x, ..., y)
{
  dots <- list(...)
  print(dots)
  print(paste0("y: ", y))
}

この場合、x, y以外の変数が全部...にマッチされている。

> f2("a", 3, z=7, y=3)
[[1]]
[1] 3

$z
[1] 7

[1] "y: 3"

R言語徹底解説にあったようなお話(data.frameに対するfilter操作)

R言語徹底解説のどこだかに載っていた気がする「data.frameに対するfilter操作」をdplyrぽく書く話も以下のように実現できる。

詳しくはoreorefilter関数のコメント参照。

oreorefilter <- function(df, ...)
{
  # 評価前に言語オブジェクトにしちゃう
  object <- as.list(substitute(list(...)))
  # どういう型として内部で保持されているかの確認
  print(str(object))
  # as.listしない場合の構造(つながった1リストになっちゃってる)
  print(str(substitute(list(...))))
  # 条件の順次評価(再帰で書いたほうがスマートか?)、一個目の要素はいらないのでループに入れない
  condition <- rep(TRUE, nrow(df))
  for(obj in object[-1]){
    # 条件を順次df環境で評価し、AND(&)条件として結ぶ
    condition <- condition & eval(obj, df)
  }
  df[condition, ]
}

実際に使ってみると、ちゃんとirisデータをFilteringできていることが分かる。

> oreorefilter(iris, Species=="setosa", Petal.Width==0.2, Sepal.Length==5.0)
List of 4
 $ : symbol list
 $ : language Species == "setosa"
 $ : language Petal.Width == 0.2
 $ : language Sepal.Length == 5
NULL
 language list(Species == "setosa", Petal.Width == 0.2, Sepal.Length == 5)
NULL
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
5             5         3.6          1.4         0.2  setosa
8             5         3.4          1.5         0.2  setosa
26            5         3.0          1.6         0.2  setosa
36            5         3.2          1.2         0.2  setosa
50            5         3.3          1.4         0.2  setosa

参考

2016-08-16

金融市場が混乱する間隔を日経平均ボラティリティー・インデックス(日経平均VI)から推定したい

| 18:48 | 金融市場が混乱する間隔を日経平均ボラティリティー・インデックス(日経平均VI)から推定したいを含むブックマーク

日経ボラティリティー・インデックス日経平均VI)祭りはまだ続く。

やりたいことは簡単で

だと考えて、"日経VIが跳ねるタイミング"の確率分布がどのような形になるのかを推定するというお話。

個人的には指数分布なんじゃね?と思っていたが、違った。

データの準備

ほぼ

でやった内容とおなじで、日経VIが跳ねたポイント(下図赤点)を手でやるのがめんどうなので、Rに同定させるためのコード。

↑のLINKとはrollapp()関数がちょっと違ってるので注意

library("Quandl")
# Sorce: https://www.quandl.com/data/NIKKEI/VLTL-Nikkei-Stock-Average-Volatility-Index
# Download Nikkei VI
vi <- Quandl("NIKKEI/VLTL")
#Decay detect
rollapp <- function(x, k, fun)
{
  w <- k/2
  sapply(seq_len(length(x)), function(i){
    index_s <- max(i-w, 1)
    index_e <- min(i+w, length(x))
    middle <- index_s + i-w
    fun(1+i-index_s, x[index_s:index_e])
  })
}
x <- vi$Close
max_point <- rollapp(x, 74, function(i, x){x[i]==max(x)})
max_point <- replace(max_point, is.na(max_point), FALSE)
plot(x, type="l")
points(seq_len(length(x))[max_point], x[max_point], col = "red", pch=16)

f:id:teramonagi:20160816185327p:image

このデータの赤点同士の時間間隔がどういう確率分布に従っているのかを評価したい。

推定(指数分布)

指数分布で確率密度関数を推定してみる

# 適当にdata.frame化
df <- dplyr::group_by(data.frame(index=cumsum(max_point)), index) %>% summarize(n=n())
lambda <- 1/mean(df$n)
# 下記でも同じ
# library("MASS")
# exponential <- coef(fitdistr(df$n, "exponential"))
hist(df$n, freq=FALSE)
curve(dexp(x, rate=lambda),add=TRUE,lty=2)

f:id:teramonagi:20160816185326p:image

う、うーん、これではない感。

推定(ワイブル分布)

によると、ワイブル分布は「時間に対する劣化現象や寿命を統計的に記述するためにも利用される」とあるので、こいつでやってみたい。

金融市場が破裂する=寿命が尽きると考えるわけだ。

library("MASS")
weibull <- coef(fitdistr(df$n, "weibull"))
hist(df$n, freq=FALSE)
curve(dweibull(x, weibull["shape"], weibull["scale"]), add=TRUE, lty=2)

f:id:teramonagi:20160816185325p:image

これだなー!

あとはこの裏側の寿命強度過程的なものがどうなっているのかを調べれば大体OKだな。

参考

日経平均ボラティリティー・インデックス(日経平均VI)の半減期が知りたい

| 15:26 | 日経平均ボラティリティー・インデックス(日経平均VI)の半減期が知りたいを含むブックマーク

引き続き日経VIネタ。

日経VIの挙動は、株価のような対数幾何分布に従うような確率過程(そうは見えないけどそうだと思っておく!)ではなく、平均回帰的であり、上昇するも下落するも急なものであるというのが(私的に)一般的に知られた事実である。

ここでは日経VIがどのくらいの速さで平均回帰するのか?、すなわち急上昇したのちにどれくらいの日数をかけて落ちるのかについて考えたい。

物理学では元の値の半分の大きさ(すなわち、金融でいうところの半値に相当)になるまでにかかる時間(日数)を半減期(とか、特徴的な長さスケール)だなんて呼ぶのでそれを計算してやる。

これがわかるとおおよそ日経VIがピークになってからその半分の値になるまでの時間(日数)がわかるのである。

データの取得

データは従前の記事同様、Quandlパッケージを活用して取得する。

library("Quandl")
# Sorce: https://www.quandl.com/data/NIKKEI/VLTL-Nikkei-Stock-Average-Volatility-Index
# Download Nikkei VI
vi <- Quandl("NIKKEI/VLTL")

ピーク位置の同定

これが今回の一番悩んだところ。

ここでは以下のアイディアに基づいて抽出しているが、これが一般的なのかどうかはさっぱりわからない*1

後から調べた感じ、

この辺を参考にパッケージ使えばよかった…。


ここではまず、時系列データの各点を中心として半径(単系列での半径なので、ただの一次元円の半径、すなわちただの距離)k/2以内にあるデータに対してfunという処理をかませるrollapp()関数を定義する。

これは、時系列の一番後ろ(時間的に昔)からデータを舐めて処理するzooパッケージにあるrollapplyとは違い、各時点の未来・過去k/2単位(今回の場合日次データなので、単位は日)にまたがったデータに対して処理を適用する関数

#Decay detect
rollapp <- function(x, k, fun)
{
  w <- k/2
  sapply(seq_len(length(x)), function(i){
    if((i-w) <= 0 || (i+w) > length(x)){
      NA
    } else{
      fun(x[(i-w):(i+w)])
    }
  })
}

こいつを使って

  • 上値ピークポイント:直近150日(未来・過去方向それぞれ75日)分のデータの中で最も大きく、かつ、日経VIの値自体は30以上
  • 下値ピークポイント:直近125日(未来・過去方向それぞれ62日)分のデータの中で最も小さい

点を同定し、それぞれmax_point, min_pointとする。

NAはよしなに削除で。

x <- vi$Close
max_point <- rollapp(x, 150, function(x){x_middle <- x[1+length(x)/2]; (x_middle==max(x)) && (x_middle > 30)})
max_point <- replace(max_point, is.na(max_point), FALSE)
min_point <- rollapp(x, 125, function(x){x_middle <- x[1+length(x)/2]; (x_middle==min(x))})
min_point <- replace(min_point, is.na(min_point), FALSE)

ここで計算したピークポイントがどの辺にあるのかというと、下図のとおり、まぁ見た感じよさそうだぞと。

(上(下)値ピークが赤(青)点、日経VIの値は黒線で表示)

(逆に、こういう見た目になるように↑の日数とか日経VIの水準を足きりしている点に注意)

plot(x, type="l")
points(seq_len(length(x))[max_point], x[max_point], col = "red", pch=16)
points(seq_len(length(x))[min_point], x[min_point], col = "blue", pch=16)

f:id:teramonagi:20160816153707p:image

上昇→下落フェーズのデータの抽出

ここで、「定義した上値・下値の間にあるデータ」のみを抽出する。

つまり上の図でいうところの、各赤点・青点の間にある日経VIのデータを抽出するということだ。

そのためには汚いが、データを全舐めして、下記のようなコードでデータを抽出する。

index変数が0以外の値で、かつ、共通の値を持つ箇所がおなじ上昇→下落局面に対応している。

index <- numeric(length(x))
counter <- 0
mode <- "max"
for(i in seq_len(length(x))){
  if(mode == "max" && max_point[i]){
    mode <- "min"
    counter <- counter + 1
    index[i] <- counter
  } else if(mode == "min" && min_point[i]){
    mode <- "max"
  } else if(mode == "min"){
    index[i] <- counter
  }
}
xs <- lapply(seq_len(max(index)), function(i){x[which(index == i)]})

結果は、各赤点・青点の間にあるデータごとのリストになっている。

半減期の同定

上で作ったリストの各要素が上昇→下落局面に対応する

 NIKKEI-VI = a*¥exp¥{-t/¥tau¥}

という回帰式に対してτを推定することがメインであり、このτが半減期、すなわち日経VIがどの程度の速度で下落するのかの指標になっている*2

推定する部分を関数にし、元データ・τ・推定結果(データ自体をlog化し、そのデータの誤差が正規分布になってほしいなと仮定して、線形回帰に直して計算させる)を返却させるようにしよう。

# x = a*exp(-t/tau)
# log(x) = log(a) - 1/tau * t
decay <- function(x) {
  df <- data.frame(t=seq_len(length(x)), logx=log(x))
  regression <- lm(logx~t, data=df)
  tau <- -1/coef(regression)["t"]
  list(x=x, tau=tau, regression=regression)
}

これを、上で作った全データに適用してやると

result <- lapply(xs, decay)

となる。

得られたτは

> tau <- sapply(result, function(x){x$tau})
> tau
        t         t         t         t         t         t         t         t         t 
 56.60291  58.95572 165.28740 395.43998 201.99652 109.37092  40.24398  33.13756  27.06593 

となり、30くらいから400くらいまでとかなり幅広だ。

各々の結果を出したデータと回帰の当てはまり具合を下記のように目視で確認すると、長い半減期の系はだいたいフィッティングがおかしいので、無視してよさそうだ。

これは、どちらかというと、データの中に小さなピークが立っていることからもわかるように、ピークの同定がうまくないからなんだろうなと。

従って、短い方を意識して大体、半減期は10-30日程度と意識して日経VI先物投資すればよいだろう(と期待したい)。

τ=56.60291

f:id:teramonagi:20160816160309p:image

τ=58.95572

f:id:teramonagi:20160816160308p:image

τ=165.28740

f:id:teramonagi:20160816160307p:image

τ=395.43998

f:id:teramonagi:20160816160306p:image

τ=201.99652

f:id:teramonagi:20160816160304p:image

τ=109.37092

f:id:teramonagi:20160816160303p:image

τ=40.24398

f:id:teramonagi:20160816160302p:image

τ=33.13756

f:id:teramonagi:20160816160301p:image

τ=27.06593

f:id:teramonagi:20160816160300p:image

*1:まじめに調べてもいない

*2:より正確にはこの値の0.7倍程度か https://ja.wikipedia.org/wiki/%E5%8D%8A%E6%B8%9B%E6%9C%9F