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.1

3.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 100

3.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.15

3.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.3426

3.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")
# ch03_19: ベクトル化によるNPVの計算

R <- 0.1
CF <- c(-100, 50, 50, 50)
years <- 0:3
PV_CF <- CF / (1 + R)^years
NPV <- sum(PV_CF)

## [1] 24.3426

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.90909

3.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)
# ch03_25: 関数内でグローバル変数の変更を試みた例

calculate_PV <- function(R) {
  PV <- 100 / (1 + R)
  R <- R + 0.01
  return(PV)
}

R <- 0.1
for (i in 1:11) print(calculate_PV(R))

## [1] 90.90909
## [1] 90.90909
## (同じ計算結果が全部で11回繰り返し表示される)

3.4 演習: IRRの計算

3.4.2 多項式の数値解

# ch03_26: polyroot()関数に基づく3次方程式の数値解

polyroot(c(60, 50, 40, -100))

## [1] -0.4082389+0.5714615i -0.4082389-0.5714615i  1.2164779-0.0000000i
# ch03_27: 数値解の中から実数解を目視で選んでIRRを計算

solutions <- polyroot(c(60, 50, 40, -100))
Y <- Re(solutions[3])
IRR <- Y - 1

## [1] 0.2164779

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 FALSE
# ch03_30: 数値的誤差を許容して実数解かどうかを判定

error_tolerance <- 1e-10

abs(Im(solutions)) # abs()関数は絶対値を返す
## [1] 5.714615e-01 5.714615e-01 7.815970e-16

abs(Im(solutions)) < error_tolerance
## [1] FALSE FALSE  TRUE
# ch03_31: 数値解の中から実数解を機械的に抽出

is_real <- (abs(Im(solutions)) < error_tolerance)
Re(solutions[is_real]) # Re()関数は複素数の実部を返す

## [1] 1.216478

3.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) 
# ch03_47: industry列をファクター型に変換

firm_data$industry <- as.factor(firm_data$industry)
str(firm_data)

## 'data.frame':    3 obs. of  3 variables:
##  $ firm_ID : num  1 2 3
##  $ name    : chr  "Firm A" "Firm B" "Firm C"
##  $ industry: Factor w/ 2 levels "Chemicals","Machinery": 2 1 2

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  TRUE

3.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", size = 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.1 tidyverseパッケージの概要

# ch04_01: tidyverseの読み込み

library(tidyverse)

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)とする
# ch04_03: firm_ID列とindustry_ID列をファクター型に変換

financial_data$firm_ID <- as.factor(financial_data$firm_ID)
financial_data$industry_ID <- as.factor(financial_data$industry_ID)

class(financial_data$firm_ID)
## [1] "factor"

class(financial_data$industry_ID)
## [1] "factor"
# ch04_04: read_csv()関数を用いるとdate列が最初から日付型に

daily_stock_return <- read_csv("ch03_daily_stock_return.csv")
class(daily_stock_return$date)

## [1] "Date"

4.3 探索的データ分析

4.3.1 データセットの概要確認

# ch04_05: 要約統計量の表示

summary(financial_data)
# ch04_06: year列の固有要素を抽出

unique(financial_data$year)

## [1] 2015 2016 2017 2018 2019 2020
# ch04_07: firm_ID,及びindustry_ID列の固有要素数をカウント

length(unique(financial_data$firm_ID))
## [1] 1515

length(unique(financial_data$industry_ID))
## [1] 10

4.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] 7919
# ch04_10: 欠損行を削除

nrow(drop_na(financial_data)) # 欠損行を削除した上で行数をカウント
## [1] 7919

financial_data <- drop_na(financial_data) # 欠損行を削除した上でデータを上書き

4.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()
# ch04_18: パイプ演算子%>%が便利な理由

financial_data$sales %>% 
  log() %>% # この時点でlog(financial_data$sales)を計算
  median()
## [1] 10.6074

median(log(financial_data$sales))
## [1] 10.6074

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 1363
# ch04_20: ベクトルの名前付け

names(N_firms_by_year) <- year_set # 各要素の名前として年度名を代入
print(N_firms_by_year)

## 2015 2016 2017 2018 2019 2020 
## 1265 1293 1319 1323 1356 1363 

4.5.2 tapply()関数を用いた集計

# ch04_21: 各年度の上場企業数のカウント (2)

N_firms_by_year <- tapply(financial_data$firm_ID, financial_data$year, length) # tapply()関数は第一引数に元データ, 第二引数にグループ分けに用いる変数, 第三引数に適用したい関数を取る
print(N_firms_by_year)

## 2015 2016 2017 2018 2019 2020 
## 1265 1293 1319 1323 1356 1363 

4.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    1363
# ch04_23: group_by()関数を用いたグループ化

financial_data %>% group_by(year) # group_by(financial_data, year)と書いても同じ
# ch04_24: summarize()関数を用いた平均値の計算

financial_data %>% summarize(mean_sales = mean(sales)) # summarize(financial_data, mean_sales = mean(sales))と書いても同じ

# A tibble: 1 x 1
## mean_sales
## <dbl>
##   1    1520.647 (本文165頁では左記のように記していますが,正しくは以下の通りです)
##   1    166027.

