今週も特にありません

進捗どうですか?

tidyquantを使った分析メモ(2)

tidyquantを使った分析メモ

今回は対数リターンの計算とその可視化までをやってみます。

github.com

masaqol.hatenablog.com

対数リターンの計算

TOPIX100構成銘柄の株価データを取得したい場合、銘柄コードを用意して、tidyquanttq_getを使うことで簡単に取得することができるところまで前回確認しました。

> TOPIX100
# A tibble: 6,100 × 8
     symbol       date   open   high    low  close  volume
      <chr>     <date>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl>
1  YJ1605.T 2017-01-04 1182.5 1201.0 1180.0 1190.0 5927100
2  YJ1605.T 2017-01-05 1186.0 1187.0 1148.5 1168.0 8302500
3  YJ1605.T 2017-01-06 1139.0 1158.5 1131.0 1153.5 7737100
4  YJ1605.T 2017-01-10 1141.5 1153.0 1136.0 1136.0 4807700
5  YJ1605.T 2017-01-11 1140.5 1156.5 1133.5 1151.5 5587400
6  YJ1605.T 2017-01-12 1143.0 1156.5 1134.5 1145.5 4831300
7  YJ1605.T 2017-01-13 1130.5 1152.5 1130.5 1150.5 5257000
8  YJ1605.T 2017-01-16 1144.5 1150.0 1127.0 1131.5 3200200
9  YJ1605.T 2017-01-17 1136.0 1138.0 1117.0 1122.0 3625500
10 YJ1605.T 2017-01-18 1116.0 1146.5 1112.0 1143.0 5458900
# ... with 6,090 more rows, and 1 more variables: adjusted <dbl>

株価そのままのデータでは、銘柄ごとの比較などが行いにくいために、対数リターンを計算します。

tidyquantは複数のファイナンス関係のパッケージの便利なラッパーとして機能します。 対数リターンなどの基本的な計算は、元のパッケージの便利な関数を使うことでtidyなデータに対して処理を行えます。

tq_mutateによって、新しい変数として対数リターンを追加します。 適用できる関数の一覧はtq_mutate_fun_optionsによって確認できます。

> tq_mutate_fun_options()
$zoo
 [1] "rollapply"          "rollapplyr"        
 [3] "rollmax"            "rollmax.default"   
 [5] "rollmaxr"           "rollmean"          
 [7] "rollmean.default"   "rollmeanr"         
 [9] "rollmedian"         "rollmedian.default"
[11] "rollmedianr"        "rollsum"           
[13] "rollsum.default"    "rollsumr"          

$xts
 [1] "apply.daily"     "apply.monthly"   "apply.quarterly"
 [4] "apply.weekly"    "apply.yearly"    "diff.xts"       
 [7] "lag.xts"         "period.apply"    "period.max"     
[10] "period.min"      "period.prod"     "period.sum"     
[13] "periodicity"     "to_period"       "to.daily"       
[16] "to.hourly"       "to.minutes"      "to.minutes10" 
...

対数リターンの計算は、quantmodperiodReturnを利用します。

log_returns <- TOPIX100 %>%
  group_by(symbol) %>%
  tq_mutate(
    select = close, 
    mutate_fun = periodReturn,
    period = "daily",
    type = "log",
    col_rename = "log_return"
  ) %>%
  select(symbol, date, close, log_return)

計算したものを確認します。

> log_returns
Source: local data frame [6,100 x 4]
Groups: symbol [100]

     symbol       date  close   log_return
      <chr>     <date>  <dbl>        <dbl>
1  YJ1605.T 2017-01-04 1190.0  0.000000000
2  YJ1605.T 2017-01-05 1168.0 -0.018660423
3  YJ1605.T 2017-01-06 1153.5 -0.012492086
4  YJ1605.T 2017-01-10 1136.0 -0.015287478
5  YJ1605.T 2017-01-11 1151.5  0.013552120
6  YJ1605.T 2017-01-12 1145.5 -0.005224217
7  YJ1605.T 2017-01-13 1150.5  0.004355408
8  YJ1605.T 2017-01-16 1131.5 -0.016652444
9  YJ1605.T 2017-01-17 1122.0 -0.008431379
10 YJ1605.T 2017-01-18 1143.0  0.018543578
# ... with 6,090 more rows

