5  株式データの取得と可視化

パッケージの読み込みとテーマ

pacman::p_load(
  tidyverse, # 便利なパッケージ群
  ggthemes, # ggplot2のテーマ集
  e1071, # 機械学習アルゴリズム
  broom, # モデルの整形
  scales # 軸のスケール調整
)

mystyle <- list(
  theme_economist_white(
    base_family = "HiraKakuProN-W3"
    ),
  scale_colour_economist(),
  theme(
    text = element_text(size = 10),
    axis.title = element_text(size = 10)
  )
)

前章では財務データを題材に、クロスセクションデータを扱うために必要なデータの取得と操作と可視化について学習しました。 特に重要な操作として，mutate()やsummarise()関数で、.by = varを指定することでグループごとの変数の値を計算する方法を学びました。 本章では株式データを題材に，株式データの理解に必要な株式の基礎知識とともに，リターンの定義や計算について学び，最後には統計的推論と線形回帰分析についても学習します。

5.0.1 株価データのダウンロードと読み込み

いままでと同じように、readrパッケージのread_csv()関数でCSVファイルを読み込みます。 ここでは、ch05_stock_data.csvを使います。

株式データの読み込み

stock_data <- read_csv("data/ch05_stock_data.csv")
print(stock_data)

# A tibble: 95,040 × 9
    year month month_ID firm_ID stock_price   DPS shares_outstanding
   <dbl> <dbl>    <dbl>   <dbl>       <dbl> <dbl>              <dbl>
 1  2015     1        1       1         954     0            2422000
 2  2015     2        2       1         960     0            2422000
 3  2015     3        3       1        1113     0            2422000
 4  2015     4        4       1        1081     0            2422000
 5  2015     5        5       1        1317     0            2422000
 6  2015     6        6       1        1366    29            2422000
 7  2015     7        7       1        1353     0            2422000
 8  2015     8        8       1        1209     0            2422000
 9  2015     9        9       1        1291     0            2422000
10  2015    10       10       1        1407     0            2422000
# ℹ 95,030 more rows
# ℹ 2 more variables: adjustment_coefficient <dbl>, R_F <dbl>

9変数，95,040行という大規模なデータが読み込まれました。 紙面の都合上、長い変数名を短く省略します。変数名を変えるには、dplyr::rename()を使います。

変数名の変更

stock_data <- stock_data |>
  rename(
    n_shares = shares_outstanding,
    adj_coef = adjustment_coefficient
  )

読み込まれたデータstock_dataの中には以下の変数が含まれています。

列名	データ型	説明
year	数値型dbl	年度
month	数値型dbl	月
month_ID	数値型dbl	月ID
firm_ID	数値型dbl	企業ID
stock_price	数値型dbl	月末時点での終値
DPS	数値型dbl	一株当たり配当額
n_shares	数値型int	月末時点での発行済株式数
adj_coef	数値型int	調整係数
R_F	数値型dbl	月次無リスク金利

year〜firm_IDまでは，時系列や企業の特定に必要なキーとなり， stock_priceからR_Fが分析する対象となる株価や配当，発行済株式数，調整係数，無リスク金利となります。 7列目のn_sharesは発行済株式数です。 発行済株式数は随時，新株発行や自社株買い，株式分割で変動するため，株式数の変化を調整するためのデータとしてadj_coefがあります。 旧株式1に対して新株式2となる株式分割が行われた場合，この調整係数は2となります。

このような大規模データを前にした場合，とりあえずデータの構造を把握することから始めます。 ここでは基本関数よりも多機能なtidyverseのdplyrパッケージに含まれるglimpse関数を使って，データの構造を確認してみます。

glimpse(stock_data)

Rows: 95,040
Columns: 9
$ year        <dbl> 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015…
$ month       <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7…
$ month_ID    <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,…
$ firm_ID     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ stock_price <dbl> 954, 960, 1113, 1081, 1317, 1366, 1353, 1209, 1291, 1407, …
$ DPS         <dbl> 0, 0, 0, 0, 0, 29, 0, 0, 0, 0, 0, 29, 0, 0, 0, 0, 0, 40, 0…
$ n_shares    <dbl> 2422000, 2422000, 2422000, 2422000, 2422000, 2422000, 2422…
$ adj_coef    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ R_F         <dbl> 6.506826e-04, 5.834099e-04, 6.114423e-04, 6.848180e-04, 7.…

データの型がdblという初めて出てきた型になってますが，これはdoubleの略で，数値型の一種です。doubleは64ビットの浮動小数点数を表し，小数点以下15桁までの精度を持ちます。

