Yamamotoの日記

金融工学、スポーツ、その他当たり障りのないことを書きます。個別の投資術や銘柄を推奨する意図はなく、所属組織とは関係ありません。

ダルビッシュ・大谷・山本・佐々木を分析する

本記事は個人の見解に基づくものであり、特定の金融商品やサービスを推奨するものではありません。また、所属組織には一切関係ありません。

10/26 時点・ドジャースが 2 年連続でワールドシリーズ進出

そして昨日は山本投手のポストシーズン完投勝利 2 試合目!Santa Monica のブリュワリーで試合を見ていて、めちゃくちゃ盛り上がりました。

筆者は LA に住んでまして、地元のチームはドジャースです。大谷・山本・佐々木がいることもあり、日本人に大人気。また、ダルビッシュは自分が高校生のときから応援していた選手で、サンディエゴ・パドレスにいてくれると気軽に見に行けてありがたいです。チームはライバルチームなので応援しているとややこしいですが。。。先日、サンディエゴに行った際にはダルビッシュの先発試合で、とても楽しかったです。

ちなみに、ドジャースのチケットはパドレスの倍くらいの値段がするのでなかなか財布に厳しいです。

僕は、この中では特にダルビッシュと山本を応援しています。二人とも、ロジカルに考えて工学的なアプローチをした上でしっかり結果を出しています。また、その結果は自己責任なんだという精神性と覚悟がとてもかっこいいと思います。新しいことをどんどん試す一方で、他領域や昔の知見もしっかり取り入れていて、自分の仕事や日常生活でも見倣って行きたいと思える選手たちです。

例えばこのインタビューではどんな覚悟でやっている選手なのかわかります。これは渡米前のインタビューですが、非常にアメリカ向きの性格だなあと感じました(日本でよくやってこれたな笑)。山本選手のインタビュー - 自分が納得していないの聞いて、もしダメでも自分がクビになるじゃないですか - 「責任取って辞める」とかいう人いますけど、その人が辞めたからといって、関係ないじゃないですか。僕の野球人生にとって。

Qiita でこんな記事を見つけました

山本由伸投手の投球分析

アメリカの MLB では全投球情報が公表されており、個人でも簡単にデータ分析をすることができます。さすがデータアナリティクス・セイバーメトリクスの国!この記事を参考に、同じことを日本人投手 4 名の 2020 年~2025 年の投球に拡張して実施してみました。

データの期間は 2020/1/1 - 2025/10/14 で、ポストシーズンも一部入った結果になっています。ドジャースの 3 人は直近で大きなフォーム変更可能性があるイベントがあったので、もっと詳しく掘り下げる際は、分析期間を適切に切った方が良さそうです。

使用データのアップデートスクリプト

Pitcher 情報 update 用スクリプト (アップデートに数時間かかります。また、先方サーバーに負荷をかけている可能性があります。ご使用は自己責任でお願いします。)

今回の記事では、可視化のみやります。

リリース位置

投手ごとにリリース位置を見てみると、ダルビッシュと大谷は結構ブレが大きく見えます。他の二人に比べて長いシーズンを投げているので、変化球ごとの違いだけでなく、経年変化も影響しているのではないかと思います。僕自身はサッカーしかやらないので野球の身体メカニズムについては詳しくないのですが、前述の記事で言及されているように、怪我の前後で(フォーム・)リリースポイントが変わることがあるようで、ダルビッシュも大谷もこの期間に怪我を経験しているので、おそらくそれも影響しているでしょう。

また、ダルビッシュの球種の多さにびびりますね。ダルビッシュのせいで(おかげで)凡例が大きくなってしまう。。。

ここでは、山本と佐々木のリリースポイントの近さに注目したいと思います。前述の記事の通り、山本は怪我の前後でリリースポイントが左下に移動しているので見かけ上は少し広く見えますが、怪我前後それぞれに集中してみると、他 2 投手と比べてリリース位置がぶれないことことがわかります。

pitch_jp |>
    ggplot() +
    aes(x = release_pos_x, y = release_pos_z, color = pitch_name) +
    geom_point(alpha = 0.3) +
    ylim(4.5, 7.0) +
    xlim(-4.0, 0.0) +
    labs(
        title = "Pitch Release Position",
        x = "Release Position X (ft)",
        y = "Release Position Z (ft)",
        color = "Pitch Name"
    ) +
    facet_wrap(~player_name)

球種ごとに平均し、1 標準偏差でエラーバーを描いてみました。大谷と山本はそれぞれ 2023 年、2024 年の怪我前後でグループわけしています(かなりわかりにくいですが、どちらも怪我を通じて左下に移動しているように見えます)。こうしてみると、ダルビッシュは球種ごとの平均差がほとんどなく、球種内の標準偏差がかなり大きいので、リリースポイントで球種を絞るのはほぼ不可能ですね。

ダルビッシュと大谷の結果はシーズンの長さに起因するものかな?と思い、2023 年以降に絞って同じことをしてみましたが、あまり結果は変わりませんでした。

佐々木 or 山本とダルビッシュを比較すると、バッターからすればいずれも球種が絞りにくいものの、その理由は微妙に違いますね。ダルビッシュは「別の球種に見える」、その他二人は「どっちかわからない」という感じでしょうか。このあたり、選手によってアプローチが違ったりするのでしょうか?詳しい人に教えていただきたいところです。

pitch_jp |>
    # filter(game_date > ymd(20230101)) |>
    mutate(
        flag_after_injury =
            (str_detect(player_name, "Yamamoto") & (game_date > ymd(20240617))) |
                (str_detect(player_name, "Ohtani") & (game_date > ymd(20230831)))
    ) |>
    group_by(player_name, pitch_name, flag_after_injury) |>
    summarize(
        ave_pos_x = mean(release_pos_x),
        sd_pos_x = sd(release_pos_x),
        ave_pos_z = mean(release_pos_z),
        sd_pos_z = sd(release_pos_z),
        .groups = "drop"
    ) |>
    ggplot() +
    aes(x = ave_pos_x, y = ave_pos_z, color = pitch_name, shape = flag_after_injury) +
    geom_point(alpha = 0.3) +
    geom_errorbar(aes(ymin = ave_pos_z - sd_pos_z, ymax = ave_pos_z + sd_pos_z), width = 0.05) +
    geom_errorbarh(aes(xmin = ave_pos_x - sd_pos_x, xmax = ave_pos_x + sd_pos_x), height = 0.05) +
    ylim(4.5, 7.0) +
    xlim(-4.0, 0.0) +
    labs(
        title = "Average Pitch Release Position by Pitch type",
        x = "Release Position X (ft)",
        y = "Release Position Z (ft)",
        color = "Pitch Name",
        shape = "After Injury?"
    ) +
    facet_wrap(~player_name)

球種ごとの結果

B はボール、S はストライク、X はバットにボールが当たったもので内訳を計算してみました。もとの記事では、山本は特にカーブが取り上げられていましたが、ストレートもかなりストライクを取っていますね。

pitch_jp |>
    group_by(player_name, pitch_name, type) |>
    summarize(num = n()) |>
    group_by(player_name, pitch_name) |>
    mutate(ratio = num / sum(num)) |>
    ggplot() +
    aes(x = pitch_name, y = ratio, fill = type) +
    facet_wrap(~player_name) +
    geom_bar(stat = "identity", position = "fill") +
    labs(
        title = "Ball, Strike, and Hit Ratio by Pitch Type",
        x = "Pitch Name",
        y = "Ratio",
        fill = "Type"
    ) +
    coord_flip()

投球コース

どの選手も投球が多すぎておもしろい分析はあまりできなさそうですが、佐々木と山本のフォークは常に低めに決まっていて素晴らしいですね。山本は、左下に投げるときはフォーク、右下に投げるときはカーブが多いようです。色の感じから、高めにはストレートが多いのがよくわかります。

pitch_jp |>
    ggplot() +
    aes(x = plate_x, y = plate_z, color = pitch_name) +
    geom_point(alpha = 0.3) +
    labs(
        title = "Pitch Course at Home Plate",
        x = "Plate Position X (ft)",
        y = "Plate Position Z (ft)",
        color = "Pitch Name"
    ) +
    annotate(
        "rect",
        xmin = -0.83, xmax = 0.83,
        ymin = 1.5, ymax = 3.5,
        alpha = 0.3, fill = "black"
    ) +
    xlim(-3, 6) +
    ylim(-2, 6) +
    facet_wrap(~player_name)

リリースポイントと同じように、平均してみました。記事がかなり長くなってきたので、コードは割愛します。こうしてみると、どの投手も低めに集めており、さすがです。

打たれたボール

これも特に価値のあることは言えなさそうですが、一応載せておきます。

pitch_jp |>
    filter(type == "X") |>
    ggplot() +
    aes(x = plate_x, y = plate_z, color = pitch_name) +
    geom_point(alpha = 0.3) +
    labs(
        title = "Hit Pitch Course at Home Plate",
        x = "Plate Position X (ft)",
        y = "Plate Position Z (ft)",
        color = "Pitch Name"
    ) +
    annotate(
        "rect",
        xmin = -0.83, xmax = 0.83,
        ymin = 1.5, ymax = 3.5,
        alpha = 0.3, fill = "black"
    ) +
    xlim(-3, 6) +
    ylim(-2, 6) +
    facet_wrap(~player_name)