データの可視化

対数リターンを時系列で可視化します。(銘柄コードが7000番台の銘柄のみ)

log_returns %>% 
  filter(grepl("^YJ7", symbol)) %>%
  ggplot(aes(x = date, y = log_return, colour = symbol)) +
    geom_line() +
    theme_tq() +
    scale_color_tq() +
    facet_wrap(~ symbol, ncol = 3, scales = "free_y")

f:id:masaqol:20170722111530p:plain

対数リターンの分布を確認します。(銘柄コードが7000番台の銘柄のみ)

log_returns %>% 
  filter(grepl("^YJ7", symbol)) %>%
  ggplot(aes(x = log_return, fill = symbol)) +
    geom_histogram(binwidth = 0.005) +
    theme_tq() +
    scale_color_tq() +
    facet_wrap(~ symbol, ncol = 3, scales = "free_y")

f:id:masaqol:20170722111649p:plain

期待リターンとボラティリティ

最後に各銘柄の期待リターンとボラティリティを計算し、高い銘柄を確認します。

> return_vol <- log_returns %>%
+     group_by(symbol) %>%
+     summarise(return = mean(log_return), vol = sd(log_return))
>
> return_vol %>%
+     arrange(desc(return)) %>% 
+     head(10)
# A tibble: 10 × 3
     symbol      return         vol
      <chr>       <dbl>       <dbl>
1  YJ6273.T 0.002274630 0.016555097
2  YJ2502.T 0.002011907 0.012576164
3  YJ6758.T 0.002002305 0.013206705
4  YJ6954.T 0.001991144 0.013172032
5  YJ4452.T 0.001708134 0.010844302
6  YJ6861.T 0.001573498 0.009784067
7  YJ4188.T 0.001551808 0.016332976
8  YJ2503.T 0.001353853 0.011863891
9  YJ7741.T 0.001352014 0.011299796
10 YJ8035.T 0.001253836 0.014759141
> 
> return_vol %>%
+   arrange(desc(vol)) %>% 
+   head(10)
# A tibble: 10 × 3
     symbol        return        vol
      <chr>         <dbl>      <dbl>
1  YJ6502.T -2.278774e-03 0.05451714
2  YJ8795.T -8.095544e-05 0.02178426
3  YJ7186.T -2.028825e-03 0.02028579
4  YJ8750.T -1.430680e-04 0.01997049
5  YJ4755.T -7.893335e-04 0.01927645
6  YJ9064.T -7.722044e-04 0.01899148
7  YJ7974.T  8.831755e-04 0.01850704
8  YJ4578.T -1.365315e-04 0.01812934
9  YJ8630.T  8.039943e-06 0.01787958
10 YJ9983.T -3.522327e-03 0.01776665

期待リターントップは、6273のSMC。ボラティリティトップは6502の東芝ということがわかりました。

tidyquantを使った分析メモ

R-bloggersでも度々登場して気になっていたtidyquantを使った分析メモ

github.com

基本的な使い方は、Vignettesが充実しています。 また、日本語での資料としてはtidyquantの使い方があります。 また、Exploratoryと組み合わせて分析した記事もあります。 qiita.com

ここでは、日本の株式銘柄をそれなりに大量取得して、データを可視化するところまでやってみます。

銘柄コード一覧を作成

まず、株価データを取得する際に必要な銘柄コード一覧を作成します。 ここでは、TOPIX100構成銘柄をまとめて取得することを考えます。

銘柄の入れ替え等があった場合に適宜更新を行ってくれると思われる 松井証券さんのページを信用できる情報として、 rvestでスクレイプして作成します。

library(tidyverse)
library(tidyquant)
library(rvest)

url <- "http://www.matsui.co.jp/service/stock/trade/rule_topix100/"
html <- read_html(url, encoding = "UTF-8")

stock_name_df <- html %>% 
  html_node("table.m-tbl") %>%
  html_table