Tip数値型の種類

数値型numericは，doubleとintのどちらかになります。 intは整数(integer)型で，32ビットの整数を表し，小数点以下は切り捨てられます。 dblはintの上位互換で，小数点以下の桁数が多い場合に使われます。

たとえば，firm_IDが74，month_IDが29から32（つまり2017年5月から8月）のデータを抽出してみます。

データの抽出

stock_data |>
  # 条件を満たすデータをfilter()で抽出
  filter(firm_ID == 74 , month_ID %in% 29:32) |>
  # 特定の変数を抽出
  select(
    year, firm_ID, month_ID,
    stock_price, n_shares, adj_coef
  )

# A tibble: 4 × 6
   year firm_ID month_ID stock_price n_shares adj_coef
  <dbl>   <dbl>    <dbl>       <dbl>    <dbl>    <dbl>
1  2017      74       29        2816  4960000        1
2  2017      74       30        1402  9920000        2
3  2017      74       31        1420  9920000        1
4  2017      74       32        1502  9920000        1

出力された結果の2行のadj_coefが2となっており、同時にn_sharesが倍になっていることから、この会社は1株を2株に株式分割していることがわかります。

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

時価総額(market capitalization)を計算するためには、株式数に株価を掛けて計算します。 新しい変数を作成するときはdplyrパッケージのmutate()関数を使います。 mutate()で時価総額を表す新しい変数MEを作成します。

時価総額の計算

stock_data <- stock_data |>
  mutate(
    ME = stock_price * n_shares # 時価総額
    )

次に時価総額のヒストグラムを作成してみます。 前節で学習した内容に加えて、いろいろ見た目の指定を増やしてみます。 scale()関数を使って軸の範囲や表記方法を細かく指定します。

時価総額のヒストグラム

# 日本語を含むグラフを作成するときは，dev = "quartz_pdf"を指定
g <- ggplot(stock_data) + aes(x = ME) +
  geom_histogram() + # ヒストグラムを描く
  scale_x_continuous( # 軸の範囲と表記方法の指定
    limits = c(0, quantile(stock_data$ME, 0.95)),
    labels = label_comma(scale = 1e-6)
  ) + xlab("時価総額") + ylab("度数") +  mystyle
print(g) # グラフを出力

5.1.1 トータル・リターンと超過リターンの計算

ある株式のt期のトータル・リターンR_tは、次式で定義されます。

R_t = \frac{(\text{株価}_t + \text{1株当り配当}_t) \times \text{調整係数}_t - \text{株価}_{t-1}}{\text{株価}_{t-1}}

たいていの場合，1株当り配当\text{1株当り配当}_tはゼロで，調整係数は1となるため，次式のように簡略化できます。

R_t = \frac{\text{株価}_t - \text{株価}_{t-1}}{\text{株価}_{t-1}}

銘柄ごとにリターンを計算するので，mutate()と.by = firm_IDを使って計算します。 ここでは，firm_IDごとにトータルリターンRと，トータルリターンから無リスク利子率を除いた超過リターンReを計算しています。

stock_data <- stock_data |>
  mutate( # 新しい変数を作成
    R = ((stock_price + DPS) * adj_coef - lag(stock_price)) / lag(stock_price),
    Re = R - R_F, # 月次超過リターン
    .by = firm_ID # firm_IDごとにグループ
  )

ここまでの処理でstock_dataに時価総額ME，トータルリターンR，超過リターンReが追加されました。 以下では，このデータを使って探索的データ分析を行います。

5.1.2 株式データの探索的データ分析

探索的データ分析という名の，とりあえず何も考えずに目の前のデータから何が分かるのかを調べ倒してみる，という分析を行います。

とりあえず，summary()で基本統計量を確認します。