今後の課題

今後は、セイバーメトリクスの基本をなぞっていくような分析をしてみたいですね。投手と打者でそれぞれ全員合計すると合計出塁数はは同じになるはずなので、そういった関係などを使ってシミュレーションを走らせ、試合の出場メンバー及び打順・先発投手から勝ち負けを予想することができそうです。

また、Yale MBA 時代にこの教授のスポーツアナリティクスの授業を受けておもしろかったので、ここに書かれている分析などもやってみたいです。 Scorecasting

清原氏「我が投資術」を検証してみる

2025/10/25 修正- コードブロックとテーブルの見た目を修正しました

個別株投資

金融関係のお勉強をした人はみんな、分散投資が好きです。私自身も個別株投資反対派です。日本の会社で働いていたときは、SP500 ETF と自社株に全部投資していました。先日、こちらのブログにて個別株投資の本が紹介されていました。この業界では有名な清原氏の著作です。 https://shoshimin.hatenablog.com/entry/2024/03/16/104359

このブログに書かれている通り、基本的には、儲かりそうな個別株が本当にあったとしても、自分には見つけることはできないと思います。自分より賢く、有利な人たちが必死に探しているなか、どうして自分が見つけられると思うのでしょう?

とはいえ、企業の特徴量をもとにポートフォリオを構築するというのは、ファクター投資の手法ですね。清原氏の投資術の中には、ファクター投資の手法を使って検証できそうな部分があります。

私自身はこの手法に従うことはありませんが、実際に有効なのか、データを使って検証してみます。データには J-Quant を利用しました。J-Quant の使い方については、以下の記事が参考になります。 https://zenn.dev/morim34/articles/ae053ce5b94419

清原氏の投資術におけるスクリーニング条件

上述のブログ記事から引用します。

  • PER10 倍以下
  • ネットキャッシュ比率 60%以上
  • 小型株(時価総額 500 億円未満)

これらの条件を見ると、それぞれ Profit、Value、Size ファクターに似ているような気もします。

検証コード

日次の株価データを用意し、毎月末に上記のスクリーニング条件でリバランスを行います。スクリーニング条件に合致したものをロングし、入らないものをショートします。各銘柄のポジションサイズは、時価総額に比例させます。各銘柄のリターンは、日次リターンを月次で合計したものを用います。

ポートフォリオのリターンを計算します。結果だけ見ると、かなりもうかりますね!

# kiyohara parameters
per_bar <- 10 # should be lower than this
net_cash_bar <- 0.6 # should be higher than this
mcap_bar <- 50 * 10^9 # 50 billion yen
# update_data(target_date)
# construct portfolio
dat_port <- read_parquet(file.path(here("data"), "dat_port.parquet")) |>
    drop_na() |>
    group_by(month = ceiling_date(date, "month") - 1) |>
    mutate(flag_rebalance = (date == max(date)))

dat_port_m <- dat_port |>
    group_by(month, code) |>
    summarise(
        ret = sum(ret),
        mcap = last(mcap),
        per = last(per),
        net_cash_ratio = last(net_cash_ratio)
    ) |>
    ungroup()

trading <- dat_port_m |>
    mutate(
        buy_sell_per = if_else((per < per_bar), +1, -1),
        buy_sell_net_cash = if_else((net_cash_ratio > net_cash_bar), +1, -1),
        buy_sell_mcap = if_else((mcap > mcap_bar), +1, -1)
    ) |>
    mutate(
        buy_sell = (buy_sell_per + buy_sell_net_cash + buy_sell_mcap) / 3,
        posi = buy_sell * mcap
    ) |>
    group_by(month) |>
    mutate(wgt = posi / sum(posi)) |>
    ungroup()

result <- trading |>
    summarize(ret = sum(ret * wgt), .by = month) |>
    mutate(cum_ret = exp(cumsum(ret)))
result |>
    ggplot() +
    aes(x = month, y = cum_ret) +
    geom_line() +
    labs(
        title = "Cumulative Return of Kiyohara Strategy",
        x = "Month",
        y = "Cumulative Return"
    )

result |>
    summarize(
        total_return = last(cum_ret) - 1,
        annualized_return = (last(cum_ret))^(12 / n()) - 1,
        annualized_volatility = sd(ret) * sqrt(12),
        sharpe_ratio = annualized_return / annualized_volatility
    ) |>
    mutate(across(where(is.numeric), ~ round(.x * 100, 2))) |>
    knitr::kable(format = "markdown", caption = "Kiyohara Strategy Performance (%)")
total_return annualized_return annualized_volatility sharpe_ratio
81.96 11.69 14.22 82.18

Kiyohara Strategy Performance (%)

リターンが良いのは、隠れたリスクを取っているからでは?

これはあくまでバックテストであり、実際にこの戦略が有効かどうかはわかりません。リターンが良いのは、隠れたリスクを取っているにもかかわらず、リスクが顕在化していなかったのかもしれません。

まずは、Fama-French の 5 ファクターとモメンタムファクターを用いて、このポートフォリオがどの程度のリスクを取っているか検証してみます。

dat_factor <- frenchdata::download_french_data("Fama/French Japanese 5 Factors")$subsets$data[[1]] |>
    rename_all(str_to_lower) |>
    rename(mkt_rf = `mkt-rf`) |>
    mutate(date = ym(date) |> ceiling_date("month") - 1) |>
    rename(month = date)
dat_wml <- frenchdata::download_french_data("Japanese Momentum Factor (Mom)")$subsets$data[[1]] |>
    rename_all(str_to_lower) |>
    mutate(date = ym(date) |> ceiling_date("month") - 1) |>
    rename(month = date)
result_ff <- result |>
    select(-cum_ret) |>
    mutate(ret = ret * 100) |>
    left_join(dat_factor, by = join_by(month)) |>
    left_join(dat_wml, by = join_by(month))

# Fama-MacBeth regression
model <- lm(ret ~ mkt_rf + smb + hml + rmw + cma + wml, data = result_ff)
tidy(model) |>
    mutate(across(where(is.numeric), ~ round(.x, 2))) |>
    knitr::kable(format = "markdown", caption = "Fama-French 5 Factors + Momentum Results")
term estimate std.error statistic p.value
(Intercept) 1.11 0.36 3.07 0.00
mkt_rf 0.52 0.10 5.34 0.00
smb 0.38 0.22 1.74 0.09
hml -0.23 0.20 -1.15 0.25
rmw -0.50 0.37 -1.34 0.18
cma -0.64 0.33 -1.91 0.06
wml 0.11 0.14 0.82 0.42

Fama-French 5 Factors + Momentum Results

glance(model) |>
    mutate(across(where(is.numeric), ~ round(.x, 2))) |>
    knitr::kable(format = "markdown", caption = "Fama-French 5 Factors + Momentum information")
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.59 0.55 2.75 14.02 0.00 6 -154.39 324.77 342.17 440.04 58 65

Fama-French 5 Factors + Momentum information

円建てのポートフォリオのリターンをドル建てリターンで回帰してしまっているので説明力が低いですね。。。しかもこの期間は円安が進んだので、ドル建てリターン(ファクターの方)はリターンが下振れしており、切片が上振れしているように思います。

とはいえ、マーケットファクター以外に優位な係数はなく、強いていうならば Size(SMB)と Investment(CMA)でしょうか。小型株と積極投資株によっていることがわかります。また、もっと Value や Profit ファクターによっているものかと思いましたが、むしろ逆ですね。

やはり、円建て vs ドル建ての検証なのでノイズが大きく、この結果ではなんとも言えなさそうです。

清原氏の投資術をファクターモデルに組み込む

清原氏の投資術が有効である(=株式の値付けの仕組みを解き明かしている)ならば、Fama-French の 5 ファクターとモメンタムファクターに加えて、清原氏の投資術によるリターンを説明変数に加えることで、株式リターンの説明力が上がるはずです。Fama-MacBeth 回帰で検証してみます。

dat_returns <- trading |>
    select(month, code, ret, mcap) |>
    left_join(result_ff |> rename(kiyohara = ret), by = join_by(month))

ts_regression <- dat_returns |>
    select(-rf) |>
    group_nest(code) |>
    mutate(
        model = map(data, ~ lm(ret ~ kiyohara + mkt_rf + smb + hml + rmw + cma + wml, data = .x)),
        beta = map(model, ~ tidy(.x) |> select(term, estimate))
    ) |>
    unnest(beta)

cx_regression <- trading |>
    select(month, code, ret) |>
    left_join(ts_regression |> select(code, term, estimate), by = join_by(code)) |>
    pivot_wider(names_from = term, values_from = estimate) |>
    group_nest(month) |>
    mutate(
        model = map(data, ~ lm(ret ~ kiyohara + mkt_rf + smb + hml + rmw + cma + wml, data = .x)),
        gamma = map(model, ~ tidy(.x) |> select(term, estimate))
    ) |>
    unnest(gamma)