htmlの中からm-tblクラスのtable部分をdata.frameとして取得しています。

> stock_name_df
    銘柄コード                                        銘柄名
1         1605                              国際石油開発帝石
2         1878                                      大東建託
3         1925                                大和ハウス工業
4         1928                                    積水ハウス
5         2502                アサヒグループホールディングス
6         2503                        キリンホールディングス
7         2802                                        味の素
8         2914                                日本たばこ産業
9         3382                 セブン&アイ・ホールディングス
10        3402                                          東レ
11        3407                                        旭化成
...

データ取得で必要なのは銘柄コードのみのため、これだけを切り出せば完了です。

stock_code <- stock_name_df %>%
  select(`銘柄コード`)

データを確認します。

> stock_code
    銘柄コード
1         1605
2         1878
3         1925
4         1928
5         2502
6         2503
7         2802
8         2914
9         3382
10        3402
11        3407
...

株価データの取得

銘柄コード一覧が作成できたら、それをもとに株価データを取得します。 ここでは、tq_get()get = "stock.prices.japan"を指定して、ヤフーファイナンスからデータを取得することにします。

この際に銘柄コードは、"YJ9984.T"のようにして問い合わせれば、データを取得できます。 (取得に時間がかかるために、3ヶ月分としています。)

TOPIX100 <- stock_code %>%
  unlist %>%
  as.character %>%
  paste0("YJ", ., ".T") %>%
  tq_get(get = "stock.prices.japan",
         from = "2017-01-01",
         to = "2017-03-31")

データを確認します。

> TOPIX100
# A tibble: 6,100 × 8
     symbol       date   open   high    low  close  volume
      <chr>     <date>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl>
1  YJ1605.T 2017-01-04 1182.5 1201.0 1180.0 1190.0 5927100
2  YJ1605.T 2017-01-05 1186.0 1187.0 1148.5 1168.0 8302500
3  YJ1605.T 2017-01-06 1139.0 1158.5 1131.0 1153.5 7737100
4  YJ1605.T 2017-01-10 1141.5 1153.0 1136.0 1136.0 4807700
5  YJ1605.T 2017-01-11 1140.5 1156.5 1133.5 1151.5 5587400
6  YJ1605.T 2017-01-12 1143.0 1156.5 1134.5 1145.5 4831300
7  YJ1605.T 2017-01-13 1130.5 1152.5 1130.5 1150.5 5257000
8  YJ1605.T 2017-01-16 1144.5 1150.0 1127.0 1131.5 3200200
9  YJ1605.T 2017-01-17 1136.0 1138.0 1117.0 1122.0 3625500
10 YJ1605.T 2017-01-18 1116.0 1146.5 1112.0 1143.0 5458900
# ... with 6,090 more rows, and 1 more variables: adjusted <dbl>

データの可視化

取得したデータを可視化します。全ての銘柄をプロットすると見にくいため、 銘柄コードが7000番台の銘柄のみの株価推移を可視化します。

TOPIX100 %>% 
  filter(grepl("^YJ7", symbol)) %>%
  ggplot(aes(x = date, y = close, colour = symbol)) +
    geom_line() + 
    theme_tq() +
    scale_color_tq() +
    facet_wrap(~ symbol, ncol = 3, scales = "free_y")

f:id:masaqol:20170521123132p:plain

日本国債金利の主成分分析で2016年振り返り

マイナス金利や長短金利操作付き量的・質的金融緩和導入など、今年も何かと金融政策の変更とそれに伴う金利の変動がニュースを賑わせた一年でした。 久々に何か分析ネタを…ということで、日本国債金利を主成分分析することで、金利の構造の主なファクターが2016年はどのように推移したのかを調べたいと思います。

データはQuandlから持ってきます。 Quandlパッケージを使えば、data.framextsなどの形でデータを使えるようになり、手軽に分析を始めることができます。

まずは、必要なライブラリを読み込みます。

library(Quandl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)

Quandlからデータを取ってくるには、データを表すコードと取得開始日を最低限入れます。 ここでは、日本国債金利データをそのままを使うことにします。 金利の変動について分析を行いたい場合は、Quandl関数のtransformdiffを指定することで取ってくることができます。