summary(stock_data)

      year          month          month_ID        firm_ID      
 Min.   :2015   Min.   : 1.00   Min.   : 1.00   Min.   :   1.0  
 1st Qu.:2016   1st Qu.: 3.75   1st Qu.:19.00   1st Qu.: 384.0  
 Median :2018   Median : 6.50   Median :37.00   Median : 760.0  
 Mean   :2018   Mean   : 6.50   Mean   :37.01   Mean   : 761.2  
 3rd Qu.:2019   3rd Qu.: 9.25   3rd Qu.:55.00   3rd Qu.:1147.0  
 Max.   :2020   Max.   :12.00   Max.   :72.00   Max.   :1515.0  
                                                                
  stock_price          DPS              n_shares            adj_coef     
 Min.   :   112   Min.   :   0.000   Min.   :3.700e+04   Min.   :0.1000  
 1st Qu.:  1417   1st Qu.:   0.000   1st Qu.:3.151e+06   1st Qu.:1.0000  
 Median :  2445   Median :   0.000   Median :1.014e+07   Median :1.0000  
 Mean   :  4685   Mean   :   6.802   Mean   :6.973e+07   Mean   :0.9999  
 3rd Qu.:  4558   3rd Qu.:   0.000   3rd Qu.:3.564e+07   3rd Qu.:1.0000  
 Max.   :622796   Max.   :1913.000   Max.   :2.111e+10   Max.   :2.0000  
                                                                         
      R_F                   ME                  R                  Re          
 Min.   :-2.329e-04   Min.   :1.459e+08   Min.   :-0.38028   Min.   :-0.38020  
 1st Qu.: 8.233e-06   1st Qu.:9.648e+09   1st Qu.:-0.04057   1st Qu.:-0.04076  
 Median : 8.203e-05   Median :2.563e+10   Median : 0.01033   Median : 0.01013  
 Mean   : 2.041e-04   Mean   :1.365e+11   Mean   : 0.01586   Mean   : 0.01565  
 3rd Qu.: 4.626e-04   3rd Qu.:7.987e+10   3rd Qu.: 0.06481   3rd Qu.: 0.06460  
 Max.   : 7.368e-04   Max.   :3.293e+13   Max.   : 0.51498   Max.   : 0.51448  
                                          NAs    :1515       NAs    :1515      

さらに，分散と標準偏差も計算します。データには欠損値が含まれているため，na.rm = TRUEオプションをつけています。 教科書とは違いますが，dplyr::summarize()関数を使って一気に複数の統計量を計算します。 ここでは，stock_data <- stock_dataとしていないため，stock_dataには新しい変数は追加されず，結果を表示するだけです。

stock_data |>
  summarise(
    var_R = var(R, na.rm = TRUE), # 総リターンの分散
    sd_R  = sd(R, na.rm = TRUE), # 総リターンの標準偏差
    var_Re = var(Re, na.rm = TRUE), # 超過リターンの分散
    sd_Re = sd(Re, na.rm = TRUE) # 超過リターンの標準偏差
  )

# A tibble: 1 × 4
    var_R   sd_R  var_Re  sd_Re
    <dbl>  <dbl>   <dbl>  <dbl>
1 0.00830 0.0911 0.00830 0.0911

データの分布の形を表す統計量である，3次のモーメントである歪度(skewness)と4次のモーメントである尖度(kurtosis)も計算してみます。 たとえば確率変数xの歪度は \text{歪度} = \frac 1n \sum _{i = 1}^n \left( \frac{x_i - \bar x}{\hat \sigma_x} \right) ^3 で表され、正の値をとるとき分布が右に歪んでいる(裾が右に長い)ことを示している。 尖度は， \text{尖度} = \frac 1n \sum _{i = 1}^n \left( \frac{x_i - \bar x}{\hat \sigma_x} \right) ^4 で定義され，尖度が3のとき，正規分布と同じ尖度を持つことを示し，3より大きいとき、正規分布よりとがった分布をしていることを示します。

歪度と尖度を計算するには，e1071パッケージのskewness()関数を使います。

stock_data |>
  summarise(
    skewness_R = skewness(R, na.rm = TRUE), # 総リターンの歪度
    kurtosis_R = kurtosis(R, na.rm = TRUE), # 総リターンの尖度
    skewness_Re = skewness(Re, na.rm = TRUE), # 超過リターンの歪度
    kurtosis_Re = kurtosis(Re, na.rm = TRUE) # 超過リターンの尖度
  )

# A tibble: 1 × 4
  skewness_R kurtosis_R skewness_Re kurtosis_Re
       <dbl>      <dbl>       <dbl>       <dbl>
1      0.507       1.27       0.507        1.27

トータルリターンの歪度が0.507 >0なので，右に裾の広い分布となっており，尖度が1.27 < 3なので正規分布よりもとがった形となっていることがわかります。

図でも確認するために，トータルリターンRのヒストグラムを書いてみます。

ggplot(stock_data) +
  aes(x = R) + geom_histogram() + # トータルリターンのヒストグラム
  xlab("月次株式リターン") + ylab("度数") + mystyle# 軸ラベル

トータルリターンのヒストグラムに正規分布のグラフを重ねてみると，次のようになります。

# トータルリターン
R <- stock_data$R
# 平均と標準偏差を計算
m <- mean(stock_data$R, na.rm = TRUE)
s <- sd(stock_data$R, na.rm = TRUE)

