S3 本書のコード一覧
このページでは，本書で利用するRコードがダウンロード可能である．コードは章別に整理してあり，右のメニューから該当ページにアクセスできる．各コードは章ごとに独立して実行できるが，各章の中では順番に実行してほしい．
各コードの右隅にあるアイコン
をクリックすればクリップボードにコピーできるので，それを自身の実行環境にペーストして各コードを実行するのが便利である．なお，本書のコードは，本節の末尾に記す実行環境で動作確認を行っている．
なお，以下で掲載する全てのコードは，本書のGitHubリポジトリから一括してダウンロード可能である．
また，本書の付録的な位置づけとして，以下のRコードと一対一で対応するPythonコードを別のGitHubリポジトリで公開している．Pythonを用いた実証会計・ファイナンスに興味のある読者は併せて参照してほしい．
S3.1 第3章で利用するコード
3.1 Rの基本的な機能
3.1.1 スカラー変数の定義
# ch03_01: 無リスク金利の値を格納する変数Rを作成
R <- 0.1
print(R)
## [1] 0.13.1.2 ベクトル変数の定義
# ch03_02: 2次元の無リスク金利ベクトルを作成
R <- c(0.1, 0.2)
## [1] 0.1 0.2
# ch03_03: 11次元の無リスク金利ベクトルを作成
R <- seq(0.1, 0.2, 0.01)
## [1] 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.20
R <- seq(from = 0.1, to = 0.2, length = 11)
## [1] 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.20
# ch03_04: ベクトルの要素へのアクセス
# 第2要素へのアクセス
R[2]
i <- 2
R[i]
## [1] 0.11
# 第2-第5要素へのアクセス
R[2:5]
## [1] 0.11 0.12 0.13 0.14
# 第2-第11要素へのアクセス
R[2:length(R)] # length(R)はRの次元11を返す
R[-1]
## [1] 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.20
# ch03_05: 1年後の確実な100万円の現在価値
PV <- 100 / (1 + R)
## [1] 90.90909 90.09009 89.28571 88.49558 87.71930 86.95652 86.20690 85.47009 84.74576 84.03361 83.33333
# ch03_06: 次元が一致するベクトル同士の割り算
CF <- seq(from = 110, to = 120, length = 11) # CFを11次元の等差数列として定義
PV <- CF / (1 + R)
## [1] 100 100 100 100 100 100 100 100 100 100 1003.1.3 基本パッケージによる作図
# ch03_07: 割引率と現在価値の関係を描画
PV <- 100 / (1 + R)
plot(PV)
# ch03_08: 割引率と現在価値の関係を表した図の見栄えを改善
PV <- 100 / (1 + R)
plot(x = R,
y = PV,
xlab = "Discount Rate",
ylab = "Present Value",
type = "l",
main = "Figure: Effect of Discount Rate on Present Value")
# ch03_09: 割引率と現在価値の関係を表した図をPNG形式で出力
png("PV.png",
pointsize = 12,
width = 20,
height = 10,
units = "cm",
res = 200)
plot(x = R,
y = PV,
xlab = "Discount Rate",
ylab = "Present Value",
type = "l",
main = "Figure: Effect of Discount Rate on Present Value")
dev.off()
# ch03_10: 3行2列の行列の作成
x <- matrix(1:6, nrow = 3, ncol = 2)
## [,1] [,2]
## [1,] 1 4
## [2,] 2 5
## [3,] 3 6
# ch03_11: 文字を要素とする行列の作成
y <- matrix(c("a", "b", "c", "d"), nrow = 2, ncol = 2)
## [,1] [,2]
## [1,] "a" "c"
## [2,] "b" "d"
# ch03_12: リストの作成
x <- list(10:11, "ABC")
x[[1]]
## [1] 10 11
x[[2]]
## [1] "ABC"
# ch03_13: データフレームの作成
ID <- c(1001, 1002, 1003)
name <- c("Firm A", "Firm B", "Firm C")
ROE <- c(0.08, 0.12, 0.15)
x <- data.frame(ID, name, ROE)
## ID name ROE
## 1 1001 Firm A 0.08
## 2 1002 Firm B 0.12
## 3 1003 Firm C 0.153.2 for文の使い方
3.2.1 NPVの計算例
# ch03_14: 割引率を0.1とする場合のプロジェクトXのNPVの計算
R <- 0.1
NPV <- -100 + 50 / (1 + R) + 50 / (1 + R)^2 + 50 / (1 + R)^3
## [1] 24.34263.2.2 for文を用いたNPVの計算
# ch03_15: for文を用いたプロジェクトXのNPVの計算 (1)
R <- 0.1
NPV <- -100
for (i in 1:3) {
NPV <- NPV + 50 / (1 + R)^i
}
print(NPV)
# ch03_16: for文を用いたプロジェクトXのNPVの計算 (2)
R <- 0.1
NPV <- -100
for (i in 1:3) {
print(NPV) # ここにprint()関数を挿入
NPV <- NPV + 50 / (1 + R)^i
}
print(NPV)
## [1] -100
## [1] -54.54545
## [1] -13.22314
## [1] 24.3426
# ch03_17: NPVに基づく投資判断
if (NPV >= 0) {
print("プロジェクトを実行！")
}3.2.3 NPVと割引率の関係の可視化
# ch03_18: 異なる割引率に対してプロジェクトXのNPVを計算した上で折れ線グラフに可視化
R <- seq(0.1, 0.2, length = 11) # 11次元のRベクトルを用意
NPV <- rep(NA, length(R)) # 対応するNPVをNAで初期化（NAは欠損値の意味）
for (i in 1:length(R)) { # 外側のfor文で特定の割引率を固定
NPV[i] <- -100
for (j in 1:3) { # 内側のfor文で各年度の将来キャッシュフローの現在価値を累積
NPV[i] <- NPV[i] + 50 / (1 + R[i])^j
}
}
plot(x = R,
y = NPV,
xlab = "Discount Rate",
ylab = "Present Value",
type = "l",
main = "Figure: Discount Rate and NPV")3.3 独自関数の定義の仕方
3.3.1 現在価値を計算する独自関数の定義
# ch03_20: 任意の割引率に対して1年後に受け取る100万円の現在価値を計算する独自関数
calculate_PV <- function(R) {
PV <- 100 / (1 + R)
return(PV)
}
calculate_PV(0.1)
## [1] 90.90909
# ch03_21: よりシンプルにcalculate_PV()を定義した場合
calculate_PV <- function(R) 100 / (1 + R)
calculate_PV(0.1)
## [1] 90.909093.3.2 現在価値を計算する関数の機能拡張
# ch03_22: 二つの割引率に対して現在価値を計算できるようcalculate_PV()を拡張
calculate_PV <- function(R1, R2 = 0) {
PV1 <- 100 / (1 + R1)
PV2 <- 100 / (1 + R2)
c(PV1, PV2)
}
calculate_PV(0.1, 0.2)
## [1] 90.90909 83.33333
calculate_PV(0.1)
## [1] 90.90909 100.00000
# ch03_23: 任意の割引率とキャッシュフローに対して現在価値を計算できるようcalculate_PV()を拡張
calculate_PV <- function(R, CF) {
N <- length(CF) # キャッシュフローの発生する回数をNとして定義
PV <- CF[1]
for (i in 2:N) {
PV <- PV + CF[i] / (1 + R)^(i - 1)
}
return(PV)
}
calculate_PV(0.1, c(-100, 40, 50, 60))
## [1] 22.76484
# ch03_24: コメントが詳細で分かりやすいコードの例
# 2021年1月1日作成
# 現在価値の計算
# 任意の割引率に対して現在価値を計算するcalculate_PV()を定義
calculate_PV <- function(R) 100 / (1 + R)
# 割引率を0.1とする例
calculate_PV(0.1)3.4 演習: IRRの計算
3.4.2 多項式の数値解
3.4.3 実数解の抽出
# ch03_28: Im()関数を用いて虚部のみ抽出
Im(solutions)
## [1] 5.714615e-01 -5.714615e-01 -7.815970e-16
# ch03_29: 等価演算子==を用いて実数解かどうかを判定
Im(solutions) == 0
## [1] FALSE FALSE FALSE3.4.4 IRRを計算する独自関数の定義
# ch03_32: 任意のCFに対してIRRを計算するcalculate_IRR()関数を定義
calculate_IRR <- function(CF) {
N <- length(CF) # キャッシュフローの発生する回数をNとして定義
solutions <- polyroot(CF[N:1]) # CFベクトルの順序を逆転させた上でpolyroot()関数に代入
error_tolerance <- 1e-10
is_real <- (abs(Im(solutions)) < error_tolerance)
Re(solutions[is_real]) - 1
}
# ch03_33: calculate_IRR()関数の実行例
CF <- c(-100, 40, 50, 60)
calculate_IRR(CF)
# [1] 0.2164779
# ch03_34: calculate_IRR()関数の完成版（エラーメッセージ付き）
calculate_IRR <- function(CF) {
N <- length(CF)
solutions <- polyroot(CF[N:1])
error_tolerance <- 1e-10
is_real <- (abs(Im(solutions)) < error_tolerance)
if (sum(is_real) != 1) stop("IRR法の適用を再考！") # !=は「左辺と右辺が等しくなければ」の意味
Re(solutions[is_real]) - 1
}
# ch03_35: エラーメッセージが表示される例
CF <- c(-100, 100, 120, -120)
result <- calculate_IRR(CF)
## calculate_IRR(CF) でエラー: IRR法の適用を再考！3.5 データフレーム入門
3.5.1 CSVファイルの読み込み
# ch03_36: CSVファイルの読み込み
daily_stock_return <- read.csv("ch03_daily_stock_return.csv") 3.5.2 データフレームに対する基本的な操作
# ch03_37: データフレームの冒頭を確認
head(daily_stock_return)
# ch03_38: データフレームの行数を確認
nrow(daily_stock_return)
## [1] 21
# ch03_39: データフレームの内部構造を確認
str(daily_stock_return)
## 'data.frame': 21 obs. of 3 variables:
## $ date : chr "2020-04-01" "2020-04-02" "2020-04-03" "2020-04-06" ...
## $ firm1: num -0.03948 0.00598 0.05579 0.04193 -0.02019 ...
## $ firm2: num 0.07696 -0.00725 -0.0173 0.00217 0.07555 ...
# ch03_40: 日次リターンの平均値を計算
mean(daily_stock_return$firm1)
## [1] 0.01416117
mean(daily_stock_return$firm2)
## [1] 0.008025261
# ch03_41: 日次リターンの標準偏差と相関を計算
sd(daily_stock_return$firm1) # sdはstandard deviation（標準偏差）の略
## [1] 0.0538396
cor(daily_stock_return$firm1, daily_stock_return$firm2) # corはcorrelation（相関）の略
## [1] 0.2334892
# ch03_42: 最も日次リターンが低い日付を抽出
worst_day_ID <- which.min(daily_stock_return$firm1) # which.min()関数で行番号を取得
## [1] 20
daily_stock_return$date[worst_day_ID] # 取得した行番号を代入して日付を取得
## [1] "2020-04-28"
# ch03_43: 日次リターンがプラスである日付のみを抽出
daily_stock_return[daily_stock_return$firm1 > 0, ]$date # 企業1のみで条件付け
# 2次元目の要素を空白とすると，全ての列が選択される
# daily_stock_return$date[daily_stock_return$firm1 > 0]と書いても同じ
## [1] "2020-04-02" "2020-04-03" "2020-04-06" "2020-04-09" "2020-04-10" "2020-04-13" "2020-04-15" "2020-04-17" "2020-04-20" "2020-04-22" "2020-04-27" "2020-04-30"
daily_stock_return[daily_stock_return$firm1 > 0 & daily_stock_return$firm2 > 0, ]$date # 企業1と2の両方で条件付け
## [1] "2020-04-06" "2020-04-09" "2020-04-10" "2020-04-13" "2020-04-15" "2020-04-17" "2020-04-20"
# ch03_44: データフレームにおける列の作成・書き換え・削除
daily_stock_return$difference <- daily_stock_return$firm1 - daily_stock_return$firm2 # 両銘柄の差分を計算して新しい列に保存
daily_stock_return$firm1 <- 1 + daily_stock_return$firm1 # 企業1のネットリターンをグロスリターンに変換
daily_stock_return$firm2 <- NULL # 企業2のデータを削除
# ch03_45: データフレームをCSVファイルに出力
write.csv(daily_stock_return, "ch03_output.csv")3.6 ファクター型と日付型
3.6.1 ファクター型入門
# ch03_46: industry列を文字型として定義
firm_ID <- c(1, 2, 3)
name <- c("Firm A", "Firm B", "Firm C")
industry <- c("Machinery", "Chemicals", "Machinery")
firm_data <- data.frame(firm_ID, name, industry) 3.6.2 日付型入門
# ch03_48: date列を文字型として定義
date <- c("2021/4/1", "2021/4/2", "2021/4/5")
stock_return <- c(0.02, -0.01, -0.02)
stock_return_data <- data.frame(date, stock_return)
str(stock_return_data)
## 'data.frame': 3 obs. of 2 variables:
## $ date : chr "2021/4/1" "2021/4/2" "2021/4/5"
## $ stock_return: num 0.02 -0.01 -0.02
# ch03_49: date列を日付型に変換
stock_return_data$date <- as.Date(stock_return_data$date)
str(stock_return_data)
## 'data.frame': 3 obs. of 2 variables:
## $ date : Date, format: "2021-04-01" "2021-04-02" "2021-04-05"
## $ stock_return: num 0.02 -0.01 -0.02
# ch03_50: 日付型の変数に対する足し算の例
stock_return_data$date <- stock_return_data$date + 1 # 文字型の変数に足し算するとエラーになる
head(stock_return_data)
## date stock_return
## 1 2021-04-02 0.02
## 2 2021-04-03 -0.01
## 3 2021-04-06 -0.02 3.7 外部パッケージの使い方
3.7.1 外部パッケージのインストール
# ch03_51: tidyverseのインストール
install.packages("tidyverse")
# ch03_52: dplyrに含まれるnear()関数
dplyr::near(c(1, 2, 1 + 1e-16), 1) # tidyverse::near()だとエラーになる
## [1] TRUE FALSE TRUE
# ch03_53: tidyverseの読み込み
library("tidyverse")
near(c(1, 2, 1 + 1e-16), 1) # dplyr::を省略できる
## [1] TRUE FALSE TRUE3.7.2 作図データの準備
# ch03_54: 異なる初期投資額に対して内部収益率をそれぞれ計算
initial_cost <- seq(80, 100, by = 1) # 21次元のinitial_costベクトルを用意
IRR <- rep(NA, length(initial_cost)) # 対応するIRRをNAで初期化
for (i in 1:length(initial_cost)) {
IRR[i] <- calculate_IRR(c(-initial_cost[i], 100)) # 各initial_costに対応してIRRを計算
}
figure_data <- data.frame(initial_cost, IRR) # 計算結果をデータフレームに格納
head(figure_data)3.7.3 ggplot2による作図
# ch03_55: ggplot2を用いて初期投資額と内部収益率の関係を描画
ggplot(data = figure_data) + geom_line(mapping = aes(x = initial_cost, y = IRR)) # ggplot()関数で描画オブジェクトを作成し, geom_line()関数で折れ線グラフを描画
# ch03_56: 折れ線グラフの上に点グラフを追加
ggplot(data = figure_data) +
geom_line(mapping = aes(x = initial_cost, y = IRR)) +
geom_point(mapping = aes(x = initial_cost, y = IRR)) # geom_point()関数で点グラフを描画
# ch03_57: ggplot2で作成した図の見栄えを改善
ggplot(data = figure_data) +
geom_line(mapping = aes(x = initial_cost, y = IRR)) +
geom_point(mapping = aes(x = initial_cost, y = IRR)) +
labs(x = "Initial Cost", y = "IRR", title = "Initial Cost and IRR") + # 両軸の変数名やタイトルを設定
theme_classic() # 背景などグラフ全体の体裁（テーマ）を設定
# ch03_58: 初期投資額が100の点を強調
ggplot(data = figure_data) +
geom_line(mapping = aes(x = initial_cost, y = IRR)) +
geom_point(mapping = aes(x = initial_cost, y = IRR)) +
labs(x = "Initial Cost", y = "IRR", title = "Initial Cost and IRR") +
theme_classic() +
annotate(geom = "text", x = 99, y = 0.05, label = "Initial Cost 100") + # 位置を指定して文字列を追加
annotate(geom = "segment", x = 100, xend = 100, y = 0.04, yend = 0.01, color = "black", linewidth = 0.5, arrow = arrow(length = unit(0.3, "cm"))) # 始点や終点などを指定して矢印を追加
ggsave("IRR.pdf", width = 20, height = 10, units = "cm") # 完成した図をPDF型式で保存S3.2 第4章で利用するコード
4.2 Rを利用した財務データ分析入門
4.2.2 財務データの読み込み
# ch04_02: CSVファイルの読み込み
financial_data <- read_csv("ch04_financial_data.csv")
nrow(financial_data) # 行数の確認
## [1] 7920
head(financial_data) # 冒頭N行を確認するにはhead(financial_data, N)とする4.3 探索的データ分析
4.3.1 データセットの概要確認
# ch04_05: 要約統計量の表示
summary(financial_data)
# ch04_06: year列の固有要素を抽出
unique(financial_data$year)
## [1] 2015 2016 2017 2018 2019 20204.3.2 欠損データの処理
# ch04_08: 各行の欠損の有無を判定
head(complete.cases(financial_data)) # head()関数を用いて冒頭6行の結果のみ表示
## [1] FALSE TRUE TRUE TRUE TRUE TRUE
# ch04_09: 欠損のない行の総数をカウント
sum(complete.cases(financial_data)) # TRUE/FALSEを1/0に変換して足し合わせる
## [1] 79194.4 データの抽出とヒストグラムによる可視化
4.4.1 条件に合うデータの抽出方法
# ch04_11: 2015年度のデータのみ抽出 (1)
financial_data_2015 <- financial_data[financial_data$year == 2015, ] # 行番号に論理値ベクトルを代入
# ch04_12: 2015年度のデータのみ抽出 (2)
financial_data_2015 <- subset(financial_data, year == 2015) # subset()関数の第二引数に条件式を代入
# ch04_13: 2015年度のデータのみ抽出 (3)
financial_data_2015 <- financial_data %>% filter(year == 2015) # tidyverseのdplyrを用いた書き方
# ch04_14: 2015年度のデータのみ抽出 (4)
financial_data_2015 <- filter(financial_data, year == 2015) # パイプ演算子%>%は第一引数への代入と同じ意味4.4.2 ヒストグラムによる売上高の可視化
# ch04_15: 2015年度の売上高をヒストグラムに図示
ggplot(financial_data_2015) + # 引数名のdataは省略可能
geom_histogram(aes(x = sales)) + # ヒストグラムを描くにはgeom_histogram()関数を用いる
scale_y_continuous(expand = c(0,0)) + # x軸とヒストグラムの間のスペースを消す
theme_classic() # グラフの全体的な体裁を設定
# ch04_16: 桁区切りのカンマをx軸の目盛りで表示
install.packages("scales") # scalesパッケージをインストール
ggplot(financial_data_2015) +
geom_histogram(aes(x = sales)) +
scale_x_continuous(labels = scales::label_comma()) + # scalesパッケージのlabel_comma()関数を利用
scale_y_continuous(expand = c(0,0)) +
theme_classic()
# ch04_17: 売上高の自然対数を取ってヒストグラムで可視化
ggplot(financial_data_2015) +
geom_histogram(aes(x = log(sales))) + # 事前にデータを加工せずともaes()関数内で自然対数が取れる
scale_y_continuous(expand = c(0,0)) +
theme_classic()4.5 データの集計と折れ線グラフによる可視化
4.5.1 for文を用いた集計
# ch04_19: 各年度の上場企業数のカウント (1)
year_set <- sort(unique(financial_data$year)) # year列の固有要素を抽出して昇順に並び替え
N_firms_by_year <- rep(NA, length(year_set)) # 結果を保存するための空ベクトルを準備
for (i in seq_along(year_set)) {
year_i <- year_set[i] # i番目の年度を抽出
N_firms_by_year[i] <- sum(financial_data$year == year_i) # i番目の年度のデータをカウント
}
print(N_firms_by_year) # 結果の表示
## [1] 1265 1293 1319 1323 1356 13634.5.3 dplyrを用いた集計
# ch04_22: 各年度の上場企業数のカウント (3)
N_firms_by_year <- financial_data %>%
group_by(year) %>% # group_by()関数を用いてyearごとにグループ化
summarize(N_firms = n()) # 各グループごとにデータ数をカウント
## # A tibble: 6 x 2
## year N_firms
## <dbl> <int>
## 1 2015 1265
## 2 2016 1293
## 3 2017 1319
## 4 2018 1323
## 5 2019 1356
## 6 2020 13634.5.4 折れ線グラフによる上場企業数の可視化
# ch04_25: 上場企業数の推移を折れ線グラフで可視化
ggplot(N_firms_by_year) +
geom_line(aes(x = year, y = N_firms)) + # 折れ線グラフを描くにはgeom_line()関数を用いる
labs(x = "Year", y = "Number of Firms") + # 両軸のラベルを設定
theme_classic() # グラフ全体の体裁を設定 4.6 変数の作成とヒストグラムによる可視化
4.6.1 新しい系列の追加
# ch04_26: 株主資本の計算 (1)
financial_data$BE <- (financial_data$OA - financial_data$OL) - (financial_data$FO - financial_data$FA) # 新しくBE列を定義4.6.2 ラグの取り方
4.6.3 ROEの計算とヒストグラムによる可視化
# ch04_30: ROEの計算とヒストグラムによる可視化
financial_data <- financial_data %>%
mutate(ROE = X / lagged_BE) # mutate(financial_data, ROE = X / lagged_BE)と書いても同じ
ggplot(financial_data) +
geom_histogram(aes(x = ROE)) +
scale_x_continuous(limits = c(-0.1, 0.5)) + # 表示するx軸の範囲を(-0.1, 0.5)に限定
scale_y_continuous(expand = c(0, 0)) +
theme_classic()4.7 グループごとの集計とランク付け
4.7.1 産業ごとのROE平均値と棒グラフによる可視化
# ch04_31: 産業別のROE平均値を計算
financial_data %>%
group_by(industry_ID) %>% # industry_IDごとにグループ化
summarize(mean_ROE = mean(ROE, na.rm = TRUE)) # 欠損値を除いてROEの平均値を計算
# ch04_32: 産業別のROE平均値を棒グラフで可視化
financial_data %>%
group_by(industry_ID) %>%
summarize(mean_ROE = mean(ROE, na.rm = TRUE)) %>%
ggplot() + # パイプ演算子%>%を用いて, 産業別のROE平均値をggplot()関数に引き渡す
geom_col(aes(x = industry_ID, y = mean_ROE)) + # 棒グラフを描くにはgeom_col()関数を用いる
labs(x = "Industry ID", y = "Industry Average ROE") +
scale_y_continuous(expand = c(0, 0)) +
theme_classic()4.7.2 産業内でのROEのランク付け
# ch04_34: 最終年度に関してfirm_ID・industry_ID・ROEの各列を抽出
financial_data %>%
filter(year == 2020) %>%
select(firm_ID, industry_ID, ROE) # 特定の列のみを抽出するにはselect()関数を用いる
# ch04_35: 産業別に各企業をROEで順序付け (1)
financial_data %>%
filter(year == 2020) %>%
select(firm_ID, industry_ID, ROE) %>%
group_by(industry_ID) %>% # industry_IDでグループ化
mutate(ROE_rank = rank(ROE)) %>% # 順序付けにはrank()関数を用いる
ungroup()
# ch04_36: 産業別に各企業をROEで順序付け (2)
financial_data %>%
filter(year == 2020) %>%
select(firm_ID, industry_ID, ROE) %>%
group_by(industry_ID) %>%
mutate(ROE_rank = rank(desc(ROE))) %>% # 降順で順序付けするにはdesc()関数を用いる
ungroup()
# ch04_37: 産業別に各企業をROEで順序付け (3)
financial_data %>%
filter(year == 2020) %>%
select(firm_ID, industry_ID, ROE) %>%
group_by(industry_ID) %>%
mutate(ROE_rank = rank(desc(ROE))) %>%
ungroup() %>%
arrange(industry_ID, ROE_rank) # 元データを特定の変数の値に基づき並び替えるにはarrange()関数を用いる
# ch04_38: 各産業でROEが最も高い企業のみを抽出 (1)
financial_data %>%
filter(year == 2020) %>%
select(firm_ID, industry_ID, ROE) %>%
group_by(industry_ID) %>%
mutate(ROE_rank = rank(desc(ROE))) %>%
ungroup() %>%
filter(ROE_rank == 1) %>% # ROEが最も高い企業のみ抽出
select(industry_ID, firm_ID, ROE) # ROE_rank列を省略した上で，industry_IDを左端に
# ch04_39: 各産業でROEが最も高い企業のみを抽出 (2)
ROE_rank_data <- financial_data %>% # 中間変数としてROE_rank_dataを定義
filter(year == 2020) %>%
select(firm_ID, industry_ID, ROE) %>%
group_by(industry_ID) %>%
mutate(ROE_rank = rank(desc(ROE))) %>%
ungroup()
ROE_rank_data %>%
filter(ROE_rank == 1) %>%
select(industry_ID, firm_ID, ROE)4.8 上級デュポン・モデルによるROEの分析とその可視化
4.8.1 分析用のデータセットの作成
# ch04_40: 上級デュポン・モデルによるROEの分析
financial_data_DuPont <- financial_data %>%
group_by(firm_ID) %>%
mutate(NOA = OA - OL,
RNOA = OX / lag(NOA),
PM = OX / sales,
ATO = sales / lag(NOA),
NFO = FO - FA,
lagged_FLEV = lag(NFO) / lagged_BE,
NBC = NFE / lag(NFO),
ROE_DuPont = RNOA + lagged_FLEV * (RNOA - NBC)) %>% # 上級デュポン・モデルに基づき計算したROE
ungroup()
# ch04_41: ROEの計算方法の比較
all.equal(financial_data_DuPont$ROE, financial_data_DuPont$ROE_DuPont)
## [1] "Mean relative difference: 4.396878e-06"4.8.2 箱ひげ図による産業別比較
# ch04_42: 箱ひげ図を用いたPMの可視化
financial_data_DuPont %>%
filter(year == 2020,
industry_ID %in% 2:6) %>% # industry_IDが2から6のデータを抽出
ggplot() +
geom_boxplot(aes(x = industry_ID, y = PM)) + # 箱ひげ図を描くにはgeom_boxplot()関数を用いる
labs(x = "Industry ID") +
theme_classic()4.8.3 散布図による産業別比較
# ch04_43: ATOとPMのトレードオフの可視化
financial_data_DuPont %>%
group_by(industry_ID) %>%
summarize(industry_median_PM = median(PM, na.rm = TRUE), # 欠損値を除いた上でPMとATOの産業別中央値を計算
industry_median_ATO = median(ATO, na.rm = TRUE)) %>%
ggplot() +
geom_point(aes(x = industry_median_ATO, y = industry_median_PM)) +
labs(x = "Industry Median ATO", y = "Industry Median PM") +
theme_classic()
# ch04_44: ATOとPMが反比例するグラフを追加
median_RNOA <- median(financial_data_DuPont$RNOA, na.rm = TRUE) # 欠損値を除く全データに関してRNOAの中央値を計算
financial_data_DuPont %>%
group_by(industry_ID) %>%
summarize(industry_median_PM = median(PM, na.rm = TRUE),
industry_median_ATO = median(ATO, na.rm = TRUE)) %>%
ggplot() +
geom_point(aes(x = industry_median_ATO, y = industry_median_PM)) +
labs(x = "Industry Median ATO", y = "Industry Median PM") +
theme_classic() +
stat_function(fun = function(x) median_RNOA / x, linetype = "longdash") # グラフに関数を書き込むにはstat_function()関数を用いる
# ch04_45: データの保存
write_csv(financial_data, "ch04_output.csv") # 基本パッケージのwrite.csv()関数に対応S3.3 第5章で利用するコード
5.1 Rを利用した株式データ分析入門
5.1.2 株価データのダウンロードと読み込み
# ch05_02: 株式データの目視
print(stock_data)5.2 時価総額とリターンの計算
5.2.1 時価総額の計算とヒストグラムによる可視化
# ch05_06: 時価総額の分布をヒストグラムで可視化
library(scales)
ggplot(stock_data) +
geom_histogram(aes(x = ME)) +
labs(x = "Market Equity", y = "Count") +
scale_x_continuous(limits = c(0, quantile(stock_data$ME, 0.95)), # x軸の上限を95%点に設定
labels = label_comma(scale = 1e-6)) + # 単位を100万円にした上で桁区切りのカンマを追加
scale_y_continuous(expand = c(0, 0)) +
theme_classic()5.2.2 トータル・リターンと超過リターンの計算
# ch05_07: 前月の株価を追加
stock_data %>%
group_by(firm_ID) %>% # firm_IDに関してグループ化
mutate(lagged_stock_price = lag(stock_price)) %>%
ungroup()5.2.3 株式データの探索的データ分析
# ch05_10: 月次リターンの要約統計量を確認
summary(stock_data$R)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -0.3803 -0.0406 0.0103 0.0159 0.0648 0.5150 1515
# ch05_11: 月次リターンの標準偏差と分散を計算
sd(stock_data$R, na.rm = TRUE) # na.rm = TRUEを忘れると計算結果はNAに
## [1] 0.09113098
var(stock_data$R, na.rm = TRUE)
## [1] 0.008304855
# ch05_12: 月次リターンの歪度を計算
skewness <- function(x) (1 / length(x)) * sum(((x - mean(x)) / sd(x))^3) # 歪度を計算する関数の定義
skewness(na.omit(stock_data$R)) # 欠損値を削除した上で代入
## [1] 0.5065672
# ch05_13: 月次リターンの尖度を計算
kurtosis <- function(x) (1 / length(x)) * sum(((x - mean(x)) / sd(x))^4) # 尖度を計算する関数の定義
kurtosis(na.omit(stock_data$R)) # 欠損値を削除した上で代入
## [1] 4.268045
# ch05_14: 月次リターンの分布をヒストグラムで可視化
ggplot(stock_data) +
geom_histogram(aes(x = R)) +
labs(x = "Monthly Stock Return", y = "Count") +
scale_y_continuous(expand = c(0, 0)) +
theme_classic()5.3 リターンの累積
5.4 株式データと財務データを組み合わせた分析
5.4.1 二つのデータセットの結合
# ch05_17: 財務データと株式データの行数を確認
financial_data <- read_csv("ch04_output.csv")
nrow(financial_data) # 年次財務データの行数
## [1] 7919
nrow(annual_stock_data) # 年次リターン・データの行数
## [1] 7920
nrow(stock_data) # 月次リターン・データの行数
## [1] 95040
# ch05_18: 年次リターン・データと財務データの結合 (1)
annual_data <- annual_stock_data %>%
full_join(financial_data, by = c("firm_ID", "year")) # firm_IDとyearのペアをキーとして設定
# ch05_19: 年次リターン・データと財務データの結合 (2)
annual_data <- annual_stock_data %>%
full_join(financial_data) # キーを省略した場合, 列名が同じ変数がキーに
## Joining, by = c("firm_ID", "year")
# ch05_20: 月次リターン・データと財務データの結合
monthly_data <- stock_data %>%
full_join(financial_data, by = c("firm_ID", "year"))
# ch05_21: full_join()関数による欠損値の処理
A <- tibble(firm_ID = c(1, 2), stock_price = c(120, 500)) # データセットの作成
B <- tibble(firm_ID = c(1, 3), DPS = c(5, 10))
full_join(A, B, by = "firm_ID") # A %>% full_join(B, by = "firm_ID")と書いても同じ
# ch05_22: inner_join()関数による欠損値の処理
inner_join(A, B, by = "firm_ID")
# ch05_23: left_join()関数，及びright_join()関数による欠損値の処理
left_join(A, B, by = "firm_ID")
right_join(A, B, by = "firm_ID")5.4.2 結合後データの探索的データ分析とバブルチャートによる可視化
# ch05_24: 年度末の時価総額を年次データに追加
annual_data <- stock_data %>%
filter(month == 12) %>% # 12月のデータのみを抽出
select(year, firm_ID, ME) %>% # 追加したい列のみ選択
full_join(annual_data, ., by = c("year", "firm_ID")) %>% # 年次データと結合
mutate(ME = ME / 1e6) # 時価総額の単位を百万円に
# ch05_25: 売上高・当期純利益・時価総額をバブルチャートで可視化
annual_data %>%
filter(year == 2015,
firm_ID %in% 2:20, # firm_IDが2から20のデータを抽出
X > 0) %>% # 対数を取るため当期純利益（X）が正のデータのみ抽出
ggplot() +
geom_point(aes(x = log(sales), y = log(X), size = ME), alpha = 0.4) + # バブルチャートを描くにはsize引数を指定
scale_size(range = c(1, 20), name = "Market Equity") + # range引数でバブルの最小・最大面積を指定
scale_x_continuous(limits = c(8, 14)) + # 両軸の範囲を指定
scale_y_continuous(limits = c(2, 11)) +
theme_classic()5.5 統計的推論入門
5.5.1 リターン・データに関する仮定
# ch05_26: firm_IDが1の企業の月次超過リターンを折れ線グラフで可視化
stock_data %>%
filter(firm_ID == 1) %>% # firm_IDが1の企業のデータのみ抽出
ggplot() +
geom_line(aes(x = month_ID, y = Re)) + # x軸にmonth_ID, y軸に月次超過リターンを表示
labs(x = "Month ID", y = "Firm 1's Monthly Excess Return") +
scale_x_continuous(expand = c(0.01, 0)) + # 折れ線グラフとy軸の間の空白を指定
theme_classic()
# ch05_27: firm_IDが1の企業の月次超過リターンをヒストグラムで可視化
stock_data %>%
filter(firm_ID == 1) %>%
ggplot() +
geom_histogram(aes(x = Re)) +
labs(x = "Firm 1's Monthly Excess Return", y = "Count") +
scale_y_continuous(expand = c(0, 0)) +
theme_classic()5.5.3 \(t\)値の計算
# ch05_29: 月次超過リターンの期待値に関するt検定 (1)
Re_firm_ID_1 <- stock_data %>%
filter(firm_ID == 1) %>% # firm_IDが1の企業のみ抽出
select(Re) %>% # 月次超過リターンのみ抽出
drop_na() %>% # 欠損値を削除
unlist() # データフレームからベクトルに変換
mu0 <- 0 # 帰無仮説を期待値0と設定
n <- length(Re_firm_ID_1) # 標本サイズ
t_value <- (mean(Re_firm_ID_1) - mu0) / sqrt(var(Re_firm_ID_1) / n) # 定義に従ってt値を計算
## [1] 2.1212965.5.4 統計的検定の考え方
# ch05_30: 月次超過リターンの期待値に関するt検定 (2)
t.test(Re_firm_ID_1) # t.test()関数で第二引数以降を省略5.6 線形回帰入門
5.6.1 線形回帰の概要
# ch05_31: 簿価時価比率と株式リターンの関係を可視化
lm_sample_data <- annual_data %>%
group_by(firm_ID) %>%
mutate(lagged_BEME = lagged_BE / lag(ME)) %>% # 企業ごとに前年度の簿価時価比率を計算
ungroup() %>%
filter(year == 2016, # firm_IDが1から10までの企業の2016年のデータを抽出
firm_ID <= 10) %>%
select(firm_ID, year, Re, lagged_BEME) %>% # 必要な列のみ抽出
drop_na() # 欠損データを削除
ggplot(lm_sample_data) +
geom_point(aes(x = lagged_BEME, y = Re)) + # x軸に簿価時価比率, y軸に超過リターン
labs(x = "BE/ME at the End of Year t", y = "Excess Return for Year t + 1") +
theme_classic()
# ch05_32: 簿価時価比率と株式リターンの散布図に回帰直線を追加
ggplot(lm_sample_data) +
geom_point(aes(x = lagged_BEME, y = Re)) +
geom_smooth(aes(x = lagged_BEME, y = Re), method = "lm", se = FALSE, color = "black") + # 回帰直線を追加するにはgeom_smooth()関数を用いる
labs(x = "BE/ME at the End of Year t", y = "Excess Return for Year t + 1") +
theme_classic()5.6.2 最小二乗推定量
# ch05_33: lm()関数を用いた線形回帰
lm_results <- lm(Re ~ lagged_BEME, data = lm_sample_data) # ~の左に従属変数, 右に独立変数を記す
names(lm_results)
## [1] "coefficients" "residuals" "effects" "rank" "fitted.values" "assign" "qr" "df.residual" "xlevels" "call" "terms" "model"
# ch05_34: lm()関数で計算した回帰係数
print(lm_results$coefficients)
## (Intercept) lagged_BEME
## 0.05513385 0.07552921 5.6.3 線形回帰モデルの解釈
# ch05_35: broomパッケージのtidy()関数で係数の推定値に関する結果を確認
install.packages("broom")
library(broom)
tidy(lm_results)
# ch05_36: broomパッケージのglance()関数でモデル全体の当てはまりを確認
glance(lm_results)S3.4 第6章で利用するコード
6.1 ファクター構築の準備
6.1.1 市場ポートフォリオの構築
# ch06_01: tidyverseとデータの読み込み
library(tidyverse)
library(broom)
monthly_data <- read_csv("ch05_output1.csv")
annual_data <- read_csv("ch05_output2.csv")
# ch06_02: 前年度の時価総額を追加
annual_data <- annual_data %>%
group_by(firm_ID) %>% # firm_IDでグループ化
mutate(lagged_ME = lag(ME)) %>%
ungroup()
# ch06_03: 市場ポートフォリオ作成のための保有比率の計算
annual_data <- annual_data %>%
group_by(year) %>% # yearでグループ化
mutate(w_M = lagged_ME / sum(lagged_ME, na.rm = TRUE)) %>% # 分子は当該銘柄の時価総額, 分母は各銘柄の時価総額の合計
ungroup()
# ch06_04: 保有比率の欠損値を0で置き換え
annual_data <- annual_data %>%
mutate(w_M = replace(w_M, year >= 2016 & is.na(w_M), 0)) # 2016年以降の欠損値は0で置き換え
# ch06_05: 各年度の保有比率の合計が1になっていることを確認
annual_data %>%
group_by(year) %>%
summarize(weight_sum = sum(w_M)) # あえてsum(w_M, na.rm = TRUE)としない
# ch06_06: 月次データに保有比率のデータを追加
monthly_data <- annual_data %>%
select(year, firm_ID, w_M) %>% # 追加に必要な情報のみ抽出
full_join(monthly_data, by = c("year", "firm_ID")) %>% # yearとfirm_IDのペアをキーに結合
select(-w_M, w_M) # w_Mを最終列に移動
# ch06_07: 市場ポートフォリオの月次リターンを計算
factor_data <- monthly_data %>%
filter(month_ID >= 13) %>% # 2016年以降のデータを抽出
group_by(month_ID) %>% # month_IDでグループ化
summarize(R_F = R_F[1], # 当該月の無リスク金利はどの銘柄でも共通なので最初の値を抽出
R_M = sum(w_M * R, na.rm = TRUE)) %>% # 各銘柄の月次リターンの加重平均を計算
mutate(R_Me = R_M - R_F) # 月次超過リターンを計算
# ch06_08: 市場ポートフォリオの月次超過リターンをヒストグラムで可視化
ggplot(factor_data) +
geom_histogram(aes(x = R_Me)) +
labs(x = "Monthly Excess Return of Market Portfolio", y = "Count") +
scale_y_continuous(expand = c(0, 0)) +
theme_classic()
# ch06_09: 市場ポートフォリオの累積リターンの可視化 (1)
factor_data %>%
mutate(gross_R_M = 1 + R_M, # ネット・リターンをグロス・リターンに換算
cumulative_gross_R_M = cumprod(gross_R_M)) %>% # グロス・リターンを累積
ggplot() +
geom_line(aes(x = month_ID, y = cumulative_gross_R_M)) +
labs(x = "Month ID", y = "Cumulative Gross Return") +
theme_classic()
# ch06_10: 市場ポートフォリオの累積リターンの可視化 (2)
factor_data %>%
mutate(gross_R_M = 1 + R_M,
cumulative_gross_R_M = cumprod(gross_R_M)) %>%
select(month_ID, cumulative_gross_R_M) %>%
rbind(c(12, 1), .) %>% # 折れ線グラフの始点を追加
ggplot() +
geom_line(aes(x = month_ID, y = cumulative_gross_R_M)) +
geom_hline(yintercept = 1, linetype = "dotted") + # 元本の水準を点線で図示
labs(x = "Month ID", y = "Cumulative Gross Return") +
scale_x_continuous(expand = c(0, 0)) +
theme_classic()6.1.2 ポートフォリオ・ソート
# ch06_11: 前年度末の時価総額に基づくポートフォリオ・ソート (1)
annual_data <- annual_data %>%
group_by(year) %>%
mutate(ME_rank10 = as.factor(ntile(lagged_ME, 10))) %>% # ntile()関数を用いて十個のグループに分類
ungroup()
# ch06_12: 各ポートフォリオに属する企業数を確認
annual_data %>%
select(year, firm_ID, ME_rank10) %>%
drop_na() %>% # 欠損行を削除
group_by(year, ME_rank10) %>% # yearとME_rank10でグループ化
summarize(N = n()) %>% # 各ポートフォリオの企業数をカウント
ungroup()
# ch06_13: 前年度末の時価総額に基づくポートフォリオ・ソート (2)
ME_sorted_portfolio <- annual_data %>%
select(year, firm_ID, ME_rank10) %>% # 年次データから追加したい情報を抽出
full_join(monthly_data, by = c("year", "firm_ID")) %>% # yearとfirm_IDをキーに月次データと結合
drop_na() %>% # 欠損行を削除
group_by(month_ID, ME_rank10) %>% # month_IDとME_rank10に関してグループ化
summarize(Re = mean(Re)) %>% # 各グループで月次超過リターンの平均値を計算
ungroup()
# ch06_14: 各ポートフォリオの平均超過リターンを可視化
ME_cross_sectional_return <- ME_sorted_portfolio %>%
group_by(ME_rank10) %>% # ME_rank10に関してグループ化
summarize(mean_Re = mean(Re)) # 月次超過リターンの平均値を計算
ggplot(ME_cross_sectional_return) +
geom_col(aes(x = ME_rank10, y = mean_Re)) + # 棒グラフを描くにはgeom_col()関数を用いる
labs(x = "ME Rank", y = "Mean Monthly Excess Return") +
scale_y_continuous(expand = c(0, 0)) +
theme_classic()
# ch06_15a: 簿価時価比率に基づくポートフォリオ・ソート（等加重の場合）
annual_data <- annual_data %>%
mutate(lagged_BEME = lagged_BE / lagged_ME) %>%
group_by(year) %>%
mutate(BEME_rank10 = as.factor(ntile(lagged_BEME, 10))) %>% # 簿価時価比率に基づいて十個のグループに分類
ungroup()
BEME_sorted_portfolio <- annual_data %>%
select(year, firm_ID, BEME_rank10, lagged_ME) %>%
full_join(monthly_data, by = c("year", "firm_ID")) %>%
drop_na() %>%
group_by(month_ID, BEME_rank10) %>%
summarize(Re = mean(Re)) %>% # 月次超過リターンの平均値を計算
ungroup()
group_by(BEME_sorted_portfolio, BEME_rank10) %>%
summarize(mean_Re = mean(Re)) %>%
ggplot() +
geom_col(aes(x = BEME_rank10, y = mean_Re)) +
geom_hline(yintercept = 0) + # y = 0の直線を追加
labs(x = "BE/ME Rank", y = "Mean Monthly Excess Return") +
scale_y_continuous(limits = c(-0.005, 0.02)) +
theme_classic()
# ch06_15b: 簿価時価比率に基づくポートフォリオ・ソート（時価総額加重の場合）
# 中盤で保有比率wと月次超過リターンReを計算している箇所を除けば, ch06_15aと全く同じ
annual_data <- annual_data %>%
mutate(lagged_BEME = lagged_BE / lagged_ME) %>%
group_by(year) %>%
mutate(BEME_rank10 = as.factor(ntile(lagged_BEME, 10))) %>%
ungroup()
BEME_sorted_portfolio <- annual_data %>%
select(year, firm_ID, BEME_rank10, lagged_ME) %>%
full_join(monthly_data, by = c("year", "firm_ID")) %>%
drop_na() %>%
group_by(month_ID, BEME_rank10) %>%
mutate(w = lagged_ME / sum(lagged_ME)) %>% # 各ポートフォリオで保有比率を計算
summarize(Re = sum(w * Re)) %>% # 時価総額加重の月次超過リターンを計算
ungroup()
group_by(BEME_sorted_portfolio, BEME_rank10) %>%
summarize(mean_Re = mean(Re)) %>%
ggplot() +
geom_col(aes(x = BEME_rank10, y = mean_Re)) +
geom_hline(yintercept = 0) +
labs(x = "BE/ME Rank", y = "Mean Monthly Excess Return") +
scale_y_continuous(limits = c(-0.005, 0.02)) +
theme_classic()6.2 CAPMの実証的な検証
6.2.2 時系列回帰
# ch06_16: 市場ポートフォリオの超過リターンを追加
ME_sorted_portfolio <- factor_data %>%
select(-R_F) %>% # 無リスク金利は重複するので結合前に削除
full_join(ME_sorted_portfolio, by = "month_ID") %>% # month_IDをキーに
select(-R_Me, R_Me) # R_Meを最終列へ移動
# ch06_17: 時系列回帰 (1)
ME_sorted_portfolio %>%
filter(ME_rank10 == 1) %>% # 時価総額が最小のポートフォリオを抽出
lm(Re ~ R_Me, data = .) %>% # .を使ってlm()関数の第二引数にデータを代入
tidy() # 線形回帰の結果をtidy()関数でデータフレームに変換
# ch06_18: 時系列回帰 (2)
ME_sorted_portfolio %>%
filter(ME_rank10 == 1) %>%
ggplot(aes(x = R_Me, y = Re)) + # aes()関数はggplot()関数の中にも代入可能
geom_point() + # geom_point()関数と次のgeom_smooth()関数で共通のaes()関数を受け取る
geom_smooth(method = "lm", color = "black") +
labs(x = "Excess Return of Market Portfolio", y = "Excess Return of Small Size Portfolio") +
theme_classic()6.2.3 ポートフォリオごとの回帰
# ch06_19: CAPMの実証的な検証 (1)
CAPM_results <- list(NA) # 推定結果を保存するために空のリストを準備
for(i in 1:10){
CAPM_results[[i]] <- ME_sorted_portfolio %>%
filter(ME_rank10 == i) %>%
lm(Re ~ R_Me, data = .) %>%
tidy() %>%
mutate(ME_rank10 = i) %>% # 推定対象のポートフォリオ名を保存
select(ME_rank10, everything()) # ME_rank10を第一列に移動
}
# ch06_20: CAPMの実証的な検証 (2)
CAPM_results <- do.call(rbind, CAPM_results) # do.call()関数を用いて複数のデータフレームから構成されるリストを一つのデータフレームに統合
# ch06_21: グループごとの線形回帰 (1) lapply()関数を使う場合
ME_sorted_portfolio_splitted <- split(ME_sorted_portfolio, ME_sorted_portfolio$ME_rank10) # 元データをME_rank10の値に応じて十個のデータフレームに分割
estimate_CAPM <- function(return_data) { # リターン・データを受け取り, CAPMの推定結果をデータフレームで返す関数を準備
lm_results <- lm(Re ~ R_Me, data = return_data)
tidied_lm_results <- tidy(lm_results)
}
CAPM_results_by_lapply <- lapply(ME_sorted_portfolio_splitted, estimate_CAPM) # lapply()関数は第一引数にリスト, 第二引数に関数を取る
# lapply()関数の返り値はリストなので，一つのデータフレームにまとめたい場合はdo.call()関数を用いる
# ch06_22: グループごとの線形回帰 (2) map()関数を使う場合
ME_sorted_portfolio %>%
group_by(ME_rank10) %>%
nest() %>% # nest()関数を用いて各グループを要素とするメタ・データフレームを作成
mutate(CAPM_regression = map(data, ~ lm(Re ~ R_Me, data = .)), # map()関数を用いて各グループを線形回帰
CAPM_summary = map(CAPM_regression, tidy)) %>% # tidy()関数を用いて線形回帰の結果を整理
select(-c(data, CAPM_regression)) %>% # 線形回帰の結果のみを抽出
unnest(cols = CAPM_summary) %>% # nest()関数による畳み込みを解除
ungroup()6.2.4 CAPMアルファ
# ch06_23: CAPMアルファの可視化
CAPM_results %>%
filter(term == "(Intercept)") %>% # 定数項に関する推定結果のみを抽出
mutate(ME_rank10 = as.factor(ME_rank10)) %>% # ME_rank10を整数型からファクター型に
ggplot() +
geom_col(aes(x = ME_rank10, y = estimate)) + # 横軸をME_rank10, 縦軸をCAPM_alphaとする棒グラフ
geom_hline(yintercept = 0) +
labs(x = "ME Rank", y = "CAPM alpha") +
scale_y_continuous(limits = c(-0.003, 0.013)) +
theme_classic()
# ch06_24: CAPMアルファの統計的な有意性を評価
CAPM_results %>%
filter(term == "(Intercept)") %>% # 定数項に関する推定結果のみを抽出
rename(CAPM_alpha = estimate, p_value = p.value) %>% # 列名を変更
mutate(significance = cut(p_value,
breaks = c(0, 0.01, 0.05, 0.1, 1),
labels = c("***", "**", "*", ""),
include.lowest = TRUE)) %>% # 統計的に有意な結果を*で強調
select(ME_rank10, CAPM_alpha, p_value, significance) # 出力したい列を指定
# ch06_25: 証券市場線の推定
ME_cross_sectional_return <- CAPM_results %>%
filter(term == "R_Me") %>% # R_Meの係数に関する推定結果のみを抽出
rename(CAPM_beta = estimate) %>% # 推定値をestimateからCAPM_betaに名称変更
select(ME_rank10, CAPM_beta) %>%
mutate(ME_rank10 = as.factor(ME_rank10)) %>% # ME_rank10を整数型からファクター型に
full_join(ME_cross_sectional_return, ., by = "ME_rank10") # 超過リターンのデータと結合
mean_R_Me <- mean(factor_data$R_Me) # 市場ポートフォリオの実現超過リターンにより証券市場線の傾きを推定
ggplot(ME_cross_sectional_return) +
geom_point(aes(x = CAPM_beta, y = mean_Re)) +
geom_abline(intercept = 0, slope = mean_R_Me) +
labs(x = "Market beta", y = "Mean Excess Return") +
scale_x_continuous(limits = c(0, 1.2), expand = c(0, 0)) +
scale_y_continuous(limits = c(0, 0.02)) +
theme_classic()6.3 Fama-Frenchの3ファクター・モデル
6.3.3 銘柄のランク付け
# ch06_26: 前年度の時価総額に基づくランク付け
annual_data <- annual_data %>%
mutate(lagged_ME = replace(lagged_ME, is.na(lagged_BEME), NA)) %>% # lagged_BEMEが欠損している場合は欠損扱いに
group_by(year) %>%
mutate(ME_rank2 = as.factor(ntile(lagged_ME, 2))) %>%
ungroup()
# ch06_27: 簿価時価比率に基づくランク付け (1)
annual_data %>%
group_by(year) %>%
mutate(BEME_percent_rank = percent_rank(lagged_BEME)) %>% # 年度ごとに簿価時価比率のパーセンタイル順位を計算
ungroup()
# ch06_28: 簿価時価比率に基づくランク付け (2)
annual_data <- annual_data %>%
group_by(year) %>%
mutate(BEME_percent_rank = percent_rank(lagged_BEME)) %>% # 年度ごとに簿価時価比率のパーセンタイル順位を計算
ungroup() %>%
mutate(BEME_rank3 = cut(BEME_percent_rank,
breaks = c(0, 0.3, 0.7, 1),
labels = c(1, 2, 3),
include.lowest = TRUE)) # BEME_percent_rankの値に応じて1から3までBEME_rank3の値を定義6.3.4 時価総額とBE/MEに基づくポートフォリオ・ソート
# ch06_29: 6 Size-BE/MEポートフォリオへの分類 (1)
annual_data <- annual_data %>%
mutate(FF_portfolio_type = interaction(ME_rank2, BEME_rank3)) # ME_rank2とBEME_rank3の組合せで, ファクター型の変数FF_portfolio_typeを定義
# ch06_30: 6 Size-BE/MEポートフォリオへの分類 (2)
annual_data <- annual_data %>%
mutate(FF_portfolio_type = fct_recode(FF_portfolio_type,
SL = "1.1",
BL = "2.1",
SN = "1.2",
BN = "2.2",
SH = "1.3",
BH = "2.3")) # FF_portfolio_typeの水準を変更
# ch06_31: 6 Size-BE/MEポートフォリオへの分類 (3)
annual_data %>%
group_by(ME_rank2, BEME_rank3) %>% # ME_rank2とBEME_rank3のペアでグループ化
summarize(FF_portfolio_type = FF_portfolio_type[1],
mean_BEME = mean(lagged_BEME),
mean_ME = mean(lagged_ME),
mean_N_stocks = n() / length(unique(year))) %>%
ungroup() %>%
drop_na() # 欠損データを削除
# ch06_32: 6 Size-BE/MEポートフォリオの構築 (1)
annual_data <- annual_data %>%
group_by(year, FF_portfolio_type) %>% # yearとFF_portfolio_typeのペアでグループ化
mutate(w = lagged_ME / sum(lagged_ME, na.rm = TRUE)) %>% # 各ポートフォリオ内で時価総額加重の保有比率を計算
ungroup()
# ch06_33: 6 Size-BE/MEポートフォリオの構築 (2)
FF_portfolio <- annual_data %>%
select(year, firm_ID, FF_portfolio_type, ME_rank2, BEME_rank3, w) %>%
full_join(monthly_data, by = c("year", "firm_ID")) %>% # 今までに準備したデータと月次データを結合
group_by(month_ID, FF_portfolio_type) %>% # month_IDとFF_portfolio_typeでグループ化
summarize(ME_rank2 = ME_rank2[1],
BEME_rank3 = BEME_rank3[1],
R = sum(w * R, na.rm = TRUE), # 各ポートフォリオの月次リターンを計算
R_F = R_F[1]) %>%
ungroup() %>%
drop_na() # 欠損データを削除
# ch06_34: 6 Size-BE/MEポートフォリオのリターンの可視化 (1)
FF_portfolio_mean_return <- FF_portfolio %>%
mutate(Re = R - R_F) %>%
group_by(FF_portfolio_type) %>% # FF_portfolio_typeでグループ化
summarize(ME_rank2 = ME_rank2[1],
BEME_rank3 = BEME_rank3[1],
mean_Re = mean(Re)) # 各ポートフォリオの超過リターンの平均値を計算
ggplot(FF_portfolio_mean_return) +
geom_col(aes(x = BEME_rank3, y = mean_Re, fill = ME_rank2), position = "dodge") + # x軸をBEME_rank3, y軸をmean_Reに, ME_rank2のサブグループで色分け
scale_fill_grey() + # 棒グラフの色をモノトーンに
labs(x = "BE/ME Rank", y = "Mean Monthly Excess Return", fill = "ME Rank") +
scale_y_continuous(expand = c(0, 0)) +
theme_classic()
# ch06_35: 6 Size-BE/MEポートフォリオのリターンの可視化 (2)
ggplot(FF_portfolio_mean_return) +
geom_col(aes(x = BEME_rank3, y = mean_Re, fill = ME_rank2), position = "dodge") +
scale_fill_grey() +
geom_text(aes(x = BEME_rank3, y = mean_Re, group = ME_rank2, label = FF_portfolio_type), # (x, y)座標を指定して各ポートフォリオの名前をグラフに挿入
vjust = -0.5, # 棒グラフが重ならないよう文字ラベルを上にずらす
position = position_dodge(width = 0.9)) + # ME_rank2のサブグループで文字ラベルが左右にずれるよう調整
labs(x = "BE/ME Rank", y = "Mean Monthly Excess Return", fill = "ME Rank") +
scale_y_continuous(expand = c(0, 0), limits = c(0, 0.015)) + # 文字ラベルがはみ出ないようy軸の範囲を指定
theme_classic()
# ch06_36: 6 Size-BE/MEポートフォリオのリターンの可視化 (3)
initial_point <- tibble(month_ID = c(12, 12), # 累積リターンの起点を定義
cumulative_gross_R = c(1, 1),
FF_portfolio_type = c("BL", "SH"))
FF_portfolio_cumulative_return <- FF_portfolio %>%
group_by(FF_portfolio_type) %>% # FF_portfolio_typeでグループ化
mutate(cumulative_gross_R = cumprod(1 + R)) %>% # グロス・リターンを累積
ungroup() %>%
filter(FF_portfolio_type %in% c("BL", "SH")) %>%
select(month_ID, cumulative_gross_R, FF_portfolio_type) %>%
rbind(initial_point, .) # initial_pointを第一行に挿入
ggplot(FF_portfolio_cumulative_return) +
geom_line(aes(x = month_ID, y = cumulative_gross_R, linetype = FF_portfolio_type)) +
scale_linetype_manual(values = c("longdash", "solid")) +
geom_hline(yintercept = 1, linetype = "dotted") +
labs(x = "Month ID", y = "Cumulative Gross Return", linetype = "") +
scale_x_continuous(expand = c(0, 0)) +
theme_classic()6.3.5 ファクター・リターンの計算
# ch06_37: SMBとHMLの構築 (1)
FF_portfolio <- FF_portfolio %>%
pivot_wider(id_cols = month_ID, names_from = FF_portfolio_type, values_from = R) # FF_portfolio_typeの値に基づく列を作成し, 縦長から横長のデータに変換
# ch06_38: SMBとHMLの構築 (2)
factor_data <- FF_portfolio %>%
mutate(SMB = (SH + SN + SL) / 3 - (BH + BN + BL) / 3, # SMBとHMLを計算
HML = (SH + BH) / 2 - (SL + BL) / 2) %>%
select(month_ID, SMB, HML) %>%
full_join(factor_data, by = "month_ID") %>% # 3ファクターの実現値をfactor_dataに集約
select(-c("SMB", "HML"), c("SMB", "HML")) # SMBとHMLを最後列に移動6.3.6 FF3アルファ
# ch06_39: FF3モデルの推定
ME_sorted_portfolio <- ME_sorted_portfolio %>%
select(-c(R_Me, R_M)) %>%
full_join(factor_data, by = "month_ID") # 3ファクターの実現値をME_sorted_portfolioに追加
FF3_results <- list(NA) # 推定結果を保存するために空のリストを準備
for(i in 1:10) {
FF3_results[[i]] <- ME_sorted_portfolio %>%
filter(ME_rank10 == i) %>%
lm(Re ~ R_Me + SMB + HML, data = .) %>% # 3ファクターの実現値を独立変数として重回帰
tidy() %>%
mutate(ME_rank10 = i) %>% # 推定対象のポートフォリオ名を保存
select(ME_rank10, everything()) # ME_rank10を第一列に移動
}
FF3_results <- do.call(rbind, FF3_results) # do.call()関数を用いて複数のデータフレームから構成されるリストを一つのデータフレームに統合
# ch06_40: FF3アルファの可視化
FF3_results %>%
filter(term == "(Intercept)") %>% # 定数項に関する推定結果のみを抽出
mutate(ME_rank10 = as.factor(ME_rank10)) %>% # ME_rank10を整数型からファクター型に
ggplot() +
geom_col(aes(x = ME_rank10, y = estimate)) + # 横軸をME_rank10, 縦軸をFF3_alphaとする棒グラフ
geom_hline(yintercept = 0) +
labs(x = "ME Rank", y = "FF3 alpha") +
scale_y_continuous(limits = c(-0.003, 0.013)) +
theme_classic()
# ch06_41: FF3アルファの統計的な有意性を評価
FF3_results %>%
filter(term == "(Intercept)") %>% # 定数項に関する推定結果のみを抽出
rename(FF3_alpha = estimate, p_value = p.value) %>% # 列名を変更
mutate(significance = cut(p_value,
breaks = c(0, 0.01, 0.05, 0.1, 1),
labels = c("***", "**", "*", ""),
include.lowest = TRUE)) %>% # 統計的に有意な結果を*で強調
select(ME_rank10, FF3_alpha, p_value, significance) # 出力したい列を指定
# ch06_42: データの保存
write_csv(factor_data, "ch06_output.csv")S3.5 第7章で利用するコード
7.1 資本コストの推定
7.1.2 FF3モデルによる株式資本コストの推定
# ch07_01: 外部パッケージの読み込みとデータの準備
library(tidyverse)
library(broom)
monthly_data <- read_csv("ch05_output1.csv")
annual_data <- read_csv("ch05_output2.csv")
factor_data <- read_csv("ch06_output.csv")
monthly_data <- monthly_data %>%
group_by(firm_ID) %>% # firm_IDでグループ化
mutate(is_public_2020 = (max(month_ID) == 72), # 2020年の年末時点で上場しているかフラグ付け
N_observations = n()) %>% # 各firm_IDのデータ数をカウント
ungroup() %>%
filter(is_public_2020 == TRUE, # 2020年の年末時点で上場しているfirm_IDを抽出
N_observations >= 36) %>% # 36ヶ月以上のデータがあるfirm_IDを抽出
select(-N_observations) # N_observations列を削除
annual_data <- monthly_data %>% # monthly_dataに含まれるyearとfirm_IDのペアを抽出
select(year, firm_ID) %>%
unique() %>%
left_join(annual_data, by = c("year", "firm_ID")) # annual_dataとmonthly_dataでyearとfirm_IDのペアを整合的に
# ch07_02: FF3モデルの推定 (1)
firm_ID_set <- unique(monthly_data$firm_ID) # firm_IDの固有要素を抽出
N_firms <- length(firm_ID_set) # firm_ID_setの要素数をカウント
monthly_data <- monthly_data %>%
select(-R_F) %>% # R_Fはmonthly_dataとfactor_dataで重複するので削除
full_join(factor_data, by = "month_ID")
FF3_results <- list(NA) # 推定結果を保存するために空のリストを準備
for (i in 1:N_firms) {
FF3_results[[i]] <- monthly_data %>%
filter(firm_ID == firm_ID_set[i]) %>%
lm(Re ~ 0 + R_Me + SMB + HML, data = .) %>% # 冒頭に"0 +"を追加して定数項をゼロに
tidy() %>%
mutate(firm_ID = firm_ID_set[i]) %>% # 推定対象のfirm_IDを保存
select(firm_ID, everything())
}
FF3_results <- do.call(rbind, FF3_results) # do.call()関数を用いて複数のデータフレームから構成されるリストを一つのデータフレームに統合
# ch07_03: FF3モデルの推定 (2)
FF3_loadings <- FF3_results %>%
pivot_wider(id_cols = "firm_ID", names_from = "term", values_from = "estimate") %>% # term列の値に応じて新しく列を定義
rename(beta_M = R_Me, beta_SMB = SMB, beta_HML = HML) # 列名を変更
# ch07_04: ファクター・リスクプレミアムの計算
factor_risk_premium <- apply(factor_data[ , 4:6], 2, mean) * 12 # 12を掛けて年次プレミアムに換算
## R_Me SMB HML
## 0.04800514 0.07858358 0.08328208
# ch07_05: Fama-Frenchの3ファクター・モデルに基づいて株式資本コストを推定
expected_R_FF3 <- as.matrix(FF3_loadings[ , 2:4]) %*% factor_risk_premium # Fama-Frenchの3ファクター・モデルに基づいて期待リターンを計算
FF3_cost_of_capital <- tibble(firm_ID = FF3_loadings$firm_ID, cost_of_capital = as.vector(expected_R_FF3)) # 行列からデータフレームに変換してfirm_ID列を追加
# ch07_06: 株式資本コストの推定値をヒストグラムで可視化
ggplot(FF3_cost_of_capital) +
geom_histogram(aes(x = cost_of_capital)) +
labs(x = "FF3 Cost of Capital", y = "Count") +
scale_y_discrete(expand = c(0, 0)) +
theme_classic()
# ch07_07: 時価総額別に株式資本コストの推定値を可視化
annual_data_2020 <- annual_data %>%
group_by(firm_ID) %>% # 企業ごとに時価総額を計算
mutate(lagged_ME = lag(ME)) %>%
ungroup() %>%
filter(year == 2020) %>% # 2020年のデータのみ抽出
mutate(ME_rank2 = as.factor(ntile(lagged_ME, 2)), # 時価総額に基づき二分割
ME_rank2 = fct_recode(ME_rank2,
Small = "1",
Large = "2")) %>%
full_join(FF3_cost_of_capital, by = "firm_ID") %>% # 株式資本コストの推定値を追加
drop_na() # 欠損値を削除
ggplot(annual_data_2020) +
geom_density(aes(x = cost_of_capital, group = ME_rank2, linetype = ME_rank2)) + # 密度関数を描くにはgeom_density()関数を用いる
labs(x = "FF3 Cost of Capital", y = "Estimated Density", linetype = "Firm Type") +
scale_y_continuous(expand = c(0, 0), breaks = NULL) + # break引数で目盛り線を非表示に指定
theme_classic()7.1.3 株式資本コストの応用:潜在配当成長率の推定
# ch07_08: 配当利回りの計算
dividend_yield_2020 <- monthly_data %>%
filter(year == 2020) %>%
group_by(firm_ID) %>% # firm_IDでグループ化
summarize(annual_dividend = sum(DPS * shares_outstanding), # 各月の配当支払総額を合計
latest_ME = ME[which.max(month_ID)]) %>% # 最新の時価総額を保存
mutate(dividend_yield = annual_dividend / latest_ME) %>% # 配当利回りの実績値を計算
filter(dividend_yield > 0) # 無配の企業を落とす7.2 平均分散ポートフォリオの構築
7.2.2 分散共分散行列の推定
# ch07_10: ファクターの分散共分散行列の推定
Sigma_FF3 <- 12 * cov(factor_data[ , 4:6]) # 12を掛けて年次データに換算
## R_Me SMB HML
## R_Me 0.020685969 -0.003901607 0.001529890
## SMB -0.003901607 0.006268377 -0.001886351
## HML 0.001529890 -0.001886351 0.003915638
# ch07_11: 投資対象企業の選定
N_portfolio_firms <- 100 # 投資対象の企業を100社に限定
investment_universe <- annual_data %>%
filter(year == 2020) %>% # 2020年度のデータを抽出
filter(rank(desc(ME)) <= N_portfolio_firms) %>% # その中で時価総額が上位100社の銘柄を抽出
select(firm_ID) %>%
unlist() # データフレームからベクトルに変換
names(investment_universe) <- NULL # firm_IDという列名を消去
## [1] 33 101 155 …（中略）… 1454 1473 1499
# ch07_12: 誤差項の分散共分散行列の推定
Sigma_epsilon <- monthly_data %>%
filter(firm_ID %in% investment_universe) %>% # 投資対象に含まれるデータのみを抽出
left_join(FF3_loadings, by = "firm_ID") %>%
mutate(R_FF3 = R_F + beta_M * R_Me + beta_SMB * SMB + beta_HML * HML, # ファクターの実現値からFama-Frenchの3ファクター・モデルに基づくリターンを計算
epsilon = R - R_FF3) %>% # 実際のリターンとの違いから誤差項を推定
group_by(firm_ID) %>%
summarize(epsilon_variance = 12 * var(epsilon, na.rm = TRUE)) %>% # firm_IDごとに誤差項の分散を推定
select(epsilon_variance) %>%
unlist() %>% # データフレームからベクトルに変換
diag() # ベクトルから対角行列を作成
# ch07_13: リターンの分散共分散行列を準備
beta <- FF3_loadings %>%
filter(firm_ID %in% investment_universe) %>% # 投資対象に含まれるデータのみを抽出
select(-firm_ID) %>% # ファクター・ローディングのみを抽出
as.matrix() # データフレームから行列に変換
Sigma <- beta %*% Sigma_FF3 %*% t(beta) + Sigma_epsilon # 分解式に基づいて分散共分散行列を計算
# ch07_14: 期待リターンのベクトルを準備
mu <- FF3_cost_of_capital %>% # 7.1.2節の推定結果をそのまま利用
filter(firm_ID %in% investment_universe) %>% # 投資対象に含まれるデータのみを抽出
select(-firm_ID) %>%
unlist() # データフレームからベクトルに変換
names(mu) <- NULL # cost_of_capitalという列名を消去
## [1] 0.203791678 -0.001750860 0.071846430 …（中略）… 0.088872719 0.153100331 0.3616627587.2.3 平均分散ポートフォリオの計算
# ch07_15: quadprogパッケージのインストール
install.packages("quadprog") # quadprogはquadratic programming（二次計画法）の略
# ch07_16: solve.QP()関数を利用した平均分散ポートフォリオの計算 (1)
library(quadprog) # quadprogの読み込み
target_return <- 0.1 # 目標期待リターンを0.1に設定
Dmat <- Sigma # 分散共分散行列
dvec <- rep(0, N_portfolio_firms) # 対応項無し
Amat <- cbind(mu, rep(1, N_portfolio_firms)) # 目標期待リターン，及び保有比率の係数
bvec <- c(target_return, 1) # 目標期待リターン，及び保有比率の合計値
MV_portfolio <- solve.QP(Dmat, dvec, Amat, bvec, meq = 2) # 等号制約の数をmeq引数で表す
# ch07_17: solve.QP()関数を利用した平均分散ポートフォリオの計算 (2)
str(MV_portfolio) # solve.QP()関数の返り値の構造を確認
# ch07_18: solve.QP()関数を利用した平均分散ポートフォリオの計算 (3)
optimal_weight <- MV_portfolio$solution # 最適保有比率を保存
print(optimal_weight[1:10]) # 最初の10社のみ表示
## [1] 0.0580726666 0.0245486385 0.0121478107 -0.0270912773 0.0123755582 -0.1703321035 0.0118779871 0.0128634656 0.0517985612 -0.0111403267
minimized_risk <- sqrt(2 * MV_portfolio$value) # 平均分散ポートフォリオのリスクを保存して表示
print(minimized_risk)
## [1] 0.083710037.2.4 平均分散フロンティアの描画
# ch07_19: 平均分散ポートフォリオの描画 (1)
target_return <- seq(-0.1, 0.4, length = 100) # -0.1から0.4の範囲を離散化して目標期待リターンを100次元のベクトルで準備
N_points <- length(target_return) # target_returnの次元100をN_pointsと定義
optimal_weight <- matrix(NA, nrow = N_points, ncol = N_portfolio_firms) # 各目標期待リターンに対して最適保有比率を保存する空行列を準備
minimized_risk <- rep(NA, N_points) # 同様に平均分散ポートフォリオのリスクを保存する空ベクトルを準備
for (i in 1:N_points) { # 平均分散ポートフォリオの計算を100回繰り返す
bvec <- c(target_return[i], 1)
MV_portfolio <- solve.QP(Dmat, dvec, Amat, bvec, meq = 2)
optimal_weight[i, ] <- MV_portfolio$solution # 結果の保存
minimized_risk[i] <- sqrt(2 * MV_portfolio$value)
}
print(optimal_weight[N_points, 1:10]) # target_return = 0.4のときの最適保有比率（最初の10社のみ）
## [1] 0.102911702 0.020711220 0.001801566 -0.046204100 0.016718819 -0.444390843 0.023600316 0.021858155 0.096069270 -0.004657669
print(minimized_risk[N_points]) # target_return = 0.4のときの平均分散ポートフォリオのリスク
## [1] 0.2148099
# ch07_20: 平均分散ポートフォリオの描画 (2)
MV_frontier <- tibble(target_return, minimized_risk)
ggplot(MV_frontier) +
geom_line(aes(x = target_return, y = minimized_risk)) + # x軸にtarget_return, y軸にminimized_riskを指定
coord_flip() + # 縦軸と横軸をひっくり返すのにcoord_flip()関数を用いる
labs(x = "Expected Return", y = "Risk") +
theme_classic()S3.6 第8章で利用するコード
8.3 決算発表のイベント・スタディを例にしたRでの実践
8.3.2 グループ分けの方法
# ch08_02: イベントIDを付与した上で予想利益サプライズを計算
event_data <- event_data %>%
mutate(event_ID = 1:nrow(event_data), # イベントIDを付与
forecast_innovation = (earnings_forecast - realized_earnings) / lagged_ME, # 予想利益サプライズを計算
year = as.integer(substr(event_date, 1, 4))) %>% # 日付データから年度を抽出
group_by(year) %>%
mutate(event_strength = as.factor(ntile(forecast_innovation, 5))) %>% # 各年度ごとにイベントの強弱を5段階に区分
ungroup() 8.3.3 相対日次の設定
# ch08_03: パラメータの準備
N_days <- nrow(market_return_data) # データに含まれる日数をカウント
N_firms <- length(unique(return_data$firm_ID)) # 企業の固有数をカウント
N_events <- nrow(event_data) # イベント数をカウント
N1 <- 100 # モデルの推定期間
N2 <- 30 # イベント前の分析期間
N3 <- 30 # イベント後の分析期間
# ch08_04: 日付IDの付与
market_return_data$date_ID <- 1:N_days # 各営業日に日付IDを順番に付与
data_ID_table <- market_return_data %>%
select(date, date_ID) # 各営業日と日付IDを紐付けるデータテーブルを作成
return_data <- return_data %>%
full_join(data_ID_table, by = "date") # リターン・データに日付IDを付与
event_data <- event_data %>% # イベント・データに日付IDを付与
left_join(data_ID_table, by = c("event_date" = "date")) %>% # event_dataとdata_ID_tableでキー列の名前が異なる点に注意
select(event_ID, date_ID, firm_ID, event_strength) %>% # 不要な列を削除
rename(event_date_ID = date_ID) # 後のステップで列名の重複を避けるために列名を変更
# ch08_05: 各データセットを統合して単一のデータフレームに (1)
return_data <- return_data %>%
select(-date) %>% # date列が重複するので削除
full_join(market_return_data, by = "date_ID") # 市場ポートフォリオのリターンと結合
# ch08_06: 各データセットを統合して単一のデータフレームに (2)
full_sample_data <- tibble(event_ID = sort(rep(1:N_events, N1 + N2 + N3 + 1)),
relative_days = rep(-(N1 + N2):N3, N_events)) # イベント・データとリターン・データを結合する上で骨格となるデータフレームを作成
# ch08_07: 各データセットを統合して単一のデータフレームに (3)
full_sample_data <- full_sample_data %>%
full_join(event_data, by = "event_ID") %>% # イベント・データと結合
mutate(date_ID = event_date_ID + relative_days + 1) %>% # event_date_IDとrelative_daysからdate_IDを逆算（取引時間終了後の決算発表を仮定して1日足す）
left_join(return_data, by = c("firm_ID", "date_ID")) %>% # リターン・データと結合
select(event_ID, event_strength, relative_days, R, R_M)8.3.4 異常リターンの算定
# ch08_08: マーケット・モデルの推定 (1)
estimation_window_data <- full_sample_data %>%
filter(relative_days < -N2) # 推定期間のデータのみ抽出
market_model_results <- list(NA) # 推定結果を保存するために空のリストを準備
for (i in 1:N_events) {
lm_results <- estimation_window_data %>%
filter(event_ID == i) %>%
lm(R ~ R_M, data = .) # マーケット・モデルの推定
tidied_lm_results <- lm_results %>%
tidy() %>%
mutate(event_ID = i) %>% # 推定対象のイベントIDを保存
select(event_ID, everything()) # event_IDを第一列に移動
tidied_lm_results$sigma_AR <- glance(lm_results)$sigma # マーケット・モデルの推定により得られた残差の標準誤差を抽出してsigma_ARと保存
market_model_results[[i]] <- tidied_lm_results # リストの一要素として推定結果を保存
}
market_model_results <- do.call(rbind, market_model_results) # do.call()関数を用いて複数のデータフレームから構成されるリストを一つのデータフレームに統合
# ch08_09: マーケット・モデルの推定 (2)
full_sample_data <- market_model_results %>%
pivot_wider(id_cols = c("event_ID", "sigma_AR"), names_from = "term", values_from = "estimate") %>% # term列の値に応じて新しく列を定義
rename(alpha = "(Intercept)", beta = R_M) %>% # 列名を変更
full_join(full_sample_data, by = "event_ID") %>% # full_sample_dataと結合
select(-c("alpha", "beta", "sigma_AR"), c("alpha", "beta", "sigma_AR")) # alpha，beta，sigma_ARを最終列に移動8.3.5 異常リターンの集計と結果の解釈
# ch08_11: イベントの強弱に応じて平均CARの推移を可視化
mean_CAR_by_event_strength <- event_window_data %>%
group_by(relative_days, event_strength) %>%
summarize(mean_CAR = mean(CAR)) # CARの平均値を計算
ggplot(mean_CAR_by_event_strength) +
geom_line(aes(x = relative_days, y = mean_CAR, linetype = event_strength)) + # event_strengthごとに線の種類を変更
scale_linetype_manual(values = c("dotted", "dotdash", "dashed", "longdash", "solid")) +
annotate("rect", xmin = -1, xmax = 0, ymin = -Inf, ymax = Inf, alpha = 0.1) + # イベント発生期間を灰色で強調
labs(x = "Relative Days", y = "Mean CAR", linetype = "Event Strength") +
scale_x_continuous(expand = c(0.02, 0)) + # y軸と折れ線グラフの左端の空白を指定
scale_y_continuous(labels = scales::label_percent()) + # y軸をパーセント表示に
theme_classic()
# ch08_12: ARを利用した統計的検定 (1)
output_table <- event_window_data %>%
filter(event_strength == 5) %>% # イベント強度が最も強いグループを抽出
group_by(relative_days) %>%
summarize(mean_AR = mean(AR),
mean_CAR = mean(CAR),
sigma_mean_AR = sqrt(sum(sigma_AR^2)) / n()) # ARとCARの平均値，及び平均ARの日次標準偏差を推定
# ch08_13: ARを利用した統計的検定 (2)
output_table <- output_table %>%
mutate(t_value = mean_AR / sigma_mean_AR, # t値を計算
p_value = (1 - pnorm(abs(t_value))) * 2) %>% # 対応するp値を計算
mutate(significance = cut(p_value,
breaks = c(0, 0.01, 0.05, 0.1, 1),
labels = c("***", "**", "*", ""),
include.lowest = TRUE)) %>% # 統計的に有意な結果を*で強調
select(relative_days, mean_AR, t_value, p_value, significance, mean_CAR) %>% # 列の順序を変更
mutate(mean_AR = round(mean_AR, 5) * 100,
t_value = round(t_value, 2),
p_value = round(p_value, 2),
mean_CAR = round(mean_CAR, 5) * 100) # 各列で表示する桁数を指定動作確認をした実行環境
install.packages("sessioninfo")
sessioninfo::session_info()
## ─ Session info ─────────────────────────────────────
## setting value
## version R version 4.0.5 (2021-03-31)
## os macOS 12.1
## system x86_64, darwin17.0
## ui RStudio
## language (EN)
## rstudio 2021.09.2+382 Ghost Orchid (desktop)
## ─ Packages ─────────────────────────────────────────
## package * version date (UTC) lib source
## broom * 0.7.12 2022-01-28 [1] CRAN (R 4.0.5)
## dplyr * 1.0.7 2021-06-18 [1] CRAN (R 4.0.2)
## forcats * 0.5.1 2021-01-27 [1] CRAN (R 4.0.2)
## ggplot2 * 3.3.5 2021-06-25 [1] CRAN (R 4.0.2)
## purrr * 0.3.4 2020-04-17 [1] CRAN (R 4.0.2)
## quadprog * 1.5-8 2019-11-20 [1] CRAN (R 4.0.2)
## readr * 2.1.2 2022-01-30 [1] CRAN (R 4.0.5)
## scales * 1.1.1 2020-05-11 [1] CRAN (R 4.0.2)
## stringr * 1.4.0 2019-02-10 [1] CRAN (R 4.0.2)
## tibble * 3.1.6 2021-11-07 [1] CRAN (R 4.0.2)
## tidyr * 1.2.0 2022-02-01 [1] CRAN (R 4.0.5)
## tidyverse * 1.3.1 2021-04-15 [1] CRAN (R 4.0.2)