cx_regression |>
    group_by(term) |>
    summarise(
        mean = mean(estimate),
        sd = sd(estimate),
        t_stat = mean / (sd / sqrt(n())),
        p_value = 2 * pt(-abs(t_stat), df = n() - 1)
    ) |>
    knitr::kable(
        format = "markdown",
        caption = "Fama-French 5 Factors + Momentum + Kiyohara Strategy Fama-MacBeth Regression Results",
        digits = 3)
term mean sd t_stat p_value
(Intercept) 0.004 0.011 2.661 0.010
cma 0.120 2.131 0.452 0.653
hml 0.448 3.807 0.949 0.346
kiyohara -0.438 4.163 -0.849 0.399
mkt_rf -0.100 4.243 -0.189 0.851
rmw -0.159 2.032 -0.630 0.531
smb -0.226 1.798 -1.014 0.314
wml 0.273 2.826 0.779 0.439

Fama-French 5 Factors + Momentum + Kiyohara Strategy Fama-MacBeth Regression Results

係数がどれも有意ではないですし切片も有意にゼロではないので、個別株のリターンを説明する上ではこのファクターは役に立たないようですね。正直、ドル建てなのでノイズが大きく、あまり意味のある検証にはなっていない気もします。とはいえ、手続き的には最後に切片の有意性を$\chi2$検定(つまり、GRS 検定)で確認するべきだと思うので、念の為最後に実施してこの記事をまとめます。

面倒だったので ChatGPT 先生に GRS 検定のコードを書いてもらいました。中身をよく考えずに丸投げしたのは今回が初めてだったのですが、逆行列問題(多重共線性問題)や変数の数の問題にもうまく対応してくれていて、質の高さにびっくりしています。ただ、コードの書き方にクセがあってわかりにくいですね。不要な変数を作りがちだったり、データ構造に依存しないコードの書き方になっていたり、2020 年以前の R コードの雰囲気を感じます。

# ---- 1) GRS 本体(前回の最小実装と同じ) ----
grs_test <- function(R, F, rf = NULL, epsilon = 1e-8) {
  R <- as.matrix(R); F <- as.matrix(F)
  stopifnot(nrow(R) == nrow(F))
  Tt <- nrow(R); N <- ncol(R); K <- ncol(F)
  if (!is.null(rf)) {
    rf <- as.numeric(rf); stopifnot(length(rf) == Tt)
    R <- R - rf
    F <- F - rf
  }
  X <- cbind(1, F)
  Bhat <- solve(t(X) %*% X, t(X) %*% R)
  alpha_hat <- Bhat[1, , drop=FALSE]
  E  <- R - X %*% Bhat
  Sigma_hat <- crossprod(E) / (Tt - K - 1)
  # Add small constant to diagonal to handle near-singularity
  diag(Sigma_hat) <- diag(Sigma_hat) + epsilon
  fbar <- colMeans(F)
  Sf   <- stats::var(F)
  c_term <- as.numeric(1 + t(fbar) %*% solve(Sf) %*% fbar)
  A <- as.numeric(alpha_hat %*% solve(Sigma_hat) %*% t(alpha_hat))
  grs <- ((Tt - N - K) / N) * (A / c_term)
  list(
    statistic = unname(grs),
    df = c(df1 = N, df2 = Tt - N - K),
    p.value = 1 - stats::pf(grs, N, Tt - N - K),
    alpha = drop(alpha_hat),
    Sigma = Sigma_hat,
    fbar = fbar
  )
}

# ---- 2) ロング型 → 行列作成 → GRS を一発で ----
# dat: ロング型(例の dat_returns)
# time_col, id_col, ret_col, rf_col, factor_cols を列名で指定
# scale_ret: ret が % 表記なら 100、実数なら 1
grs_from_long <- function(
  dat,
  time_col   = month,
  id_col     = code,
  ret_col    = ret,
  rf_col     = kiyohara,
  factor_cols = c("mkt_rf","smb","hml","rmw","cma","wml"),
  scale_ret  = 100,
  max_assets = 30  # New parameter to limit number of assets
) {
  time_col <- enquo(time_col)
  id_col   <- enquo(id_col)
  ret_col  <- enquo(ret_col)
  rf_col   <- enquo(rf_col)

  # 1) 必要列だけ & ソート
  df <- dat %>%
    select(!!time_col, !!id_col, !!ret_col, !!rf_col, all_of(factor_cols)) %>%
    arrange(!!time_col, !!id_col)

  # 2) リターンを超過化(ret - rf)。% → 実数化もここで実施
  df <- df %>%
    mutate(
      ret_excess = (!!ret_col)/scale_ret - (!!rf_col)/scale_ret
    )

  # 3) ファクターは時点×Kのユニーク行(銘柄に重複しているので distinct)
  F_mat <- df %>%
    select(!!time_col, all_of(factor_cols)) %>%
    distinct() %>%
    arrange(!!time_col)

  # 4) アセットの超過リターンをワイド化(時点×銘柄)
  R_wide <- df %>%
    select(!!time_col, !!id_col, ret_excess) %>%
    pivot_wider(names_from = !!id_col, values_from = ret_excess) %>%
    arrange(!!time_col)

  # 5) 時系列を揃える(F, R_wide ともに同じ時点だけを残す)
  common_dates <- intersect(pull(F_mat, !!time_col), pull(R_wide, !!time_col))
  F_mat  <- F_mat  %>% filter(!!time_col %in% common_dates) %>% arrange(!!time_col)
  R_wide <- R_wide %>% filter(!!time_col %in% common_dates) %>% arrange(!!time_col)

  # 6) Select top assets by average market cap
  avg_mcap <- dat %>%
    group_by(!!id_col) %>%
    summarise(avg_mcap = mean(mcap, na.rm = TRUE)) %>%
    arrange(desc(avg_mcap)) %>%
    slice_head(n = max_assets) %>%
    pull(!!id_col)

  # 7) 欠損のある銘柄は落とす(完全観測の列のみ採用)
  R_mat_all <- R_wide %>% select(-!!time_col)
  keep_cols <- names(R_mat_all)[colSums(is.na(R_mat_all)) == 0]
  keep_cols <- intersect(keep_cols, as.character(avg_mcap))  # Only keep top assets by market cap
  R_mat <- R_mat_all %>% select(all_of(keep_cols)) %>% as.matrix()

  # 8) F 行列
  F_m <- F_mat %>% select(all_of(factor_cols)) %>% as.matrix()

  # 8) GRS 実行
  out <- grs_test(R = R_mat, F = F_m)

  # 9) 付帯情報(どの銘柄が入ったか等)を整形して返す
  list(
    grs = tibble(
      statistic = out$statistic,
      df1 = out$df["df1"],
      df2 = out$df["df2"],
      p.value = out$p.value,
      T = nrow(F_m),
      N = ncol(R_mat),
      K = ncol(F_m)
    ),
    alpha = tibble(code = keep_cols, alpha = out$alpha),
    info = list(
      used_codes = keep_cols,
      used_dates = common_dates,
      factors = factor_cols
    )
  )
}



res <- grs_from_long(
  dat_returns,
  time_col = month,
  id_col   = code,
  ret_col  = ret,
  rf_col   = rf,
  factor_cols = c("mkt_rf","smb","hml","rmw","cma","wml", "kiyohara"),  # 必要に応じて増減
  scale_ret = 100  # ret, rf が%表示なら100。実数なら1
)

res$grs |>     # GRS統計量・自由度・p値など(1行tibble)
knitr::kable(format = "markdown", caption = "GRS Test Including Kiyohara Factor Results")
statistic df1 df2 p.value T N K
11.62003 4 54 7e-07 65 4 7

GRS Test Including Kiyohara Factor Results

res$alpha  |>  # 各銘柄のα推定値(超過月次リターンの実数単位)
    ggplot() + aes(x =  alpha) +
    geom_histogram(bins = 30) +
        labs(
            title = "Estimated Alphas from GRS Test Including Kiyohara Factor",
            x = "Estimated Alpha",
            y = "Frequency"
        )

ということで、p.value が非常に小さいので、「すべての銘柄でα=0」という帰無仮説は棄却されます。つまり、このファクターモデルは銘柄をうまくプライシングできておらず、清原氏の投資術は市場のクロスセクション方向においてミスプライシングしているようです。

まとめ・今後の課題

お手元の手軽な検証では、清原氏の投資術に基づくポートフォリオはバックテスト上では高いリターンを示すことがわかりました。しかし、ファクターモデルを用いた検証ではノイズが大きすぎ、GRS 検定では有意にミスプライシングしていることがわかりました。やはり、個人的にはこの要件を使って投資しようとは思わないですね。

FF ファクターがドル建てなのでノイズが大きすぎて、後半の検証はほとんど意味がないようにも思えます。今後は、J-Quant のデータを使って円建てファクターを作成し、検証してみたいところです。

給料は何で決まるのか?

給料は何できまるのか

給料を上げたい!