4.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列を定義
# ch04_27: 株主資本の計算 (2)

financial_data <- financial_data %>%
  mutate(BE = (OA - OL) - (FO - FA)) # mutate(financial_data, BE = (OA - OL) - (FO - FA))と書いても同じ

4.6.2 ラグの取り方

# ch04_28: lag()関数を用いて前期の株主資本を取得 (1)

head(lag(financial_data$BE)) # head()関数で冒頭6行の結果のみ表示

## [1]       NA 10013.82 10426.33 10842.01 11074.95 11593.58
# ch04_29: lag()関数を用いて前期の株主資本を取得 (2)

financial_data <- financial_data %>%
  group_by(firm_ID) %>% # firm_IDごとにグループ化
  mutate(lagged_BE = lag(BE)) %>% # lag()関数をfirm_IDごとに適用
  ungroup() # グループ化を解除

head(financial_data$lagged_BE) # head()関数で冒頭6行のlagged_BEを表示

## [1]       NA 10013.82 10426.33 10842.01 11074.95       NA

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_33: 最終年度のデータのみを抽出

financial_data %>% 
  filter(year == 2020)
# 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_01: tidyverseと株式データの読み込み

library(tidyverse)
stock_data <- read_csv("ch05_stock_data.csv")
# ch05_02: 株式データの目視

print(stock_data)
# ch05_03: 配当支払いや発行済株式数変化の例

stock_data %>% 
  filter(firm_ID == 1 & month_ID %in% 27:30)
# ch05_04: 調整係数が1以外の値を取る例

stock_data %>% 
  filter(firm_ID == 74 & month_ID %in% 29:32)

5.2 時価総額とリターンの計算

5.2.1 時価総額の計算とヒストグラムによる可視化

# ch05_05: 時価総額の追加

stock_data <- stock_data %>% 
  mutate(ME = stock_price * shares_outstanding)
# 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()
# ch05_08: 月次リターンの追加

stock_data <- stock_data %>% 
  group_by(firm_ID) %>% 
  mutate(lagged_stock_price = lag(stock_price)) %>% 
  ungroup() %>% 
  mutate(R = ((stock_price + DPS) * adjustment_coefficient - lagged_stock_price) / lagged_stock_price) # (5.1)式に従って月次リターンを計算
# ch05_09: 月次超過リターンの追加

stock_data <- stock_data %>%
  mutate(Re = R - R_F) # 月次超過リターンを計算

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.3.2 年次リターンの計算

# ch05_15: 月次リターンを累積して年次リターンを計算 (1)

annual_stock_data <- stock_data %>% 
  group_by(firm_ID, year) %>% # firm_IDとyearのペアでグループ化
  summarize(R = prod(1 + R) - 1, # バイ・アンド・ホールドの年次リターン
            R_F = prod(1 + R_F) - 1) %>%
  mutate(Re = R - R_F) %>%
  select(firm_ID, year, R, Re, R_F) %>%
  ungroup()
# ch05_16: 月次リターンを累積して年次リターンを計算 (2)

stock_data %>% 
  group_by(firm_ID, year) %>% 
  summarize(simple_cumulative_R = sum(R)) %>% # 元本が一定となるよう毎月リバランスした場合の年次リターン
  ungroup()

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.2 推定量と推定値の違い

# ch05_28: firm_IDが1の企業の平均月次超過リターンを計算

stock_data %>% 
  filter(firm_ID == 1) %>% 
  select(Re) %>%
  unlist() %>% # データフレームからベクトルに変換
  mean(na.rm = TRUE) # 第一引数は月次超過リターンのベクトル

## [1] 0.02906058

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.121296

5.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)

5.6.4 対数回帰モデル

# ch05_37: 線形・対数モデルによる推定

tidy(lm(Re ~ log(lagged_BEME), data = lm_sample_data)) # 右辺のみlog()関数で自然対数を取る
# ch05_38: データの保存

write_csv(monthly_data, "ch05_output1.csv")
write_csv(annual_data, "ch05_output2.csv")

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) # 無配の企業を落とす
# ch07_09: 潜在配当成長率の推定

dividend_yield_2020 %>% 
  left_join(FF3_cost_of_capital, by = "firm_ID") %>% # 株式資本コストの推定値を追加
  mutate(implied_dividend_growth = (cost_of_capital - dividend_yield) / (1 + dividend_yield)) # 配当割引モデルから潜在配当成長率を逆算

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.361662758

7.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.08371003

7.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.1 イベントの特定

# ch08_01: 外部パッケージとデータの読み込み

library(tidyverse)
library(broom)

return_data <- read_csv("ch08_return_data.csv")
market_return_data <- read_csv("ch08_market_return_data.csv")
event_data <- read_csv("ch08_event_data.csv")

head(return_data) # 各データの冒頭を確認
head(market_return_data)
head(event_data)

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を最終列に移動
# ch08_10: 累積異常リターン(CAR)の計算

event_window_data <- full_sample_data %>% 
  filter(relative_days >= -N2) %>% # イベント期間のデータを抽出
  mutate(R_normal = alpha + beta * R_M,
         AR = R - R_normal) %>%
  group_by(event_ID) %>% 
  mutate(CAR = cumsum(AR)) %>%
  ungroup() 

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)