ggplot(stock_data) +
  aes(x = R) +
  # y軸を密度(density)に設定
  geom_histogram(aes(y = after_stat(density)), bins = 100, alpha = 0.3, fill = "gray50") +
  geom_density(color = "blue", linewidth = 1) +  # 密度曲線を追加
  # 正規分布の曲線を追加
  stat_function(
    fun = dnorm, args = list(mean = m, sd = s),
    color = "red", linetype = "dashed", linewidth = 1) +
  xlab("月次株式リターン") + ylab("密度") + mystyle

歪度と尖度の結果と整合的に、実際のトータル・リターンのデータの確率密度(青い線)と正規分布(赤い点線)を比べてみると、実際のデータは右に裾が長く，尖度も正規分布よりもとがった形となっていることが分かります。

より厳密にデータが正規分布にしたがっているかを確認するために，Q-Qプロットを描いてみます。

Q-Qプロット

ggplot(stock_data) +
  aes(sample = R) +
  stat_qq() +
  stat_qq_line(color = "red") +
  labs(x = "理論上の分位点", y = "サンプルの分位点") +
  mystyle

5.2 リターンの累積

5.2.1 バイ・アンド・ホールド・リターンの考え方

バイ・アンド・ホールド・リターン(buy-and-hold-return)は，ある時点で株式を購入し，そのまま保有し続けたときのリターンのことを指します。 たとえば，12月末に株式を購入し，3月末に売却したときの、バイ・アンド・ホールド・リターンの累積は，次式で計算できます。ただし配当は無視します。 \begin{aligned} \left ( \frac{W_{\text{3月}}}{W_{\text{12月}}} \right) = \underbrace{\left ( \frac{W_{\text{1月}}}{W_{\text{12月}}} \right)}_{1 + R_{\text{1月}}} \times \underbrace{\left ( \frac{W_{\text{2月}}}{W_{\text{1月}}} \right)}_{1 + R_{\text{2月}}} \times \underbrace{\left ( \frac{W_{\text{3月}}}{W_{\text{2月}}} \right)}_{1 + R_{\text{3月}}} \times \end{aligned}

一般的に、月次でt時点からT時点までの年次リターンは、

\left ( \frac{W_{\text{翌年12月末}}}{W_{\text{12月末}}} \right) = \prod_{t = \text{1月}}^{\text{12月}} (1 + R_t)

と計算できます。

5.3 年次リターンの計算

では、stock_dataをもとに年次リターンを計算してみましょう。

annual_stock_data <- stock_data |>
  summarise(
    annual_R = prod(1 + R) - 1, # B&H年次リターン
    annual_R_F = prod(1 + R_F) - 1, # 年次超過リターン
    .by = c(firm_ID, year) # 企業と年度でグループ化
  ) |>
  mutate(
    annual_Re = annual_R - annual_R_F # 年次超過リターン
    )
head(annual_stock_data)

# A tibble: 6 × 5
  firm_ID  year annual_R annual_R_F annual_Re
    <dbl> <dbl>    <dbl>      <dbl>     <dbl>
1       1  2015   NA      0.00743      NA    
2       1  2016    0.997  0.000565      0.997
3       1  2017    0.688  0.0000488     0.688
4       1  2018   -0.214  0.00579      -0.219
5       1  2019    0.647 -0.000770      0.648
6       1  2020   -0.284  0.000380     -0.285

5.4 株式データと財務データを組み合わせた分析

ここでは，株式データと財務データを組み合わせて分析を行います。

5.4.1 データの結合

複数のデータフレームを結合する際に重要なところは、

1. データの頻度の違い
2. タイミングの一致

となる。 今まで使ってきた財務データは年次データであるのに対して，株式データは月次データであるため，データの頻度が異なります。 確認してみましょう。

financial_data <- readr::read_csv("data/ch04_output.csv")
nrow(financial_data) # 年次財務データの行数

[1] 7919

nrow(annual_stock_data) # 年次リターン・データの行数

[1] 7920

nrow(stock_data) # 月次リターン・データの行数

[1] 95040

月次リターンの行数が上の年次データとは大きく異なっていることが分かります。

先読みバイアス(look-ahead bias)とは、ある時点での情報を使って、その時点よりも未来の情報を使っていることを指します。12月末決算の会社のディスクロジャージャーは、最速で決算日後45日以内に出される決算短信か、3ヶ月以内に出される有価証券報告書があります。 このため、年次データを使って同時期の年次リターンを計算すると、年次データの発表後に出される有価証券報告書の情報を使っていることになります。これを先読みバイアスといいます。