僕も含め、多くの人が毎日(ではなくとも毎月は)思っているであろうことですね。思っていない人はこの記事を読むのは時間の無駄なのでここで読むのを止めてください。 研究活動の関係で、東洋経済データサービスの企業データを手に入れたので、給与関連情報を使って給料の高い会社・高い理由を探っていきましょう。私自身が30代ということもあり、30歳時点での給与に強い興味があるのでまずはそこから分析していきます。

https://biz.toyokeizai.net/data/

30歳時点での平均給与分析

データ分析のセオリーに従い、まずはデータの内容を把握しましょう。時系列で見てみたり、ランキングをつけてみます。

ここ数年を見ると、30歳時点での平均給与はなんとなく上昇傾向に見えますね。

# dir_path: 東洋経済のデータへのパス
dat_social <-
  list.files(dir_path, full.names = TRUE) |> 
  map(~ read_csv(.,
                 locale = locale(encoding = "CP932"),
                 show_col_types = FALSE)
      ) |>
  bind_rows() |> 
  select(-contains("その他"),
         -starts_with("注記"),
         -starts_with("概要")) |> 
  arrange(年版)

dat_target <- 
  dat_social |>  
  select(fy = 年版, code = 証券コード,
         company = 社名,
         industry = 業種名,
         leave_male = `離職者・男性`,
         leave_female = `離職者・女性`,
         leave_total = `離職者・合計`,
         salary = `平均年間給与/前期`, 
         age_male = `平均年齢男性/前期`,
         age_female = `平均年齢女性/前期`,
         salary_30y = `30歳平均賃金`,
         salary_30y_max = `30歳最高賃金`,
         salary_30y_min = `30歳最低賃金`,
         hours_annual = 年間総労働時間,
         hire_grad = `新卒採用・今年4月入社・全体合計`,
         hire_exp = `中途採用・昨年4月~今年3月入社・全体合計`,
         employees = `従業員数合計/前期`
         ) |> 
  mutate(hours_annual = if_else(fy == 2020 & code == 4004,
                                hours_annual/10,
                                hours_annual)) |> 
  drop_na()

dat_target |> 
  group_by(fy) |> 
  summarise(salary_30y = mean(salary_30y),
            salary_30y_max = mean(salary_30y_max),
            salary_30y_min = mean(salary_30y_min)) |> 
  ggplot() + aes(x = fy, y = salary_30y ) +
  geom_line() +
  geom_errorbar(aes(ymin = salary_30y_min, ymax = salary_30y_max, width = .1)) +
  scale_y_continuous(labels = scales::comma, limits = c(0, 5 * 10^5)) +
  labs(x = "", y = "", title = "30歳時点での平均給与・最大額の平均・最低額の平均")

対象期間の平均給与でランキングをつけてみました。金額感的に、これはどうやら、月次かつボーナスが入っていない基本給のようですね。とはいえ、博報堂、野村HD、商社や海運会社など、高給でよく知られる会社がでてきます。

dat_target |> 
  mutate(salary_30y = if_else(salary_30y > 10^6, salary_30y/10, salary_30y),
         salary_30y_max = if_else(salary_30y_max > 10^6, salary_30y_max/10, salary_30y_max)) |> 
  group_by(company) |> 
  summarise(salary_30y = mean(salary_30y),
            salary_30y_max = mean(salary_30y_max)) |> 
  arrange(desc(salary_30y)) |> head(20) |> 
  kable(digits = c(0,0),
        format.args = list(big.mark = ","))
company salary_30y salary_30y_max
大成温調 591,773 691,487
博報堂DYホールディングス 547,021 643,333
MS&Consulting 542,400 663,500
テノックス 526,580 567,292
GA technologies 520,804 177,000
ケー・エフ・シー 507,693 615,914
マクアケ 501,000 760,000
商船三井 493,200 511,827
住友商事 491,400 541,200
アイエックス・ナレッジ 490,891 603,587
LINE 483,750 650,000
第一生命ホールディングス 480,457 648,093
大和証券グループ本社 478,666 647,850
フォーバル 477,940 503,659
三菱商事 474,587 552,500
岩崎通信機 453,742 492,604
双日 449,885 543,680
ディジタルメディアプロフェッショナル 438,250 489,067
野村ホールディングス 437,676 317,692
ストライク 425,836 484,200

上述の通り、これは基本給なのでボーナスの比率が大きい会社は実態を表してなさそうです。商社では給料の半分以上がボーナスでしたし。。。ここからは、ボーナスも含んだ金額が開示されているはずの全社平均給与を対象にして分析します。

全社平均給与の分析

平均給与を見ると、上昇傾向はあまり見られないですね。。。

dat_target |> 
  group_by(fy) |> 
  summarise(ave_salary = mean(salary),
            salary_sd = sd(salary))|> 
  ggplot() + aes(x = fy, y = ave_salary ) +
  geom_line() +
  geom_errorbar(aes(ymin = ave_salary - salary_sd, ymax = ave_salary + salary_sd, width = .1)) +
  labs(x = "", y = "", title = "平均給与±標準偏差") +
  scale_y_continuous(labels = scales::comma, limits = c(0, 9*10^6))

30歳時点と同じように、高い順に並べてみましょう。商社がツートップですね。三井物産伊藤忠・丸紅が出てこない理由が少し気になりますが。。。

dat_target |> 
  group_by(company) |> 
  summarise(salary = mean(salary)) |> 
  arrange(desc(salary)) |> head(20) |> 
  kable(digits = 0,
        format.args = list(big.mark = ","))
company salary
三菱商事 16,150,024
住友商事 13,688,822
ストライク 13,506,500
電通 12,477,350
電通グループ 12,284,489
野村総合研究所 12,052,333
博報堂DYホールディングス 11,094,068
第一三共 11,080,665
双日 11,032,042
豊田通商 10,768,813
大和証券グループ本社 10,613,868
大気社 10,560,577
WOWOW 10,484,500
野村ホールディングス 10,452,100
中外製薬 10,274,714
商船三井 10,001,674
野村不動産ホールディングス 9,939,427
大成建設 9,912,430
アドバンテスト 9,831,321
みずほフィナンシャルグループ 9,822,833

業界で給与が決まる?

とはよく言われますがので、業界別に給与を見ていきます。業界内での標準誤差は低く見えますが、よくよく見ると、下3分の一と中3分の一の中では平均の差に有意差がなさそうです。業界内でのばらつきが大きいからですね。高給ネタで必ず出てくる総合商社は卸売業に入りますし、三菱地所などの不動産業も意外と上位にこないので、わりと例外的な存在ということがわかります。平均的に高いのは海運、証券、鉱業や石油などで、このあたりは確かに「業界が給与を決める」と言ってもよさそうです。

よく知られた話ですが、米国では圧倒的な高給を誇る情報・通信業が低いのはなんだか不思議ですよね。この時代に。。。

dat_target |> 
  group_by(industry) |>
  summarise(ave_salary = mean(salary),
            se = sd(salary)/sqrt(n())) |> 
  ggplot() + aes(x = reorder(x = industry, X = ave_salary)) +
  geom_bar(aes(y = ave_salary), stat = "identity") +
  geom_errorbar(aes(ymax = ave_salary + se, ymin = ave_salary - se)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "", y = "", title = "平均年収および標準誤差") 

また、回帰分析にて確認してみると、いくつかの業界のダミー変数は統計的に有意なので、その業界については「業界で給料が決まる」と言ってよさそうです。例えば証券・石油・金属製品、陸運などですね。

lm(salary~industry, data = dat_target) |> 
  tidy() |> 
  kable(digits = c(0,0,0,2,2), format.args = list(big.mark = ","))
term estimate std.error statistic p.value
(Intercept) 7,263,126 189,747 38.28 0.00
industryゴム製品 -1,007,237 265,591 -3.79 0.00
industryサービス業 -1,472,662 214,516 -6.87 0.00
industryその他製品 -823,156 237,993 -3.46 0.00
industryその他金融業 -254,135 269,798 -0.94 0.35
industryパルプ・紙 -583,346 314,130 -1.86 0.06
industry不動産業 70,519 285,900 0.25 0.81
industry保険業 -1,223,889 307,174 -3.98 0.00
industry倉庫・運輸関連業 -914,539 420,738 -2.17 0.03
industry化学 -152,239 205,544 -0.74 0.46
industry医薬品 1,241,478 252,478 4.92 0.00
industry卸売業 16,034 211,702 0.08 0.94
industry小売業 -1,669,603 214,108 -7.80 0.00
industry建設業 719,364 214,937 3.35 0.00
industry情報・通信業 -661,024 211,921 -3.12 0.00
industry機械 -297,344 212,258 -1.40 0.16
industry水産・農林業 339,269 435,706 0.78 0.44
industry海運業 2,738,549 563,947 4.86 0.00
industry石油・石炭製品 1,401,947 527,016 2.66 0.01
industry空運業 -443,279 563,947 -0.79 0.43
industry精密機器 -314,144 281,769 -1.11 0.26
industry繊維製品 -1,274,301 263,029 -4.84 0.00
industry証券、商品先物取引 3,354,519 563,947 5.95 0.00
industry証券・商品先物 1,852,047 420,738 4.40 0.00
industry輸送用機器 -825,932 215,226 -3.84 0.00
industry金属製品 -1,048,838 257,336 -4.08 0.00
industry鉄鋼 -886,931 326,362 -2.72 0.01
industry鉱業 1,956,470 1,314,609 1.49 0.14
industry銀行業 379,046 310,547 1.22 0.22
industry陸運業 -847,511 285,900 -2.96 0.00
industry電気・ガス業 576,623 677,534 0.85 0.39
industry電気機器 -178,898 204,848 -0.87 0.38
industry非鉄金属 -458,276 253,392 -1.81 0.07
industry食料品 -290,545 217,280 -1.34 0.18

