一般化線形モデル周りの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 Implementation in R | 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"))
その他、色々と付け足していけば、思いどおりのプロットができるでしょう。しかし、本題はデータを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
matplotならばどうかというと、xtsプロットの良い部分は完全に失われます。
> matplot(rate_xts, type = "l")
xtsはzooクラスとしても扱えるのだから、plot.zooで良い感じにプロットされるのではないかと思い、やってみるとmatplotとほぼ変わりません…
> plot.zoo(rate_xts, plot.type = "single")
というところまで試してみて、CRANには上がっていないxtsExtraを使うのが一番良さそうということが分かり、使ってみるとxtsプロットの良い部分もあり、matplotのような複数系列のプロットが簡単にできました。
library(devtools) devtools::install_github("joshuaulrich/xtsExtra") library(xtsExtra) plot(rate_xts, col = 1:15)
ただ、デフォルトで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）
CIRモデルのパラメータ推定
CIRモデルについてもVasicekモデルと同様にシンプルに最尤推定でモデルのパラメータを推定できる関数を用意しておく
# maximum likelihood estimation for CIR model # # rate : interest rate # param : model parameter # delta : time step # initparam : initial parameter # CIRloglik <- function(rate, param, delta) { n <- length(rate) c <- 4 * param[1] / (param[3] ^ 2 * (1 - exp(-param[1] * delta))) df <- 4 * param[1] * param[2] / (param[3] ^ 2) ncp <- c * rate[1:(n - 1)] * exp(-param[1] * delta) sum(dchisq(c * rate[2:n], df = df, ncp = ncp, log = TRUE) + log(c)) } CIRMLE <- function(rate, initparam, delta) { opt <- optim(initparam, CIRloglik, rate = rate, delta = delta, control = list(fnscale = -1)) c(kappa = opt$par[1], mu = opt$par[2], sigma = opt$par[3], loglik = opt$value, convergence = opt$convergence) }
パラメータを設定して、パスを発生させる
library(sde) r0 <- 0.01 n <- 500 delta <- 1/52 kappa <- 0.2 mu <- 0.03 sigma <- 0.08 initparam <- c(kappa, mu, sigma) set.seed(12345) rate <- sde.sim(X0 = r0, N = n, delta = delta, theta = c(kappa * mu, kappa, sigma), model = "CIR")
推定結果を確認する
> CIRMLE(rate, initparam, delta) kappa mu sigma loglik convergence 2.789964e-01 6.105776e-02 8.106622e-02 2.376850e+03 0.000000e+00
初期値の与え方にもよるが、Vasicekモデルと比べて多少パラメータの推定値にバイアスが出ることになります。 特に平均回帰パラメータkappaについては上方に推定する傾向にあります
Rからアメリカ雇用統計データ収集
FXをやっている人なら誰しも気になる経済統計といえば、アメリカの雇用統計。 【米国】雇用統計 - 経済指標 - Yahoo!ファイナンスでも直近の数年間のデータを詳しく見ることができますが、簡単にRから直接データを収集して、分析したいという人はblsAPIパッケージを使ってみるのが良さそうです
こちらmikeasilva/blsAPI · GitHubにもサンプルがありますが、発表時に一番注目される非農業部門雇用者数と失業率のデータを使った例がないようなので、そのデータを取得してから可視化するところまでをメモ
使い方は簡単で、パッケージに唯一含まれているblsAPIという関数に欲しいデータのIDコードを与えてやればデータを取って来れます。その他、list形式でパラメータを指定してやれば、いろいろと自由にデータ収集可能です
library(rjson) library(blsAPI) library(dplyr) payload <- list( 'seriesid' = c('CES0000000001', 'LNS14000000'), 'startyear' = 2007, 'endyear' = 2015) response <- blsAPI(payload) %>% fromJSON()
これだけで、2007年から2015年までの雇用者数と失業率のデータが取って来れます。 responseの中に2系列入っており、年月についても整形してdata.frameに変換してやります
series1 <- response$Results$series[[1]]$data series2 <- response$Results$series[[2]]$data date_str <- 0 value1 <- value2 <- 0 for (i in 1:length(series1)) { date_str[i] <- paste(series1[[i]]$year, sub('M', '', series1[[i]]$period), "01", sep = "-") value1[i] <- series1[[i]]$value value2[i] <- series2[[i]]$value } date <- as.Date(date_str[order(date_str)])[-1] data1 <- diff(as.numeric(value1[order(date_str)])) data2 <- as.numeric(value2[order(date_str)])[-1] df <- data.frame(date, data1, data2)
データを確認します
> head(df) date data1 data2 1 2007-02-01 88 4.5 2 2007-03-01 188 4.4 3 2007-04-01 78 4.5 4 2007-05-01 145 4.4 5 2007-06-01 71 4.6 6 2007-07-01 -34 4.7 > tail(df) date data1 data2 95 2014-12-01 329 5.6 96 2015-01-01 201 5.7 97 2015-02-01 266 5.5 98 2015-03-01 119 5.5 99 2015-04-01 221 5.4 100 2015-05-01 280 5.5
前月比で雇用者数の変化を棒グラフにする
library(ggplot2) library(scales) df %>% ggplot(aes(x = date)) + geom_bar(aes(y = data1), stat = "identity", fill = "blue") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + scale_x_date(labels = date_format("%m/%y"), breaks = date_breaks("6 months")) + scale_y_continuous(breaks = seq(-1000, 1000, by = 200)) + ylab("1-Month Net Change")
失業率もプロットする
df %>% ggplot(aes(x = date)) + geom_line(aes(y = data2), size = 1, colour = "red") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + scale_x_date(labels = date_format("%m/%y"), breaks = date_breaks("6 months")) + scale_y_continuous(breaks = seq(0, 12, by = 1)) + ylab("Unemployment Rate")
Vasicekモデルのパラメータ推定
SMFI5パッケージでもVasicekモデルのパラメータ推定はできるが、参考程度にパラメータを推定する際、シンプルに最尤推定できるものを用意しておいて、簡単に使い回せるようにしておきたい
# maximum likelihood estimation for Vasicek model # # rate : interest rate # param : Vasicek model parameter # delta : time step # initparam : initial parameter # Vasicekloglik <- function(rate, param, delta) { n <- length(rate) m <- param[2] + (rate[1:(n-1)] - param[2]) * exp(-param[1] * delta) v <- param[3]^2 * (1 - exp(-2 * param[1] * delta)) / (2 * param[1]) sum(dnorm(rate[2:n], mean = m, sd = sqrt(v), log = TRUE)) } VasicekMLE <- function(rate, initparam, delta) { opt <- optim(initparam, Vasicekloglik, rate = rate, delta = delta, control = list(fnscale = -1)) c(kappa = opt$par[1], mu = opt$par[2], sigma = opt$par[3], loglik = opt$value, convergence = opt$convergence) }
optimのhessian = TRUEを足したり、２つの関数を１つにしてしまったりということもできる
sdeパッケージを使ってパスを発生させ、パラメータ推定する
library(sde) r0 <- 0.01 n <- 500 delta <- 1/52 kappa <- 0.2 mu <- 0.03 sigma <- 0.01 initparam <- c(kappa, mu, sigma) set.seed(123) rate <- sde.sim(X0 = r0, N = n, delta = delta, theta = c(kappa * mu, kappa, sigma), model = "OU")
結果を確認する
> VasicekMLE(rate, initparam, delta) kappa mu sigma loglik convergence 5.670395e-01 3.114274e-02 9.740629e-03 2.596687e+03 0.000000e+00
最適化の初期値に真の値を入れていることもあり、良い値に推定できています
金利の平均回帰速度kappaはだいぶ真の値よりも上方に推定されていますが、推定バイアスを改善させましょうというような論文がいくつかありますので、この辺をじっくり読むというのも大変面白いです
Rの金利期間構造パッケージ ycinterextra
金利の期間構造関連パッケージを使ってみるシリーズ。今回はycinterextraを調査。このパッケージもこれといった情報がないので、普通にパッケージの説明書を読んで、日本の金利データに適用してみる
パッケージの中身は、CRAN Task View: Empirical Financeにも書かれている通り、Nelson-Siegelモデル等の様々なモデルでイールドカーブの補間と補外を行うことができるパッケージだと言うこと
主な関数はycinterとycextraで、その他ユーティリティ的な関数が用意されている
ycinter(yM = NULL, p = NULL, matsin, matsout, method = c("NS", "SV", "SW", "HCSPL"), typeres = c("rates", "prices")) ycextra(yM = NULL, p = NULL, matsin, matsout, method = c("NS", "SV", "SW"), typeres = c("rates", "prices"), UFR, T_UFR = NULL)
ycinterとycextraともに、ゼロクーポン債価格とゼロクーポンイールドを対象としているので、クーポン債に適用したい場合にはその前に処理しておく必要がある
jgbcme <- read.csv("http://www.mof.go.jp/english/jgbs/reference/interest_rate/jgbcme.csv", skip = 1) JGB_rate <- as.numeric(jgbcme[1, -1]) / 100 maturity <- c(1:10, 15, 20, 25, 30, 40)
Nelson-SiegelモデルとSvenssonモデルをデータに適用していく
library(ycinterextra) NSyc <- ycinter(yM = JGB_rate, matsin = maturity, matsout = 1:40, method = "NS", typeres = "rates") SVyc <- ycinter(yM = JGB_rate, matsin = maturity, matsout = 1:40, method = "SV", typeres = "rates")
求まったパラメータの値は次の通り
> coeffs(NSyc) [1] 0.01847065 -0.01415590 -0.03796077 3.00000000 > coeffs(SVyc) [1] 0.02019232 -0.02168073 0.04297268 -0.07961120 2.01990134 3.00000000
3.00000000と出ていることからも上手く最適化できていない可能性が高そうです…
ycplotを使って、当てはまりの診断を見てみる
> ycplot(NSyc) > ycplot(SVyc)
Nelson-Siegelモデル
Svenssonモデル
残差が本当に正規なのというところは疑問が残るところですが、それなりに補間は上手く行っているようです
（method = "HCSPL"を指定すれば、Hermite cubic splineによる補間もできるのだが、上と同じようにしてもエラーが出てなぜか上手く適用できなかった…）
同じように補外の方も1年から60年までの設定で試してみる
ここで、UFRの値は上のパラメータの値から0.02を設定している
NSyc <- ycextra(yM = JGB_rate, matsin = maturity, matsout = 1:60, method = "NS", typeres = "rates", UFR = 0.02) SVyc <- ycextra(yM = JGB_rate, matsin = maturity, matsout = 1:60, method = "SV", typeres = "rates", UFR = 0.02)
求まったパラメータの値は次の通り
> coeffs(NSyc) [1] 0.02000000 -0.01457810 -0.04355021 3.00000000 > coeffs(SVyc) [1] 0.02000000 -0.02162734 0.03623047 -0.07322875 1.90287193 3.00000000
やはり、上手く最適化できていないようです…
こちらもycplotで確認
> ycplot(NSyc) > ycplot(SVyc)
Nelson-Siegelモデル
Svenssonモデル
少々、簡単でしたがycinterextraパッケージの関数を使って、実際データに適用してみました。今の日本の金利データがこれらのモデルに適用するというのが難しい面もありそうですが、パラメータを求めるあたりはあまり上手くいっていないようです。このパッケージでは、これらのモデル以外にも、Smith-Wilson methodも使えるのでお手軽に試してみたいという時は良さそうです