irate <- Quandl("MOFJ/INTEREST_RATE_JAPAN", order = "asc", start_date = "2016-01-01")

取ってきたデータを確認します。

> head(irate)
        Date     1Y     2Y     3Y     4Y    5Y    6Y
1 2016-01-04 -0.042 -0.028 -0.014  0.003 0.026 0.036
2 2016-01-05 -0.043 -0.030 -0.014  0.003 0.026 0.036
3 2016-01-06 -0.044 -0.029 -0.014  0.002 0.022 0.032
4 2016-01-07 -0.045 -0.028 -0.014 -0.003 0.016 0.027
5 2016-01-08 -0.045 -0.031 -0.018 -0.003 0.016 0.022
6 2016-01-12 -0.048 -0.028 -0.018 -0.003 0.017 0.016
     7Y    8Y    9Y   10Y   15Y   20Y   25Y   30Y   40Y
1 0.059 0.120 0.186 0.264 0.599 0.986 1.150 1.282 1.402
2 0.059 0.120 0.181 0.259 0.588 0.971 1.135 1.264 1.386
3 0.055 0.115 0.176 0.253 0.583 0.967 1.131 1.260 1.382
4 0.049 0.110 0.171 0.245 0.573 0.958 1.116 1.238 1.363
5 0.039 0.101 0.161 0.230 0.557 0.939 1.096 1.220 1.346
6 0.034 0.097 0.157 0.221 0.548 0.930 1.086 1.207 1.338

1年から40年満期までの金利データが取れました。 早速、データを可視化します。

irate %>% 
  gather(Maturity, Rate, -Date) %>%
  mutate(Maturity = factor(.$Maturity, levels = paste0(c(1:10, 15, 20, 25, 30, 40), "Y"))) %>%
    ggplot(aes(x = Date, y = Rate)) +
    geom_line(aes(colour = Maturity), size = 1) + 
    scale_x_date(labels = date_format("%y/%m"), date_breaks = "3 month") +
    guides(colour = guide_legend(reverse = TRUE))

f:id:masaqol:20161231123705p:plain

1月のマイナス金利導入以降、イールドカーブ全体的に低下傾向となり、2月終盤からは10年満期以下の国債金利が全てマイナス圏で推移するようになります。 7月には15年満期の国債金利もマイナスという状態も見られるようになり、年間を通して一番金利が低下した時期であったことがわかります。 それ以降は全体的に上昇し、アメリカの大統領選挙後には、10年満期の国債金利もマイナス圏から脱し、金利の上昇が加速しているように見えます。

それでは、このデータに対して主成分分析を行います。 これらの分析については、様々なところで行われているため詳細な説明は省きますが、金利は水準(Level)、傾き(Slope)、曲率(Curveture)の3つのファクターでほぼその構造が説明できると言われています。主成分分析によって、これらのファクターが導出されるのかを確認し、その時系列推移を可視化してみます。

特にひねりを加えるわけでもなく、Rのprcomp金利データに適用します。

irate_pc <- irate %>% 
  select(-Date) %>%
  prcomp(scale. = TRUE)

主成分分析を適用した結果を見てみます。

> summary(irate_pc)
Importance of components:
                         PC1     PC2     PC3     PC4     PC5
Standard deviation     3.654 0.99449 0.75082 0.24082 0.13649
Proportion of Variance 0.890 0.06593 0.03758 0.00387 0.00124
Cumulative Proportion  0.890 0.95589 0.99347 0.99734 0.99858
                           PC6     PC7     PC8     PC9    PC10
Standard deviation     0.09270 0.07582 0.05226 0.03636 0.03167
Proportion of Variance 0.00057 0.00038 0.00018 0.00009 0.00007
Cumulative Proportion  0.99915 0.99953 0.99972 0.99980 0.99987
                          PC11    PC12    PC13    PC14    PC15
Standard deviation     0.02533 0.02314 0.01826 0.01654 0.01175
Proportion of Variance 0.00004 0.00004 0.00002 0.00002 0.00001
Cumulative Proportion  0.99991 0.99995 0.99997 0.99999 1.00000