なぜ、業界で給与が決まるのか?

ということで、「高い給与を得るためには給与の高い業界に入るか給与の高い会社にはいる」という身も蓋もない結論が得られてしまったのですが、それではおもしろくないのでもう一歩踏み込んでみます。なぜ業界で給与に違いがうまれるのでしょうか?

よく話題に登る仮説をいくつか挙げた上で検証してみます。 - 労働時間が違う - 労働分配率が高い(設備投資が不要で費用における人件費の割合が高い)

労働時間

業界ごとに、平均労働時間を横軸、平均給与を縦軸にして散布図を描いてみましたが、全然関係なさそうですね。。。

dat_target |> 
  group_by(industry) |> 
  summarise(ave_work = mean(hours_annual),
            se = sd(hours_annual)/sqrt(n()),
            ave_salary = mean(salary)) |> 
  ggplot() + aes(x = ave_work, y = ave_salary) +
  geom_point() +
  geom_smooth(method = "lm" )
`geom_smooth()` using formula = 'y ~ x'

ちなみに、業界と平均労働時間を交差項とした回帰分析の結果はこちらです。ほとんどの業界で有意な効果は見られませんでした。長く働くことに意味はないと。。。?令和っぽいですね。

lm(salary~hours_annual * industry, data = dat_target) |> 
  tidy() |> 
  kable(digits = c(0,0,0,2,2), format.args = list(big.mark = ","))
term estimate std.error statistic p.value
(Intercept) 9,549,113 3,559,003 2.68 0.01
hours_annual -1,167 1,814 -0.64 0.52
industryゴム製品 -6,831,435 4,526,935 -1.51 0.13
industryサービス業 -2,554,652 3,659,612 -0.70 0.49
industryその他製品 1,201,953 4,238,259 0.28 0.78
industryその他金融業 1,739,996 4,305,483 0.40 0.69
industryパルプ・紙 6,143,304 6,332,777 0.97 0.33
industry不動産業 -5,931,649 4,635,521 -1.28 0.20
industry保険業 -9,714,217 7,590,638 -1.28 0.20
industry倉庫・運輸関連業 -7,017,403 7,709,819 -0.91 0.36
industry化学 -4,202,161 3,766,416 -1.12 0.26
industry医薬品 -2,379,256 5,016,693 -0.47 0.64
industry卸売業 -6,754,739 3,836,077 -1.76 0.08
industry小売業 -2,709,238 3,682,057 -0.74 0.46
industry建設業 -6,166,105 3,950,943 -1.56 0.12
industry情報・通信業 -3,281,993 3,905,310 -0.84 0.40
industry機械 -2,028,945 3,934,614 -0.52 0.61
industry水産・農林業 -3,762,085 11,699,980 -0.32 0.75
industry海運業 -12,857,786 25,907,976 -0.50 0.62
industry石油・石炭製品 -13,599,715 38,660,008 -0.35 0.73
industry空運業 2,993,346 50,091,997 0.06 0.95
industry精密機器 -2,467,099 5,207,777 -0.47 0.64
industry繊維製品 -11,619,687 5,330,193 -2.18 0.03
industry証券、商品先物取引 1,256,336 10,013,076 0.13 0.90
industry証券・商品先物 -3,320,686 10,113,171 -0.33 0.74
industry輸送用機器 -1,986,864 3,938,208 -0.50 0.61
industry金属製品 -4,711,350 5,137,677 -0.92 0.36
industry鉄鋼 -10,925,031 6,209,315 -1.76 0.08
industry鉱業 2,002,642 1,311,534 1.53 0.13
industry銀行業 -15,176,825 6,530,912 -2.32 0.02
industry陸運業 -1,006,942 4,229,554 -0.24 0.81
industry電気・ガス業 -686,324 12,847,035 -0.05 0.96
industry電気機器 -2,396,426 3,680,756 -0.65 0.52
industry非鉄金属 1,668,785 5,313,531 0.31 0.75
industry食料品 -3,342,836 3,868,759 -0.86 0.39
hours_annual:industryゴム製品 2,931 2,287 1.28 0.20
hours_annual:industryサービス業 574 1,862 0.31 0.76
hours_annual:industryその他製品 -1,018 2,156 -0.47 0.64
hours_annual:industryその他金融業 -1,001 2,188 -0.46 0.65
hours_annual:industryパルプ・紙 -3,443 3,233 -1.06 0.29
hours_annual:industry不動産業 2,975 2,317 1.28 0.20
hours_annual:industry保険業 4,534 4,063 1.12 0.26
hours_annual:industry倉庫・運輸関連業 3,134 3,959 0.79 0.43
hours_annual:industry化学 2,086 1,924 1.08 0.28
hours_annual:industry医薬品 1,882 2,621 0.72 0.47
hours_annual:industry卸売業 3,479 1,958 1.78 0.08
hours_annual:industry小売業 554 1,872 0.30 0.77
hours_annual:industry建設業 3,333 1,986 1.68 0.09
hours_annual:industry情報・通信業 1,338 1,991 0.67 0.50
hours_annual:industry機械 886 2,003 0.44 0.66
hours_annual:industry水産・農林業 2,084 5,914 0.35 0.72
hours_annual:industry海運業 7,614 12,559 0.61 0.54
hours_annual:industry石油・石炭製品 7,783 20,110 0.39 0.70
hours_annual:industry空運業 -1,607 24,287 -0.07 0.95
hours_annual:industry精密機器 1,097 2,692 0.41 0.68
hours_annual:industry繊維製品 5,525 2,808 1.97 0.05
hours_annual:industry証券、商品先物取引 1,072 5,047 0.21 0.83
hours_annual:industry証券・商品先物 2,657 5,209 0.51 0.61
hours_annual:industry輸送用機器 617 1,992 0.31 0.76
hours_annual:industry金属製品 1,839 2,560 0.72 0.47
hours_annual:industry鉄鋼 5,062 3,132 1.62 0.11
hours_annual:industry鉱業 NA NA NA NA
hours_annual:industry銀行業 7,919 3,323 2.38 0.02
hours_annual:industry陸運業 149 2,116 0.07 0.94
hours_annual:industry電気・ガス業 642 6,585 0.10 0.92
hours_annual:industry電気機器 1,132 1,878 0.60 0.55
hours_annual:industry非鉄金属 -1,027 2,671 -0.38 0.70
hours_annual:industry食料品 1,562 1,975 0.79 0.43

労働分配率

元経済学者の卵としては、なんらかのかたちで労働分配率が関係しているはずだと信じています。以下のようなコブダグラス型生産関数を想定すると、均衡状態では労働分配率と設備分配率はそれぞれの相対的な生産性に収束するはずです。つまり、設備投資がより必要な製造業では労働分配率が低くなるはずです。逆に、労働生産性が設備生産性に対して高い金融業が高くなるはずということです。

Y = A**LαK1 − α

ということで、別のプロジェクトで使った財務データを借用して、実際にどんなもんかと見てみました。稼いだキャッシュに占める設備投資の割合をみるため、減価償却費/Gross Profitを見ていきます。本当は、CAPEX/Gross Profitを見たいところですが、CAPEXの情報がなかったので、均されているはずの減価償却費を使用しました。

業界ごとに箱ヒゲ図を描いて見ると、想定通り、かなりばらつきがありますね。証券先物が低めなのはイメージ通りですが、権益への投資が大きいはずの石油・石炭が低いのはなぜなんでしょうか。。。?おそらく、IFRS基準であり、投資に伴うのれんが償却されていないからということでしょう。この点は、純粋な投資業の減価償却費が低いことで、想定どおりの計測基準になっているようです。

dat_financials <- 
  read_csv(fin_path) |> 
  mutate(date = as.Date(date), 
         code = str_replace(code, '.T', '') |> as.numeric())
target_name <- 
  c('Total Revenue',
    'Gross Profit',
    'EBITDA',
    'EBIT'
  )

fin_filtered <- 
  dat_financials |> 
  filter(name %in% target_name) |> 
  pivot_wider(names_from = name, values_from = value) |> 
  rename(gross_profit = 'Gross Profit',
         total_revenue = "Total Revenue" ) |> 
  filter(!is.na(code)) |> 
  unnest(everything()) |> 
  mutate(amoti = EBITDA - EBIT,
         invest_ratio = amoti/gross_profit) |> 
  mutate(fy = year(date)) |> 
  select(fy, code, invest_ratio)