annual_data <- annual_stock_data |>
  full_join(
    financial_data,
    by = c("firm_ID", "year")
  ) # firm_IDとyearのペアをキーとして設定

annual_data <- annual_stock_data |>
  full_join(financial_data) # キーを省略した場合，列名が同じ変数がキーになる

monthly_data <- stock_data |>
  full_join(financial_data, by = c("firm_ID", "year")) # firm_IDとyearのペアをキーとして設定

5.5 バブルチャート

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)

year2015 <- annual_data |>
   filter(
      year == 2015, # 2015年のみ
      firm_ID %in% 2:20, # firm_IDが2から20のデータを抽出
      X > 0 # 対数を取るため当期純利益が正のデータのみ抽出
      )

ggplot(year2015) +
  aes(x = log(sales), y = log(X), size = ME, alpha = 0.4) +
  geom_point() + # バブルチャートを描くには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)) + mystyle

5.6 統計的推論入門

5.6.1 リターン・データに関する仮定

統計的推論の内容に入る前に、firm_IDが1の企業の超過リターンReのデータを眺めてみます。

stock_data_1_month <- stock_data |>
  filter(firm_ID == 1) |>
  select(month_ID, Re, R)
ggplot(stock_data_1_month) +
  aes(x = month_ID, y = Re) + # 軸の設定
  geom_line() + mystyle # 折れ線グラフ

このようなデータは時系列データと呼ばれ、ある個体(ここではfirm_IDが1の企業)の一定期間にわたって観測したデータのことを指します。

Important仮定

月次超過リターンの時系列データは、何らかの確率分布(モデル)から独立に生成されている。

1. 株価そのものでは無く変化率であるリターンをモデル化する理由

株価それ自体は株式分割などの要因でも変化するし、会社の成長に応じて上昇するため、その成長率をモデル化する方が、投資のリスクに対するリターンをより明確に評価することができるからである。

2. 月次リターンではなく月次超過リターンをモデルかする理由

投資リスクを取って得られる追加的なリターンであるリスクプレミアムを明確に評価するために、無リスク金利を差し引いた超過リターンをモデル化します。

3. 月次超過リターンが独立であるとする理由

系列相関（または自己相関）とは、同じ時系列データにおいて、時点tのリターンと時点t-kのリターン（k期前のリターン）の間に相関がある状態のことです。例えば、前月のリターンが高ければ今月のリターンも高くなる傾向、というような関係が系列相関です。

効率的市場仮説（Efficient Market Hypothesis; EMH）に基づく半強度の効率的市場であれば、過去に公開された全ての情報はすでに株価に織り込まれているとされています。したがって、過去の超過リターンは将来の超過リターンと無関係であり、超過リターンに系列相関（自己相関）がないと仮定できます。

しかし、現実には以下のような理由で系列相関が観測されることが知られています：

• モメンタム効果：過去のリターンが高かった銘柄は、その後もリターンが高い傾向
• 反転現象（mean reversion）：極端に高いリターンの後には低いリターンが続く傾向
• アノマリー：市場が有効でない期間における異常な価格パターン

これらのため、現実には系列相関があると考えられます。しかし、系列相関を考慮するとモデルが著しく複雑になるため、ここでは最初の学習段階として、独立性の仮定を置いてモデルを単純化します。より高度な分析では、ARMA（自己回帰移動平均）モデルなどで系列相関を明示的に扱う必要があります。

4. 月次超過リターンが正規分布に従うとする理由

正規分布を仮定するといろいろ計算が楽になる、という理由とともに、 株価の動きは多くの独立した事象の結果であると考えられ、中心極限定理により、超過リターンは正規分布に従うと考えられます。

これらの仮定を受け入れると、超過リターンの確率分布は次式で表されます。

Re_t \sim N(\mu, \sigma^2)

これにより、母集団分布に対して統計的推論が可能になります。 firm 1の超過リターンのヒストグラムを書いてみます。

ggplot(stock_data_1_month) + aes(x = Re) + # データと変数
  geom_histogram()

サンプルサイズが小さいため凸凹しているけれど、もっとサンプルを増やせば、正規分布に近い形をとるはずです。

5.7 推定量と推定値の違い

推定量(estimator)とは、母集団分布の母数(パラメータ)を推定するために使われる統計量のことです。 推定値(estimates)とは、推定量に実際のデータを代入して計算した値のことです。 たとえば母集団分布が正規分布N(\mu, \sigma^2)に従うと仮定すると、母数は\muと\sigma^2の2つです。この母数を推定するために使われる統計量が、それぞれ標本平均\bar xと標本分散s^2です。 このときの推定値とは、実際に観察された標本から計算された標本平均\bar xと標本分散s^2のことを指します。