第3主成分までで累積寄与率が99%となり、ほぼこの3つでイールドカーブの構造を説明できていることがわかります。 主成分得点は以下のようになりました。

> head(irate_pc$x)
           PC1        PC2       PC3         PC4        PC5
[1,] -9.534255 -0.3139355 0.4521807 -0.16670843 -0.2093294
[2,] -9.425388 -0.2323307 0.4851745 -0.16565373 -0.1661761
[3,] -9.329534 -0.2554974 0.5167301 -0.11327598 -0.1747176
[4,] -9.152626 -0.2405698 0.5779134 -0.08359126 -0.2015301
[5,] -8.915227 -0.2231325 0.6567429 -0.03116720 -0.1449082
[6,] -8.801091 -0.1916643 0.7010752  0.05046791 -0.1497915
               PC6          PC7          PC8         PC9         PC10
[1,] -0.0451359411  0.057984040 -0.006210947  0.03531776  0.006980638
[2,] -0.0351578833  0.054768411  0.009829469  0.03076594  0.008573947
[3,] -0.0347377299  0.033557155  0.012547679  0.03020654  0.010608943
[4,] -0.0531471021 -0.004868424  0.018739106  0.05401370 -0.006719837
[5,] -0.0005787395  0.023142080  0.005714631  0.02306513 -0.002223937
[6,] -0.0028323655  0.009862543 -0.014780631 -0.01505679 -0.007297795
            PC11        PC12          PC13         PC14       PC15
[1,] -0.02586290 -0.03433084 -0.0009471546 -0.016075984 0.02848080
[2,] -0.03508559 -0.03784112  0.0147360705 -0.012732611 0.02126039
[3,] -0.02178643 -0.04119171  0.0143781817 -0.010653086 0.02193993
[4,] -0.02087492 -0.02625881  0.0161335068 -0.015374304 0.01088626
[5,] -0.01278054  0.01291703  0.0168199370 -0.008567272 0.02197490
[6,] -0.02024915  0.05001019  0.0143638963 -0.011216030 0.01494438

そして、各満期に対するファクターローディングスを見てみます。

data.frame(irate_pc$rotation[, 1:3]) %>%
  mutate(Maturity = factor(rownames(.), levels = rownames(.))) %>%
    gather(PC, FactorLodings, -Maturity) %>%
    ggplot(aes(x = Maturity, y = FactorLodings)) + 
    geom_line(aes(colour = PC, group = PC), size = 1) 

f:id:masaqol:20161231124810p:plain

第1から第3主成分までが金利の水準、曲率、長短期の傾きをそれぞれ表していることがわかります。 (横軸の間隔が満期に対応して一定ではないので、形がちょっと変な感じですが…)

さらに、これらのファクターの時系列推移も可視化して、一年を通してどのように変動したのかを見てみます。

data.frame(Date = irate$Date, irate_pc$x[, 1:3]) %>% 
  gather(PC, Factor, -Date) %>%
    ggplot(aes(x = Date, y = Factor)) + 
    geom_line(aes(colour = PC), size = 1) +
    scale_x_date(labels = date_format("%y/%m"), date_breaks = "3 month")

f:id:masaqol:20161231125156p:plain

1月のマイナス金利導入以降、金利の水準が大きく下落し、傾きや曲率も変動したことがわかります。 そして、年後半には再び水準が上昇したことがわかります。 年間を通して、金利が下落した後に上昇して戻ってきたと言っても、 水準と傾き、曲率の3つのファクターを見ると、年初と年末ではイールドカーブの構造が異なっていることがわかります。

ただ単純に主成分分析しただけのものなので、ガチなフラットナー/スティープナーといった金利関連の取引に役立つわけではないと思いますが、データの取得から可視化、分析手法適用までRを使って簡単に始められますというものでした。

来年も世界各国の政治イベントで大きな動きがありそうということで、勝手にデータ分析の材料を提供してくれると思うと楽しめそうです。

Rにおける最適化アルゴリズムをこれ一つで?