dat_invest <- 
  dat_target |> 
  filter(fy > 2021) |> 
  left_join(fin_filtered, by = c("fy", "code")) |> 
  filter(between(invest_ratio, 0,1))

dat_invest |> 
  group_by(industry) |> 
  mutate(ave_salary = mean(salary)) |> 
  ggplot() + aes(x = reorder(x = industry, X = ave_salary), y = invest_ratio) + 
  geom_boxplot() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "業界(平均年収低い順)", y = "設備投資比率")

さて、それでは設備投資比率と平均年収の関係性を見てみましょう。横軸に設備投資比率、縦軸に平均年収をとって散布図を描いてみました。設備投資比率が高い業界は年収が低いということがよくわかります。これは、設備投資比率が高い業界は、設備投資による生産性向上が労働生産性向上よりも大きいため、労働分配率が低くなるということを示しています。逆に、労働生産性向上が大きい業界は、労働分配率が高くなるため、平均年収が高くなるということですね。

理論がデータで実証されるとうれしいですね。

dat_invest |> 
  group_by(industry) |> 
  summarise(invest_ratio = mean(invest_ratio), 
            ave_salary = mean(salary)) |> 
  filter(invest_ratio != max(invest_ratio)) |> 
  ggplot() +
  aes(x = invest_ratio, y = ave_salary) +
  geom_point() +
  geom_smooth(method = "lm") + 
  scale_y_continuous(labels = scales::comma)
`geom_smooth()` using formula = 'y ~ x'

なお、海運業については外れ値となっていたため、この分析からは除外しました。おそらく船が高いせいで投資比率がめちゃくちゃ高い一方で、業界の特殊事情により平均年収も高い、ということかと。

dat_invest |> 
  group_by(industry) |> 
  summarise(invest_ratio = mean(invest_ratio), 
            ave_salary = mean(salary)) |> 
  filter(invest_ratio == max(invest_ratio))
# A tibble: 1 × 3
  industry invest_ratio ave_salary
  <chr>           <dbl>      <dbl>
1 海運業          0.485   10494227

まとめ

ということで、給料を上げるには業界と会社に気をつけること、そして業界で給料が決まる理由は設備投資比率ということがわかりました。

PythonでRのgroup_nestを実現する方法

PythonのPolarsでRのgroup_nestを実現する

Polarsでgroup_nestがしたい

R or Python

僕がRで分析するときはいつもこちらの記事を参考にしています。 https://qiita.com/kilometer/items/7184904765fbf0f33f04

tidyverseでデータを読み込んだ後は、なんらかの変数でグループ化して、ネストして、各データセットに対してmapを通じて処理をかけていくと。Rはスクリプト言語でありベクトル化された処理に最適化されているので、大容量のデータを扱う場合はこちらの方が速いですし、何より読みやすいです。

Rはユーザーの活動が活発でTokyoRとかのカンファレンスも多いですし、ネット上でも色々な記事があるので僕にとっては第一(プログラミング)言語なのですが、どうも、金融業界ではそういうわけではないらしいです。Yale大学でも、ファイナンスの教授はみんなPython。一人、統計学出身ぽいBryan Kellyという教授はかつてRを使っていたようなのですが、授業ではPythonを使うよう指定されました。また、彼も含め教授陣の何人かが副業で働いているAQRというヘッジファンドでも、システムはすべてPythonで構築されているそうです(どっちが本業なんでしょうか。。。?)。

個人的には、システムの信頼性やスピードなどを考えると本番トレーディングのシステムはC++かRustで書かれているのではないかと推測していますが、 少なくとも、トレーディング戦略を構築するリサーチ環境はPythonがほとんどだと思います。 と、思いましたがトレーディングシステムを開発するエンジニアもPythonが技術要件になっているので本当にPythonなのかも。。。求人(11/17アクセス)

ちなみに、弊社は偉い人たちの意向でFortranC#が使われています。この業界ではだいぶ珍しいようです。僕としては業界標準に近づけたいので、新しいシステム構築やレポート作成時にはpythonを使うようにしています。

Python: Polars

というわけで、就職してからはPythonに乗り換えようとしてきたわけですが、Pythonはやはりデータ分析というよりはシステム構築がメインなのか、使い勝手が悪いです。その中でも、Polarsというパッケージがtidyverseに似た操作感なので使ってみました。

簡単な使い方は様々なサイトで解説されているので、日本語の記事がほとんどない、group_nestについて記事を記載します。

簡単に解説すると、group_byでグルーピングし、aggにてpl.listとpl.struct型の組み合わせに変換することでデータフレームっぽく扱います。

from palmerpenguins import load_penguins
import polars as pl
df = pl.from_pandas(load_penguins())
# Grouping by code
df.group_by("species").agg(pl.all())
species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g sex year
str list[str] list[f64] list[f64] list[f64] list[f64] list[str] list[i64]
"Gentoo" ["Biscoe", "Biscoe", … "Biscoe"] [46.1, 50.0, … 49.9] [13.2, 16.3, … 16.1] [211.0, 230.0, … 213.0] [4500.0, 5700.0, … 5400.0] ["female", "male", … "male"] [2007, 2007, … 2009]
"Chinstrap" ["Dream", "Dream", … "Dream"] [46.5, 50.0, … 50.2] [17.9, 19.5, … 18.7] [192.0, 196.0, … 198.0] [3500.0, 3900.0, … 3775.0] ["female", "male", … "female"] [2007, 2007, … 2009]
"Adelie" ["Torgersen", "Torgersen", … "Dream"] [39.1, 39.5, … 41.5] [18.7, 17.4, … 18.5] [181.0, 186.0, … 201.0] [3750.0, 3800.0, … 4000.0] ["male", "female", … "male"] [2007, 2007, … 2009]
nested_df = df.group_by("species").agg([pl.struct(pl.all()).alias("nested_data")])
nested_df
species nested_data
str list[struct[8]]
"Adelie" [{"Adelie","Torgersen",39.1,18.7,181.0,3750.0,"male",2007}, {"Adelie","Torgersen",39.5,17.4,186.0,3800.0,"female",2007}, … {"Adelie","Dream",41.5,18.5,201.0,4000.0,"male",2009}]
"Chinstrap" [{"Chinstrap","Dream",46.5,17.9,192.0,3500.0,"female",2007}, {"Chinstrap","Dream",50.0,19.5,196.0,3900.0,"male",2007}, … {"Chinstrap","Dream",50.2,18.7,198.0,3775.0,"female",2009}]
"Gentoo" [{"Gentoo","Biscoe",46.1,13.2,211.0,4500.0,"female",2007}, {"Gentoo","Biscoe",50.0,16.3,230.0,5700.0,"male",2007}, … {"Gentoo","Biscoe",49.9,16.1,213.0,5400.0,"male",2009}]

nestされたデータフレームを元に戻すには、explodeでリストを解消した後、unnestでstruct型を紐解いてあげればOKです。

nested_df.select(pl.col("nested_data")).explode("nested_data").unnest("nested_data")
species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g sex year
str str f64 f64 f64 f64 str i64
"Adelie" "Torgersen" 39.1 18.7 181.0 3750.0 "male" 2007
"Adelie" "Torgersen" 39.5 17.4 186.0 3800.0 "female" 2007
"Adelie" "Torgersen" 40.3 18.0 195.0 3250.0 "female" 2007
"Adelie" "Torgersen" null null null null null 2007
"Adelie" "Torgersen" 36.7 19.3 193.0 3450.0 "female" 2007
"Gentoo" "Biscoe" null null null null null 2009
"Gentoo" "Biscoe" 46.8 14.3 215.0 4850.0 "female" 2009
"Gentoo" "Biscoe" 50.4 15.7 222.0 5750.0 "male" 2009
"Gentoo" "Biscoe" 45.2 14.8 212.0 5200.0 "female" 2009
"Gentoo" "Biscoe" 49.9 16.1 213.0 5400.0 "male" 2009

Palmer Penguinsのインストールに思ったより手間がかかったので、今日はここまで。実際の分析作業は次回の記事でやってみます。

余談

Pythonってめんどくさいですね。。。Rならこんな感じで一瞬なのに。。。

library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(palmerpenguins)

dat_nested  <- 
    penguins |> 
    group_nest(species)

dat_nested
# A tibble: 3 × 2
  species                 data
  <fct>     <list<tibble[,7]>>
1 Adelie             [152 × 7]
2 Chinstrap           [68 × 7]
3 Gentoo             [124 × 7]
dat_nested |> 
    unnest(data)
# A tibble: 344 × 8
   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
   <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
 1 Adelie  Torgersen           39.1          18.7               181        3750
 2 Adelie  Torgersen           39.5          17.4               186        3800
 3 Adelie  Torgersen           40.3          18                 195        3250
 4 Adelie  Torgersen           NA            NA                  NA          NA
 5 Adelie  Torgersen           36.7          19.3               193        3450
 6 Adelie  Torgersen           39.3          20.6               190        3650
 7 Adelie  Torgersen           38.9          17.8               181        3625
 8 Adelie  Torgersen           39.2          19.6               195        4675
 9 Adelie  Torgersen           34.1          18.1               193        3475