次の問題を考えます。

Important問題

firm_IDが1の銘柄の月次リターンR_{i,t}^eは、期待値の意味で、ゼロより大きいのだろうか？

ここで「期待値の意味で」というのは平均的に、と言い換えても問題ないです。 企業1の月次超過リターンの平均は0より大きい、ということは、企業1の月次リターンは平均的に無リスク利子率よりも大きいかどうか、を比べるということです。

この問題を解くために、まずは母集団分布の母数である期待値\muを推定する必要があります。 期待値\muを推定するために使われる推定量は標本平均\bar xなので、ここで標本平均を計算します。

stock_data_1_month |>
  summarise(
    mean_Re = mean(Re, na.rm = TRUE)
  )

# A tibble: 1 × 1
  mean_Re
    <dbl>
1  0.0291

企業1の平均月次超過リターンmean_Reが0.0291となりました。 これだけみるとゼロより大きい値ですが、これは母集団から抽出された1つのサンプルの平均ですので、他のサンプルの平均がゼロを超えるかどうかは分かりません。

5.7.1 大数の法則

母集団から無限個の標本(sample)を抽出して、それぞれの標本平均(sample mean)を計算すると、その標本平均の平均は母集団の期待値\muに一致することが知られています。これを大数の法則(law of large number)といいます。

数式よりも前にシミュレーションで確認してみましょう。 まずは、母集団分布を平均が10、標準偏差が2の正規分布N(10, 4)として、母集団から100のデータを抽出して標本を1つ作ります。そしてその標本の平均を計算します。

set.seed(123) # 乱数の種を固定
size <- 100 # 標本サイズ
# 母集団のパラメータの設定
mu <- 10 # 平均
sigma <- 4 # 標準偏差
population <- rnorm(size, mu, sigma) # 母集団からサンプルを抽出
mean(population) # 標本平均

[1] 10.36162

平均は10.3616236となりました。 母集団の平均10とは異なる数値になっています。 もう一度別の標本でやってみると、

population <- rnorm(size, mu, sigma) # 母集団からサンプルを抽出
mean(population) # 標本平均

[1] 9.569813

また違う平均が計算されました。 この作業を何度も何度も繰り返し、標本平均をたくさん計算します。 ここでは、100個のデータをもつ標本を100個作って、それぞれの標本平均を計算します。

n_sample <- 100 # サンプル数
sample_mean <- rep(NA, n_sample) # 標本平均を格納するベクトル
sample_sd <- rep(NA, n_sample) # 標本分散を格納するベクトル
for(i in 1:n_sample){
  population <- rnorm(size, mu, sigma) # 母集団からサンプルを抽出
  sample_mean[i] <- mean(population) # 標本平均を計算
  sample_sd[i] <- sd(population)
}

サンプルの平均が100個計算できたので、ヒストグラムを書いてみます。

ggplot() + aes(x = sample_mean) + geom_histogram() + mystyle

あまり正規分布のようには見えません。 ではサンプルの平均の平均を計算してみます。

mean(sample_mean) # 標本平均の平均

[1] 9.99003

母集団の平均10に近い値になっていることが分かります。もっとサンプルの数を増やしてみます。 データの数が100のサンプルを1,000,000個(100万個)抽出して、それぞれのサンプルの平均値を計算します。

n_sample <- 10^6
# 標本平均のシミュレーション
sample_mean <- replicate(n_sample, mean(rnorm(size, mu, sigma)))
# 標本標準偏差のシミュレーション（必要な場合）
sample_sd <- replicate(n_sample, sd(rnorm(size, mu, sigma)))
# 作図
ggplot() + aes(x = sample_mean) + geom_histogram() + mystyle

ほぼ正規分布のような形をしていることが分かります。ではサンプルの平均の平均を計算してみます。

mean(sample_mean) # 標本平均の平均

[1] 10.00063

mean(sample_sd) # 標本分散の平均

[1] 3.989758

このように標本の数を無限に大きくしたとき、サンプルの平均の平均は母集団の平均10に一致するし，標本標準偏差は母集団の標準偏差4に一致する，というのが大数の法則です。

5.8 t値の計算

先ほど計算したfirm_IDが1の企業の超過リターンの平均値はNAでした。 これがゼロより大きいかどうかはすぐ分かりますが，この値はたまたま今手元にある1つのサンプルから計算された平均値なので，他の標本ではどうなるか分かりません。 このように推定量にばらつきがある場合には，その推定量の分布を考える必要があります。 ここでは，その分布をt分布(t distribution)と仮定して，t値(t-value)を計算してみます。 t値は次のように定義されます。