ほぼ使われているところを見ないoptimxパッケージを使ってみたメモ

データ分析を行っている方は、日々パラメータの推定などで最適化に励んでいると思います

Rで目的関数を自分で定義して最適化しようとする場合、optimやnlminb(一変数の場合はoptimize)等を用いて最適化を行います。 optimで最適化する場合、とりあえずデフォルトで実行されるNelder-Mead法を試して、この最適化が上手くいかなかった場合に、BFGS法などのその他の方法を試すという状況がよくあると思います

optimxパッケージは、Rで実行可能な最適化手法をまとめて一通り実行して、その結果を比較できるものです

尤度関数を定義して、そのパラメータを推定する例をやってみます

library(optimx)

loglik <- function(x, param) {
  -sum(dnbinom(x, size = param[1], mu = param[2], log = TRUE))
}

set.seed(71)
data <- rnbinom(500, size = 10, mu = 5)
initparam <- c(10, 5)
optimx(par = initparam, fn = loglik, x = data, control = list(all.methods = TRUE))

いくつかの警告文が出てきますが、最適化結果は以下のようになります

                   p1       p2         value fevals gevals niter convcode kkt1 kkt2 xtimes
BFGS         9.889487 5.038007  1.190725e+03     14      5    NA        0   NA   NA   0.01
CG           9.892441 5.037995  1.190725e+03    305    101    NA        1   NA   NA   0.27
Nelder-Mead  9.892127 5.038049  1.190725e+03     41     NA    NA        0   NA   NA   0.02
L-BFGS-B     9.889511 5.038000  1.190725e+03     13     13    NA        0   NA   NA   0.03
nlm          9.889480 5.037999  1.190725e+03     NA     NA    11        0   NA   NA   0.01
nlminb       9.889512 5.038000  1.190725e+03      8     18     5        0   NA   NA   0.00
spg          9.889768 5.038000  1.190725e+03     24     NA    18        0   NA   NA   0.13
ucminf       9.889507 5.037997  1.190725e+03      9      9    NA        0   NA   NA   0.02
Rcgmin             NA       NA 8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
Rvmmin             NA       NA 8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
newuoa       9.889511 5.038000  1.190725e+03     34     NA    NA        0   NA   NA   0.02
bobyqa       9.889509 5.038000  1.190725e+03     27     NA    NA        0   NA   NA   0.01
nmkb         9.890341 5.037991  1.190725e+03     69     NA    NA        0   NA   NA   0.03
hjkb        10.000000 5.000000  1.190775e+03      1     NA     0     9999   NA   NA   0.02

関数自体の使い方は普通のoptimとほとんど変わりません。 all.methods=TRUEを指定することですべての手法を実行します。 p1、p2が最適化したパラメータの値、valueがその時の目的関数の値で、 convcodeが0以外のものは最適化が正常に終了していないことを表しています

今回の場合では、結果が大きく変わることはないですが、 より複雑な場合だとそれぞれのアルゴリズムを一覧で比較できるところは使いどころかもしれません

より詳しい説明はJournal of statistical softwareにもあります。 CRAN Task View: Optimization and Mathematical Programmingを見ると、知らないパッケージがたくさんあるなー、まだまだRはいろいろとできるのだなーと感じます

一般化線形モデル周りのRのパッケージと関数まとめ

(GLM、GLMM周りを中心に)統計モデリングのためのパッケージと関数のメモ 使い方の詳しい説明、利用事例等を見つけたら、順次追加予定

lm {stats} Fitting Linear Models

線形モデル

nls {stats} Nonlinear Least Squares

非線形モデル

glm {stats} Fitting Generalized Linear Models

一般化線形モデル。線形回帰、ロジスティクス回帰、ポアソン回帰等の推定が可能。 cv.glm {boot}でクロスバリデーションを実行可能

glm.nb {MASS} Fit a Negative Binomial Generalized Linear Model

負の二項分布を仮定する場合の一般化線形モデル

gam {mgcv} Generalized additive models with integrated smoothness estimation

一般化加法モデル

gamm {mgcv} Generalized Additive Mixed Models