10 Adelie  Torgersen           42            20.2               190        4250
# ℹ 334 more rows
# ℹ 2 more variables: sex <fct>, year <int>

CarMaxで車を買いたい

車が欲しい

今、夫婦ふたりでLAのSanta Monicaという街に住んでいます。オフィスもSanta Monica内にあるので、(米国ではほぼありえない)徒歩通勤を実現しており、車は夫婦で一台しか買っていません。 とはいえ、LAで車がないのはさすがに厳しいです。ということで、来年あたりにどんな車を買おうか考えているところです。

そこで、価格の観点からどの車を買うべきかを見ていきましょう。米国の中古車販売大手のCarMaxのAPIを活用し、9月頭時点での在庫情報を取得してきました。

データ確認

データを見る前に、当地でよく言われる「日本車は値段が下がりにくい」という話を思い出しておきましょう。これは、日本車の耐久性が高いため、走行距離が伸びても価格が下がりにくいためとか言われています。正直、データを見たことがないのでほんまかな?と疑っていました。今回の分析を通じてこの仮説を検証できればうれしいですね。

まずはデータの全体像を確認します。横軸に走行距離(マイル)、縦軸に価格(ドル)とし、メーカーごとに色分けします。当然、走行距離がながければ値段は下がりますね。すでによく知られている常識がデータで裏付けられると、とても穏やかな気持ちになります。

dat |> 
  rename(price = basePrice) |> 
  ggplot() +
  aes(x = mileage, y = price, color = make) + 
  geom_point() +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels = scales::comma)

ちなみに我が家の一台目は11,000マイル、\$18,000のNissan Rogueを購入したのですが、ボックスプロットを見るとかなり安めの層だったことがわかります。

dat |> 
  rename(price = basePrice) |> 
  ggplot() +
  aes(x = make, y = price, fill = make) +
  geom_boxplot() +
  scale_y_continuous(labels = scales::comma) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1))

走行距離1マイルあたり減価

日本車

もう少しデータを眺めてみましょう。分析対象が多すぎるとごちゃごちゃするので、日本車に限って線形回帰をしてみます。日本車はNissan, Toyota, Subaru, Honda, Lexus, Mazda, Mitsubishi, Infiniti, Acuraの9メーカーです。

ぱっと見で、ベース価格(切片)と走行距離あたりで価格が下がる割合に結構違いがあるようです。それにしても、ToyotaLexusはどの走行距離でも高い金額で安定していますね。

make_JP <- c("Nissan", "Toyota", "Subaru", "Honda", "Lexus",
             "Mazda", "Mitsubishi", "Infiniti", "Acura" )

dat_jp <- 
  dat |> 
  filter(make %in% make_JP)

dat_jp |> 
  ggplot() +
  aes(x = mileage, y = basePrice, color = make) + 
  geom_point() +
  geom_smooth(method = "lm") +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels = scales::comma)

回帰分析

日本車だけでなく、他のメーカーも入れた上で回帰分析の元の数値を見てみましょう。何の条件もない場合、新車のベース価格は\$44,000程度のようです。走行距離1マイルあたり\$0.14 価格が落ちますが、走行距離よりも販売してからの年数の方が大きな影響があるようです。1年で\$1,500も落ちるんですね。型落ちの走行距離短めの車を買うのが良さそうです。

$R2$を見ると、これらの情報でだいたい50%以上の価格変動が説明できるようです。思ったよりモデルの精度がよさそうです。これだけデータの関係性をうまく捉えられているようでしたら、機械学習モデルを活用するともっと精度があがるかもしれませんね。

dat_jp_lm <- 
  dat_jp |> 
  select( vin, basePrice, year, make, model, mileage, storeName, features,
          transferType, exteriorColor, interiorColor, driveTrain, engineType ) |> 
  mutate(features = str_extract_all(features, "(?<=')[^']+(?=')")) |> 
  unnest(features) |> 
  filter(features != ", ") 

target_dat <- dat |> 
       mutate(year_past = 2024 - year) |> 
       select(-c(vin, features, year)) 
# cross term
model <-
  lm(basePrice ~ year_past + make * mileage,
     data = target_dat)

model |> modelsummary::msummary(#output = "markdown",
                                align = "lr",
                                fmt =  \(x){formatC(x, digits = 2, big.mark = ",", format = "f")},
                                coef_omit = "make*",
                                title = "Base parameters")
Base parameters
(Intercept) 44,065.67
(520.25)
year_past -1,515.87
(16.80)
mileage -0.14
(0.01)
Num.Obs. 49042
R2 0.534
R2 Adj. 0.533
AIC 1005489.5
BIC 1006123.1
Log.Lik. -502672.754
F 800.516
RMSE 6842.47

Base parameters

メーカー×走行距離の交差項

金銭負担の観点からのメーカー選択には、同じ走行距離を走った場合にどれだけ価格が安くなるかを分析することが重要です。 交差項を追加して、あるメーカーを選んだ場合に走行距離1マイルあたりの価格低下率がどう影響するかを見てみましょう。

このグラフを見ると、上位陣には日産、マツダ、スバル、三菱、ヒュンダイ、ホンダなどが上がってきます。これらのメーカーは走行距離が伸びても価格が下がりにくいということですね。逆に、BMWアウディメルセデスベンツジャガー、ランドローバーなどは走行距離が伸びると価格が下がりやすいです。ポルシェなどは最悪ですが、ポルシェを買う人は中古価格など気にしないでしょう。

model |> 
  broom::tidy() |> 
  mutate(
    effect = case_when(estimate > 0 ~ "positive",
                       .default = "negative"),
    term = str_remove_all(term, "model|make|mileage") ) |> 
  filter(str_detect(term, ":")) |> 
  ggplot() +
  aes(y = estimate, x = reorder(x = term, X = estimate)) +
  facet_wrap(effect ~ ., scales = "free") +
  geom_point() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(y = "Coefs of cross term", x = "terms")

走行距離&年式

ここまで書いて、「走行距離というより、年式に影響を受けている可能性はないか?」と気付いたので、年式×メーカーの交差項を追加して分析してみます。

ベースとして、価格は\$45,000、1年古くなるごとに\$2,200安くなってしまうということですね。$R2$ adjusted がそこまで変わらないので、やはり情報が重複していたようです。後で見る通り、年式と走行距離の相関係数が高いので多重共線性が気になりますが、検定力に影響があるだけなので、いったん気にせず進めましょう。

model <-
  lm(basePrice ~ make * year_past + make * mileage,
     data = target_dat )
model |> modelsummary::msummary(#output = "markdown",
                                align = "lr",
                                fmt =  \(x){formatC(x, digits = 2, big.mark = ",", format = "f")},
                                coef_omit = "make*",
                                title = "Base parameters")
Base parameters
(Intercept) 45,071.76
(558.22)
year_past -2,228.42
(153.06)
mileage -0.10
(0.01)
Num.Obs. 49042
R2 0.543
R2 Adj. 0.542
AIC 1004571.2
BIC 1005504.0
Log.Lik. -502179.579
F 558.933
RMSE 6774.00

Base parameters

メーカー×年式で見ても、走行距離とそこまで大きな差はありません。一方で、ポルシェやミニが上位(価格が下がりにくい)に来るのは意外ですね。ポルシェは古くなっても価値が下がりにくいと。確かに、古くても価値は相楽なさそうです。

temp_model <- 
  model |> 
  broom::tidy() |> 
  mutate(
    effect = case_when(estimate > 0 ~ "positive",
                       .default = "negative"))
temp_model |> 
  filter(str_detect(term, "year")) |> 
  mutate(term = str_remove_all(term, "make|year_past")) |> 
  ggplot() +
  aes(y = estimate, x = reorder(x = term, X = estimate)) +
  facet_wrap(effect ~ ., scales = "free") +
  geom_point() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  labs(y = "Coefs of cross term", x = "terms",
       title = "make * Year past")

さて、年数の影響を排除した後の走行距離と価格の交差項を見ていきます。これは少し様子が変わりましたね。三菱やマツダは依然コスパが良いですが、リンカーンやキアも出てきました。ただ、影響の規模は半分以下になっています。

temp_model |> 
  filter(str_detect(term, "mileage")) |> 
  mutate(term = str_remove_all(term, "make|mileage")) |> 
  ggplot() +
  aes(y = estimate, x = reorder(x = term, X = estimate)) +
  facet_wrap(effect ~ ., scales = "free") +
  geom_point() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  labs(y = "Coefs of cross term", x = "terms",
       title = "make * mileage")

いまさらではありますが、念の為相関係数も見ておきます。やはり高いですね。

target_dat |> summarise(correlation = cor(year_past, mileage))
## # A tibble: 1 × 1
##   correlation
##         <dbl>
## 1       0.721

結論

係数の大きさを考えると、日本車の価格が下がりにくいのは走行距離というより、時間が経っても価格が下がりにくいことの影響が大きそうです。ただ、最初のグラフに戻ってみると、ある水準以下にはさがらなさそう=関係は完全な線形ではなさそうなので、機械学習モデルを使ってより精度の高い予測をしてみたいですね。