t = \frac{\bar{X} - \mu_0}{\sqrt{s^2 / n} } \stackrel{d}{\approx} N(0,1)

ここで，\bar Xは標本平均，\mu_0は帰無仮説(null hypothesis)の値で，ここでは\mu_0 = 0とします。s^2は標本分散，nは標本サイズです。 分子に注目すると，標本平均と帰無仮説の値の差となっており，もし標本平均が0に近いなら，t値は0に近い値になります。

Re_firm_ID_1 <- stock_data |>
  filter(firm_ID == 1) |> # firm_IDが1の企業のデータを抽出
  select(Re) |> #超過リターンのみ選択
  drop_na() |> # 欠損値を除去
  unlist() # データフレームをベクトルに変換

mu0 <- 0 # 帰無仮説の値
n <- length(Re_firm_ID_1) # 標本サイズ

t_value <- (mean(Re_firm_ID_1) - mu0) / sqrt(var(Re_firm_ID_1) / n) # $t$値の計算
print(t_value)

[1] 2.121296

5.8.1 統計的検定の考え方

あなたがサイコロを投げるゲームをしていて、あるプレイヤーが非常に幸運だと主張しています。彼は6回サイコロを投げて、5回も「6」が出たとします。これはただの偶然なのか、それとも何か他の要因（例えば、サイコロがいかさまであるとか）が関与しているのでしょうか？

この問いに答えるために、我々は統計的検定(statistical test)を用いることができます。 まず帰無仮説(null hypothesis)を設定します。 この例では、帰無仮説は「サイコロは公正であり、すべての出目が等確率（1/6）で出る」とすることが適切です。

次に、この帰無仮説が真(true)である場合に、我々が観察した結果（5回の「6」）がどれほどあり得ないかを計算します。これがp値(p value)です。この場合、6回投げて5回「6」が出る確率を計算します。

これを計算すると、p値は非常に小さいことが分かり（つまり、この結果は帰無仮説の下ではほぼあり得ない），帰無仮説が棄却されます。 帰無仮説が棄却されるとは，帰無仮説が正しいと仮定したときに，観測されたデータが得られる確率が小さいことを意味します。 したがって、我々はこの結果が偶然生じたとは考えにくく、その代わりにサイコロがいかさまである、または何か他の非ランダムな要因が作用している可能性を強く疑うことになります。これがp値を用いて統計的検定の考え方です。

p値が0.05より小さい場合，帰無仮説は棄却され，対立仮説が採択される，というケースが多いです。 この場合，有意水準5%で帰無仮説は棄却されます。 有意水準は，帰無仮説が正しいと仮定したときに，観測されたデータが得られる確率が小さいと判断する基準値です。 有意水準は，5%や1%がよく使われます。

Rでは，t.test()関数を使って，t値とp値を計算することができます。 ここでは，t.test()関数を使って，firm_IDが1の企業の超過リターンがゼロなのかどうなのか，を検定するために，t値とp値を計算してみましょう。

t.test(Re_firm_ID_1)

    One Sample t-test

data:  Re_firm_ID_1
t = 2.1213, df = 70, p-value = 0.03744
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.001737894 0.056383256
sample estimates:
 mean of x 
0.02906058 

5.9 線形回帰入門

年次財務データannual_dataを使って線形回帰(linear regression)について学習します。 線形(linear)とは，Xと\betaの積が線形であることを意味します。回帰(regression)とは，Xと\betaの積がYを説明することを意味します。 回帰は，2変数間の間で、一方の変数が他方の変数に対して影響を与えるという関係を想定できる場合に用いられます。 しかし因果関係を直接調べているのではないことに注意しましょう。

影響を与える変数を説明変数(explanatory variable)や独立変数(independent variable)といい、 影響を与えられる変数を目的変数(response variable)や従属変数(dependent variable)といいます。 分野などでどっちを使うのかは異なりますが、ここでは説明変数と目的変数を使います。

線形関係を仮定した関係を式にすると，次のようになります。 Y = \alpha + \beta X ここで， \alphaは定数項(constant term)と呼ばれ，切片(intercept)とも呼ばれます。 \betaは回帰係数(regression coefficient)と呼ばれ，傾き(slope)とも呼ばれます。

実際のデータが上のモデルを満たす、つまりすべてのデータが直線上に並んでいるわけではないのです。 実際のデータとモデルとの間にはずれがあることを考慮して、上のモデルを次のように書き換えます。