一般化加法混合モデル

plm {plm} Panel Data Estimators

パネルデータに対する線形モデル

polr {MASS} Ordered Logistic or Probit Regression

順序ロジット。プロビットモデルも可。polrは比例オッズモデルから来ている

mlogit {mlogit} Multinomial logit model

多項ロジットモデル。プロビットモデルも利用可能。(Kenneth Train’s exercises using the mlogit package for R

mnlogit {mnlogit} Fast estimation of multinomial logit models

多項ロジットモデル。パッケージ作者による詳しい説明あり(Fast Estimation of Multinomial Logit Models: R Package mnlogit

npmlt {mixcat} Mixed effects cumulative link and logistic regression models

多項ロジットモデル

multinom {nnet} Fit Multinomial Log-linear Models

多項ロジットモデル

gmnl {gmnl} Estimate Multinomial Logit Models with Observed and Unobserved Individual Heterogeneity.

多項ロジットモデル。詳しい説明はこちら(gmnl package in R - Mauricio Sarrias

lmer {lme4} Fit Linear Mixed-Effects Models

線形混合モデル。ランダム効果を複数指定可能。Journal of Statistical Softwareにパッケージの詳しい説明あり(Fitting Linear Mixed-Effects Models Using lme4 | Bates | Journal of Statistical Software

glmer {lme4} Fitting Generalized Linear Mixed-Effects Models

一般化線形混合モデル。ランダム効果を複数指定可

nlmer {lme4} Fitting Nonlinear Mixed-Effects Models

非線形混合モデル

glmer.nb {lme4} Fitting GLMM’s for Negative Binomial

負の二項分布を仮定する場合の一般化線形混合モデル

glmmML {glmmML} Generalized Linear Models with random intercept

一般化線形混合モデル。ランダム効果は一つのみ指定可

glmnet {glmnet} fit a GLM with lasso or elasticnet regularization

lasso、ridge、elastic netによる正則化付き一般化線形モデル。線形回帰、ポアソン回帰、(多項)ロジスティクス回帰等の推定が可能。group lassoなども追加されてきている。cv.glmnet {glmnet}でクロスバリデーションを実行可能

glmreg {mpath} fit a GLM with lasso (or elastic net), snet or mnet regularization

正則化付き一般化線形モデル。lasso、scad、mcpなどを用いることができる
cv.glmreg {mpath}でクロスバリデーションを実行可能

zeroinfl {pscl} Zero-inflated Count Data Regression

Zero-inflatedモデル

vglm {VGAM} Fitting Vector Generalized Linear Models

ベクトル一般化線形モデル。Version 1.0-1の時点で「Currently only fixed-effects models are implemented, i.e., no random-effects models.」ということで、ランダム効果は考慮できない。Springerの本をもとにしている(Vector Generalized Linear and Additive Models - With an | Thomas W. Yee | Springer

vgam {VGAM} Fitting Vector Generalized Additive Models

ベクトル一般化加法モデル

xtsでmatplotのようなプロット

最近はdplyrやtidyrとの相性で、ほぼggplot2を使うようになり、めっきりmatplotを使う機会が無くなっている

しかし、時系列データを扱う際にdata.frameよりもxtsを使いたい場合もあり、その時にxtsのままでmatplotのような感じのプロットするにはどうするのがいいのか?と思っていたので調べたメモ。結局、xtsとxtsExtraをセットで使うのが良さそう

まず、データとして、複数時系列データをQuandlから取ってくる。この時に、Quandl関数でtypeにxtsを指定すれば、xtsクラスでデータを取って来れるが、一旦、ggplot2でプロットするためにdata.frameで取って来る

library(Quandl)

rate_df <- Quandl("MOFJ/INTEREST_RATE_JAPAN", trim_start = "2015-01-01")

このデータをggplot2でプロットする

library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)

rate_df %>% gather(maturity, rate, -Date) %>% 
  ggplot(aes(x = Date, y = rate, colour = maturity)) + geom_line(size = 1) +
  scale_x_date(labels = date_format("%m/%d"))

f:id:masaqol:20151206233544p:plain

その他、色々と付け足していけば、思いどおりのプロットができるだろう。しかし、本題はデータをdata.frameではなくxtsクラスで扱いたい場合。まわりくどいがQuandlの中でやっている処理と同じようにdata.frameからxtsに変換する

rate_xts <- xts(rate_df[, -1], order.by = rate_df[, 1])

ここで、そのままplotすれば複数系列を一辺にプロットしてくれると思いきや、警告メッセージとともに残念なプロットがされる

> plot(rate_xts)
Warning message:
In plot.xts(rate_xts) : only the univariate series will be plotted

f:id:masaqol:20151206234012p:plain

matplotならばどうかというと、xtsプロットの良い部分は完全に失われる

> matplot(rate_xts, type = "l")

f:id:masaqol:20151206234802p:plain

xtsはzooクラスとしても扱えるのだから、plot.zooで良い感じにプロットされるのではないかと思い、やってみるとmatplotとほぼ変わらず…

> plot.zoo(rate_xts, plot.type = "single")

f:id:masaqol:20151206235744p:plain

というところまで試してみて、CRANには上がっていないxtsExtraを使うのが一番良さそうということが分かり、使ってみるとxtsプロットの良い部分もあり、matplotのような複数系列のプロットが簡単にできた

library(devtools)
devtools::install_github("joshuaulrich/xtsExtra")
library(xtsExtra)

plot(rate_xts, col = 1:15)

f:id:masaqol:20151207004738p:plain

ただ、デフォルトでplotするだけでは12系列までしかプロットされないようなので、色の指定を系列数分に増やす、または、addSeriesを使って追加する必要があるようです。matplotしたものの軸ラベル等をxtsプロット風に頑張って変更するよりも、断然xtsExtraを使うのが良さそうです

pforeachで最尤推定した漸近分布を確認する

最尤推定で推定したパラメータの漸近分布をシミュレーションによって確認する
繰り返し推定を行う際にpforeach(hoxo-m/pforeach · GitHub)を使う

CIRモデルのパラメータ推定と関数もパラメータの設定も同じとする

r0 <- 0.01
n <- 1000
delta <- 1/52
kappa <- 0.2
mu <- 0.03
sigma <- 0.08
initparam <- c(kappa, mu, sigma)

パスを発生させ、パラメータ推定を繰り返す。収束しなかったものは取り除く

ggplot2で推定したパラメータの分布をプロットしたいので、最終的にdata.frameで返ってくるようにする

.c = data.frameとして、data.frameのまま結合してもパラメータの名前がついたdata.frameが返って来なかったため、rbindで結果を結合したあとにdata.frameに変換している(これが最適な方法なのかについては…)

library(sde)
library(dplyr)
library(pforeach)

est <- pforeach(i = 1:1000, .c = rbind)({
  sde.sim(X0 = r0, N = n, delta = delta, 
          theta = c(kappa * mu, kappa, sigma), model = "CIR") %>%
  CIRMLE(initparam, delta)
}) %>% as.data.frame %>% filter(convergence != 1)

foreachの少し込み入った使い方に面倒くささを感じていたが、pforeachを使えば簡単、簡潔に記述できる

replicateを使う場合には、データを発生させて、パラメータを推定してという一連の処理を書いた関数を書いておいて、その関数を繰り返すということをしなければならないが、pforeachではループの中にザッと処理を書いておくだけで繰り返し推定が行える。返り値もあるので、その後一度に加工して扱える。素晴らしい

推定したパラメータの分布をプロットする

library(ggplot2)

ggplot(est, aes(x = kappa)) + geom_histogram(binwidth = 0.1)
ggplot(est, aes(x = mu)) + geom_histogram(binwidth = 0.005)
ggplot(est, aes(x = sigma)) + geom_histogram(binwidth = 0.0001)

平均回帰パラメータkappaの分布のみここでは確認する。推定された値の平均を計算すると真値よりも大きい値となる。そして、プロットしても上方に推定するバイアスが生じることが確認できる(TRUE = 0.2)

f:id:masaqol:20150831010738p:plain