で、何を買うべきかですが、そうはいっても日本人なので日本車が良いです。走行距離と年式の両方で上位に入ってくる三菱かマツダが候補になりますが、三菱はリコールとかいろいろあったしあまり好きな車がないので、マツダにすべきですね。筆者はもともとマツダファンなので、結論ありきっぽくなっちゃいました。

RでTorch

だいぶ久しぶりの記事です。前回記事から、MBA合格・渡米・入学・就職といろいろありました。

RでもNerural Networkがしたい!

論文でNeural Networkでモデルの推定をするにあたり、Torchを使いたい。しかし、PyTorchは(当然ながら)Pythonしか対応していません。授業の課題で書いたコードも全部Python。授業のコードを見ると、Pythonだと変数が多くてコードが読みにくいんですよね。

しかしなんと、最近はRでもtorchが使えるらしいじゃないですか

英語では結構文献があります。パッケージのオフィシャルサイトとCRANへのリンクを貼っておきます。 CRAN Official Site Cheat Sheet Turotrial

日本語解説がほとんどない!

みなさんご存知TJOさんは記事を書いてくれていますが、日本語の情報がほとんどないです。 https://tjo.hatenablog.com/entry/2020/10/06/090000

とにかく動かしてみるだけならこれで十分ですが、後述の通り、カスタマイズしようとするとつまづきます。こういうのを扱う人はみんなpython派なんですかねえ。。。書けるけども。。。

あと、Torchっていろいろな書き方があるんですね。よくわからない状態でいろいろググってるとさらに混乱します。

本記事では、基本的な使い方の説明は上記記事に譲って、より実践的な使い方について書きます。

実際に使ってみるときに困ること:損失関数のカスタマイズ

今回一番困ったのは損失関数の扱いです。このパッケージは基本的に必要なものをすべて用意してくれているので、L1損失関数だのRMSEだの、自分で書く必要はありません。しかし、全部用意してくれているがゆえに、自分で損失関数を定義したい場合は少々やっかいです。

僕が参照した論文はNeural Networkの損失関数をカスタマイズしていて、テンプレでは実装できませんでした。僕はこの著者のBryan Kellyという先生にこのあたりの高度な統計学を学んだんですが、この人は結構無邪気に非標準なことを要求してくる傾向があるっぽいんですよね。具体的に言うと、この論文ではRMSEにL1 penalization項を追加しろと言ってるんですが、そんなのパッケージでは言及されていないですし。 https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3159577

僕が使った損失関数の具体的な定義は以下の通りです。単純にL1 Penalty Termを追加しただけなので数学的にはとても簡単に見えますが、torchにはこの損失関数が定義されていないので、自分で作る必要があります。

loss function

カスタマイズ方法

結論だけ書くと、nn_n1_lossに倣ってnn_lossを継承した新しいオブジェクトを定義すれば良いです。この関数にテンソルを入れるとMSEとL1をあわせた損失を計算してくれます。 https://github.com/mlverse/torch/blob/main/R/nn-loss.R#L69-#L78

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(torch)
nn_l1_loss_custom <- nn_module(
# https://cran.r-project.org/web/packages/torch/torch.pdf
# https://github.com/mlverse/torch/blob/main/R/nn-loss.R
  "nn_l1_loss_custom",
  inherit = nn_l1_loss,
  initialize = function(reduction = "sum", lambda = 2){
    self$reduction <- reduction
    self$lambda <- lambda
  },
  forward = function(input, target) {
    nnf_mse_loss(input, target) + 
      self$lambda * nnf_l1_loss(input, target, reduction = self$reduction)
  }
)

具体的に実装してみましょう。まずはデータセットを用意します。みんなだいすきパーマーペンギン。

library(palmerpenguins)
dat <- palmerpenguins::penguins |> drop_na()
(dat_num <- dat |> select_if(is.numeric) |> select(-year))
## # A tibble: 333 × 4
##    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
##             <dbl>         <dbl>             <int>       <int>
##  1           39.1          18.7               181        3750
##  2           39.5          17.4               186        3800
##  3           40.3          18                 195        3250
##  4           36.7          19.3               193        3450
##  5           39.3          20.6               190        3650
##  6           38.9          17.8               181        3625
##  7           39.2          19.6               195        4675
##  8           41.1          17.6               182        3200
##  9           38.6          21.2               191        3800
## 10           34.6          21.1               198        4400
## # ℹ 323 more rows

データ分割などの細かいことはおいておき、とにかくtorchにデータを渡します。体重を予測してみましょう。シンプルに3層のみで、面倒なのでEarly stoppingも実装しません。パラメータについてはコードを見ればわかるので省略します。

こちらが一番シンプルっぽいので、こちらを参考にさくっと学習させます。 https://anderfernandez.com/en/blog/how-to-create-neural-networks-with-torch-in-r/

y <- dat_num |> 
    select(body_mass_g) |>
    as.matrix() |> 
    torch_tensor(dtype = torch_float())

x <- dat_num |> 
    select(-body_mass_g) |> 
    as.matrix() |> 
    torch_tensor(dtype = torch_float())

simple_model <- nn_sequential(
    nn_linear(x$size()[2], 32),
    nn_relu(),
    nn_linear(32, 16),
    nn_relu(),
    nn_linear(16, 1)
  )
  
epochs <- 10^4
lambda <- 3
learning_rate <- 0.1

loss_func <- nn_l1_loss_custom(reduction = "sum", lambda = lambda)
opt <- optim_adam(simple_model$parameters, lr = learning_rate)

for (i in 1:epochs) {
  opt$zero_grad()
  y_pred <- simple_model(x)
  loss <- loss_func(y_pred, y)
  loss$backward()
  opt$step()
  
  if (i %% 2000 == 0) {
    cat(" Epoch:", i, "Loss: ", loss$item(), "\n")
  }
}
##  Epoch: 2000 Loss:  568910.9 
##  Epoch: 4000 Loss:  473118.9 
##  Epoch: 6000 Loss:  473842.8 
##  Epoch: 8000 Loss:  465308.3 
##  Epoch: 10000 Loss:  503881.8
cat("learning, done!")
## learning, done!

torchの損失関数による損失計算の結果

y_pred <- simple_model(x) 
loss_func(y, y_pred)$item()
## [1] 513573.4

を手計算で再現

dat_num |> 
    bind_cols(simple_model(x) |> 
                  as_array() |> 
                  as.data.frame() |> 
                  tibble()) |> 
    rename(y_pred = V1) |> 
    summarise(loss_value = mean((y_pred - body_mass_g)^2) + (lambda * sum(abs(y_pred - body_mass_g))))
## # A tibble: 1 × 1
##   loss_value
##        <dbl>
## 1    513573.

完全一致!損失関数を自作できました!

Quantpedia記事翻訳:Man vs. Machine: Stock Analysis

基本情報

Authors

Sean Cao, Wei Jiang, Junbo L. Wang and Baozhong Yang

Link

papers.ssrn.com quantpedia.com

まえがき

近年、機械学習ベースの取引戦術や市場分析が増えている。しかし、機械が我々に取って代わることは可能なのだろうか?アルゴリズムがbig dataを解釈するより高い能力を持っていることは疑いようもない。しかしながら、市場の常ではあるが、全てが合理的なわけではないのだ。Cao et al.(2021)は、株式市場におけるこのテーマを深く分析した。目標価格と利益予想は重要なポイントで、トレーダーや投資家が活用する数字である。最新の研究では、目標価格の予想能力に関して人間の能力と機械の能力を比較した。全体として、AI-based analysisは人間による分析を凌駕したが、ことはそう単純ではない。AIは多くのデータから学習することができるが、人間がそう簡単に取って代わられることはないだろう。人間の特殊性が価値を持つ場面が多くある。例えば、流動性が低く小さな会社や、アセットライトなビジネスモデルを持つ会社だ。さらに言えば、AIと人間は競合しているというよりは補完しあっている。

要旨・内容

  • 企業財務情報、定性開示情報、マクロ経済情報を学習したAI-analystは、価格予測に関してhuman-analystを凌駕し、超過リターンを生み出した。複雑な企業の分析において、透明性が高く量の多い情報が多次元に渡り手に入る場合はAI-analystの優位性が大きかった。
  • 重要な情報の解釈に制度的な知識(無形資産の内容等)が必要な場合はhuman-analystの方が優れていた。
  • human-analystがalternative dataや企業内AIへのアクセスを手に入れるにつれ、AIの優位性は時間と共に減衰していった。
  • AI-analystとHuman-analystを合わせた分析はAI-analystを大きく凌駕し、機械と人間の融合の可能性を見せた。
    f:id:maxonblog:20210808172104p:plain
    Figre-1

f:id:maxonblog:20210808172129p:plain
Figure-2

f:id:maxonblog:20210808172149p:plain
Figure-3

特筆すべきポイント

  • 使用したAIモデルは次の通り: Elastic-Net, Support Vector Machines, Random Forest, Gradient Boosting, and Long Short-Term Memory Neural Networks
  • Man + Machineのpredictionについては、合成方法の記載なし。