Y = \alpha + \beta X + u ここでuは誤差項(error term)とか観察不可能項(unobservable term)と呼ばれます。

線形回帰の目的は，XとYの関係を表す\alphaと\betaを推定することです。 推定するために，観測できるデータを使います。 観測できるデータはXとYのペアで(X_i, Y_i)と表記します。 このペアを観測値(observed value)と呼びます。 この観測値は，XとYの関係を表す\alphaと\betaを推定するために使われます。 推定するために使われる\alphaと\betaを推定値(estimated value)と呼びます。 推定値は，\hat{\alpha}と\hat{\beta}と表記して，観察値と区別します。

5.10 OLS推定

推定方法として最も有名なものが最小二乗法(ordinary least squares; OLS)です。 最小二乗法では，観測値と推定値の差の二乗の和が最小になるように\alphaと\betaを推定します。 このとき，観測値と推定値の差を残差(residual)と呼びます。 最小二乗法では，残差の二乗の和である残差平方和(sum of squared residuals; SSR)が最小になるように\alphaと\betaを推定します。

\min _{\alpha, \beta} \sum _{i=1}^n (Y_i - \alpha - \beta X_i)^2

# 線形回帰分析のための事前準備
# ch05_24ですでにME（時価総額）が結合されていることが前提

annual_data <- annual_data |>
  arrange(firm_ID, year) |> # 時系列順に並べ替え（lag計算に必須）
  mutate(
    lag_ME = lag(ME),
    lagged_BEME = lagged_BE / lag_ME, # ここでBEMEを計算
    .by = firm_ID
  )

lm_sample_data <- annual_data |>
  filter(year == 2016, firm_ID <= 10) |>
  # lagged_BEMEは計算済みなので、selectで選ぶだけにする
  select(firm_ID, year, Re = annual_Re, lagged_BEME) |>
  drop_na() # 欠損値を除去

ggplot(lm_sample_data) +
  geom_point(aes(x = lagged_BEME, y = Re)) + # 散布図
  xlab("簿価時価比率") + ylab("超過リターン") + mystyle

# 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()

# ch05_33: lm()関数を用いた線形回帰

lm_results <- lm(Re ~ lagged_BEME, data = lm_sample_data) # 従属変数 ~ 独立変数
names(lm_results)

 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model"        

print(lm_results$coefficients)

(Intercept) lagged_BEME 
 0.05513385  0.07552921 

# ch05_35: broomパッケージのtidy()関数で係数の推定値に関する結果を確認
tidy(lm_results)

# A tibble: 2 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)   0.0551    0.156      0.354   0.736
2 lagged_BEME   0.0755    0.0844     0.895   0.405

glance(lm_results)

# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>
1     0.118       -0.0294 0.223     0.800   0.405     1   1.81  2.37  2.61
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

5.11 対数回帰モデル

独立変数Xが変化したときの従属変数Yへの影響は一定(つまり傾きが一定)と仮定してきましたが、 実際には傾きが一定でない場合もあります。 回帰式が非線形であることが想定される場合、対処法として

• 多項式回帰(polynomial regression) ：独立変数にX^2とかX^3を加える
• 対数回帰(logarithmic regression) ：独立変数や従属変数の対数をとる

たとえば、独立変数を対数変換した場合は、次のようなモデルになる。

Y_i = \beta_0 + \beta_1 \log (X_i) + \varepsilon_i

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

tidy(lm(Re ~ log(lagged_BEME), data = lm_sample_data)) # 右辺のみlog()関数で自然対数を取る

# A tibble: 2 × 5
  term             estimate std.error statistic p.value
  <chr>               <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)        0.167     0.0868     1.93    0.102
2 log(lagged_BEME)   0.0355    0.109      0.326   0.756

作成したデータフレームをcsvファイルとして保存するには，write_csv()関数を用います。 前処理が終わった後や新しい変数を作った後に、データを保存しておくと便利です。 6章以降では、以下のデータを継続して使うので、csvファイルとして保存しておきます。

データの保存

dir.create("data", showWarnings = FALSE, recursive = TRUE)

safe_write_csv <- function(df, path) {
  tryCatch(
    write_csv(df, path),
    error = function(e) {
      warning(paste0("CSVの保存に失敗しました: ", path, " / ", conditionMessage(e)))
      invisible(NULL)
    }
  )
}
# safe_write_csv()は、
safe_write_csv(monthly_data, "data/ch05_output1.csv")
safe_write_csv(annual_data, "data/ch05_output2.csv")