市場全体の動きとの連動性を表すベータ値$\beta$は，個別銘柄のリスクを検討するときよく参照される指標です．ベータの推計にあたっては，CAPM(Capital Asset Pricing Model)のフレームワークをもとに，時刻$t$における銘柄$i$のリターン$y_{it}$が以下のように与えられるものと仮定します．
y_{it} = \alpha_i + \beta_i (x_{t} - r_f) + \epsilon_{it}
ただし$x_t$はマーケットポートフォリオ(典型的にはTOPIX)のリターン，$r_f$はリスクフリーレート(長期国債利回りなど)，$\epsilon_{it}$は適当な確率変数です．素朴に考えると$y_{it}$を$x_{t}$に回帰させることで$\beta_i$を予測できますが，ベイズ修正を用いることでその精度を高めることができます．
本記事では，まずVasicekによるモデルを示し，その特別な例として(Merrill Lynchによるとも言われる)以下のベイズ修正を導出します．
$$
\mathbb{E} [\beta_i \mid {\hat \beta}_i] = \dfrac{1}{3} + \dfrac{2}{3} {\hat \beta}_i
$$
さらにより一般的な形で再度定式化し，ハイパーパラメータをエビデンス最大化により求めます．
今回の分析で使ったスクリプトは，以下に掲載しています．
https://github.com/yotapoon/Adjusted_Beta
Vasicekモデル
企業$i = 1, \cdots, M$の真のベータ$\beta_i$は，平均$\mu_{\beta}$，分散$\sigma_{\beta}^2$の事前分布から確率的に生成されるものとします．さらに過去のリターンから算出されるヒストリカルベータ$\hat{\beta}_i$は，真のベータ$\beta_i$にノイズ$\sigma_i$を伴っているものと仮定します．ここでは，具体的な確率分布として以下を想定します．
p(\beta_i) = \dfrac{1}{\sqrt{2 \pi \sigma_\beta^2}} \exp \left[ - \dfrac{(\beta_i - \mu_\beta)^2}{2 \sigma_\beta^2} \right] \\
p({\hat \beta}_i \mid \beta_i) = \dfrac{1}{\sqrt{2 \pi \sigma_i^2}} \exp \left[ - \dfrac{({\hat \beta}_i - \beta_i)^2}{2 \sigma_i^2} \right]
ベイズの枠組みに従って計算すると，以下の事後分布が求まります(Vasicekモデル)．
p(\beta_i \mid {\hat \beta}_i) \propto \exp \left[ -\dfrac{1}{2} \left(\dfrac{1}{\sigma_\beta^2} +
\dfrac{1}{\sigma_i^2} \right) \left(\beta_i - \dfrac{\sigma_i^2 \mu_\beta + \sigma_\beta^2 {\hat \beta}_i}{\sigma_i^2 + \sigma_\beta^2}\right) \right]
以上から，$\beta_i$の推計値として以下が得られます．
\mathbb{E} [\beta_i \mid {\hat \beta}_i] = \dfrac{\sigma_i^2}{\sigma_i^2 + \sigma_\beta^2} \mu_\beta + \dfrac{\sigma_\beta^2}{\sigma_i^2 + \sigma_\beta^2} {\hat \beta}_i
特に$\mu_{\beta} = 1$かつ$\sigma_i^2 : \sigma_\beta^2 = 1 : 2$となることを仮定すると，よく用いられるベータの修正値
\mathbb{E} [\beta_i \mid {\hat \beta}_i] = \dfrac{1}{3} + \dfrac{2}{3} {\hat \beta}_i
が得られます．特に$\mu_{\beta} = 1$は，マーケットポートフォリオのウエイトによるベータの加重平均が$1$となることから決め打ちされているものと考えられます．
一般化されたBayes修正
Bayes修正を，より一般化された形で整理します．企業$i$のリターン
${\bf y} _ i = (y _ {i 1} , \cdots, y _ {i N})$を，定数項およびマーケットポートフォリオのリターンを含む$n$個の以下のファクター$X$で回帰することが目的です．
X =
\begin{pmatrix}
x_{11} & \dots & x_{1j} & \dots & x_{1n} \\
\vdots & & \vdots & & \vdots \\
x_{t1} & \dots & x_{tj} & \dots & x_{tn} \\
\vdots & & \vdots & & \vdots \\
x_{N1} & \dots & x_{Nj} & \dots & x_{Nn}
\end{pmatrix}
事後分布によるベータの推計
企業$i$の(一般化された)ベータを${\boldsymbol \beta}_i$とし，以下が成り立つことを仮定します．
p({\bf y}_i \mid {\boldsymbol \beta}_i, X, \sigma_i^2) = \dfrac{1}{(2 \pi)^{N/2} \sigma_i^N} \exp \left[ -\dfrac{1}{2 \sigma_i^2} ({\bf y}_i - X {\boldsymbol \beta}_i )\right]
さらにハイパーパラメータ$\Lambda$，${\boldsymbol \mu}$のもとで${\boldsymbol \beta}_i$の事前分布が
p({\boldsymbol \beta}_i \mid {\boldsymbol \mu}, \Lambda) = \dfrac{|\Lambda|^{1/2}}{(2 \pi)^{n/2}}
\exp \left[ -\dfrac{1}{2} ({\boldsymbol \beta}_i - {\boldsymbol \mu})^t \Lambda ({\boldsymbol \beta}_i - {\boldsymbol \mu}) \right]
であるものとすると，事後分布は
p({\boldsymbol \beta}_i \mid {\boldsymbol \mu}, \Lambda, {\bf y}_i ) \propto
\exp \left[ -\dfrac{1}{2} ({\boldsymbol \beta}_i - A_i^{-1} {\bf z}_i)^t A_i ({\boldsymbol \beta}_i - A_i^{-1} {\bf z}_i) \right]
のように求められます．ただし，以下のパラメータに注意します．
\alpha_i = \dfrac{1}{\sigma_i^2} \\
A_i = \Lambda + \alpha_i X^T X \\
{\bf z} _ i = \Lambda {\boldsymbol \mu} + \alpha_i X^t {\bf y} _ i
以上からベータの推計には事後分布の期待値
\hat{\boldsymbol \beta}_i = A_i^{-1} {\bf z}_i = ( \Lambda + \alpha_i X^T X)^{-1} ( \Lambda {\boldsymbol \mu} + \alpha_i X^t {\bf y} _ i)
を用いればよいとわかります．ただし，この期待値にはハイパーパラメータ$\Lambda$，${\boldsymbol \mu}$，${\boldsymbol \alpha} = (\alpha_1, \cdots, \alpha_M) = (\sigma_1^{-2}, \cdots, \sigma_M^{-2})$が含まれているため，エビデンス$p(Y \mid \Lambda, {\boldsymbol \mu}, {\boldsymbol \alpha})$の最大化を条件としてハイパーパラメータを適切に設定します．
エビデンス近似によるハイパーパラメータ推計
今回のモデルはベイズ線形回帰(を若干拡張したもの)であるため，以下のエビデンス$ p(Y \mid {\boldsymbol \mu}, \Lambda, {\boldsymbol \alpha})$を解析的に計算できます：
p(Y \mid {\boldsymbol \mu}, \Lambda, {\boldsymbol \alpha}) = \prod_i p({\bf y} _ i \mid {\boldsymbol \mu}, \Lambda, {\boldsymbol \alpha})
= \prod_i \int p({\bf y} _ i \mid {\boldsymbol \beta} _ i, X, \alpha_i) p({\boldsymbol \beta} _ i \mid {\boldsymbol \mu}, \Lambda) d {\boldsymbol \beta} _ i
積分を実行し，定数項を無視することで，ハイパーパラメータ最適化の目的関数は以下の対数エビデンスとなります．
p(Y \mid {\boldsymbol \mu}, \Lambda, {\boldsymbol \alpha})
= \dfrac{M}{2} \log |\Lambda| + \dfrac{N}{2} \sum_{i=1}^{M} \log \alpha_i - \dfrac{1}{2} \sum_{i=1}^{M} \log | A_i |
+ \dfrac{1}{2} \sum_{i=1}^M {\bf z} _ i ^ t A_i^{-1} {\bf z} _ i
- \dfrac{M}{2} {\boldsymbol \mu}^t \Lambda {\boldsymbol \mu} - \dfrac{1}{2} \sum_{i=1}^M \alpha_i {\bf y}_i^t {\bf y}_i + {\rm const.}
対数エビデンスを最大化するハイパーパラメータ${\boldsymbol \mu}, \Lambda, {\boldsymbol \alpha}$を求めるためには，上式を(数値的に)微分する必要があります．
以上をまとめると，次のような手続きとなります．
1. ハイーパーパラメータの初期値の設定
何も工夫をしない場合，次のように設定するという方法があります：${\boldsymbol \mu} = {\bf 0},~ \Lambda = 1_n,~ {\boldsymbol \alpha} = 1$
ナイーブな線形回帰の結果をもとにすると，以下のように設定する方がより適切でしょう．
${\boldsymbol \mu}$：全銘柄について，回帰係数$\hat{\boldsymbol \beta}_i$を平均したベクトルとする
$\Lambda$：全銘柄について，回帰係数$\hat{\boldsymbol \beta}_i$の共分散行列の逆行列とする
${\boldsymbol \alpha}$：線形回帰による推定値の誤差分散の逆数とする
2. 数値微分によりエビデンス最大化
エビデンス$ p(Y \mid {\boldsymbol \mu}, \Lambda, {\boldsymbol \alpha})$を数値微分することで，適切な${\boldsymbol \mu}, \Lambda, {\boldsymbol \alpha}$を求めます．
3. 数値微分によりエビデンス最大化
ベータの事後分布$p({\boldsymbol \beta}_i \mid {\boldsymbol \mu}, \Lambda, {\boldsymbol \alpha})$の期待値$\hat{\boldsymbol \beta}_i$が，ベータの推計値となります：
\hat{\boldsymbol \beta}_i = ( \Lambda + \alpha_i X^T X)^{-1} ( \Lambda {\boldsymbol \mu} + \alpha_i X^t {\bf y} _ i)
なお，事後分布に基づくベータの分散は$A_i^{-1} = (\Lambda + \alpha_i X^T X)^{-1}$の対角成分として求められます．
Pythonによる実装
一般化されたベイズ修正について，Pythonでは以下のように実装することができます．
class BayesLinearRegressionBeta:
def __init__(self, mu, Lambda, alpha): # Yは(N, n_stocks)
self.mu, self.Lambda, self.alpha = mu, Lambda, alpha
# 入力サイズで結果が不安定にならないように，全体を値を(N*n_stocks)でスケールしている
def evidence(self, mu, Lambda, alpha):
evidence_value = (0.5/self.N) * np.log(np.linalg.det(Lambda))
evidence_value += -(0.5/self.N) * mu.dot(Lambda).dot(mu)
evidence_value += 0.5*np.log(alpha).mean()
A = [[] for _ in range(self.n_stocks)]
z_vec = [[] for _ in range(self.n_stocks)]
for i in range(self.n_stocks):
y_i = self.Y[:, i]
A[i] = Lambda + alpha[i] * self.X.T.dot(self.X)
z_vec[i] = Lambda.dot(mu) + alpha[i] * self.X.T.dot(y_i)
evidence_value += -0.5*np.log(np.linalg.det(A[i])) / (self.n_stocks*self.N)
evidence_value += 0.5*z_vec[i].dot(np.linalg.inv(A[i])).dot(z_vec[i]) / (self.n_stocks*self.N)
evidence_value += -0.5*alpha[i]*y_i.dot(y_i) / (self.n_stocks*self.N)
return evidence_value
def fit(self, Y, X, eta = 0.01, epsilon = 0.001, T = 100):
self.X, self.Y = X.copy(), Y.copy()
self.N, self.n_factors, self.n_stocks = len(X), len(X[0]), len(Y[0])
self.eta = eta
self.evidence_list = [[] for _ in range(T)]
for t in range(T):
# エビデンスを計算→格納
self.evidence_t = self.evidence(self.mu, self.Lambda, self.alpha)
self.evidence_list[t] = self.evidence_t
print(t, self.evidence_t, end = "\n")
if math.isnan(self.evidence_t):
print("!!!ERROR!!!")
break
self.update_mu(epsilon)
self.update_alpha(epsilon)
self.update_Lambda(epsilon)
def estimate(self):
self.var_beta = np.zeros((self.n_factors, self.n_stocks))
self.W = np.zeros((self.n_factors, self.n_stocks))
for i in range(self.n_stocks):
y_i = self.Y[:, i]
A_i = self.Lambda + self.alpha[i] * self.X.T.dot(self.X)
z_vec_i = self.Lambda.dot(self.mu) + self.alpha[i] * self.X.T.dot(y_i)
self.W[:, i] = np.linalg.inv(A_i).dot(z_vec_i)
self.var_beta[:, i] = np.linalg.inv(A_i).diagonal()
def update_mu(self, epsilon):
grad = np.zeros(self.n_factors)
for i in range(self.n_factors):
mu_new = self.mu.copy()
mu_new[i] += epsilon
evidence_t_new = self.evidence(mu_new, self.Lambda, self.alpha)
grad[i] = (evidence_t_new - self.evidence_t) / epsilon
self.mu += self.eta * grad
def update_alpha(self, epsilon):
grad = np.zeros(self.n_stocks)
for i in range(self.n_factors):
alpha_new = self.alpha.copy()
alpha_new[i] += epsilon
evidence_t_new = self.evidence(self.mu, self.Lambda, alpha_new)
grad[i] = (evidence_t_new - self.evidence_t) / epsilon
self.alpha += self.eta * grad
def update_Lambda(self, epsilon):
grad = np.zeros((self.n_factors, self.n_factors))
for i in range(self.n_factors):
for j in range(self.n_factors):
Lambda_new = self.Lambda.copy()
Lambda_new[i, j] += epsilon
Lambda_new[j, i] += epsilon
evidence_t_new = self.evidence(self.mu, Lambda_new, self.alpha)
grad[i, j] = (evidence_t_new - self.evidence_t) / epsilon
self.Lambda += self.eta * grad
class LinearRegressionBeta:
def __init__(self, lambda_ridge = 0.0):
self.lambda_ridge = lambda_ridge
def fit(self, Y, X): # Yは(N, n_stocks)
N, n_factors = len(X), len(X[0])
self.Y = Y.copy()
self.X = X.copy()
Lambda_ridge = self.lambda_ridge*np.eye(n_factors)
self.W = np.linalg.inv(self.X.T.dot(self.X) + Lambda_ridge).dot(self.X.T).dot(self.Y)
self.Y_pred = self.X.dot(self.W)
self.sigma2 = ((self.Y_pred - self.Y) ** 2.0).sum(axis = 0) / (N - n_factors) # 標準誤差(の2乗)
self.cov_beta = np.linalg.inv(self.X.T.dot(self.X) + Lambda_ridge) * self.sigma2.reshape(-1, 1, 1)
self.var_beta = (self.cov_beta * np.eye(n_factors)).sum(axis = 1).T # Wと同じ形でbetaの分散
model_LR_beta = LinearRegressionBeta()
model_LR_beta.fit(Y, X)
mu_initial = model_LR_beta.W.mean(axis = 1)
Lambda_initial = np.linalg.inv(np.cov(model_LR_beta.W))
alpha_initial = 1.0 / model_LR_beta.sigma2
model_bayes = BayesLinearRegressionBeta(mu_initial, Lambda_initial, alpha_initial)
model_bayes.fit(Y, X) # ハイパーパラメータ最適化
model_bayes.estimate() # ベータの推計
数値シミュレーション
素朴な線形回帰やVasicekモデルと比較して，エビデンス最大化に基づくベイズ修正によってベータの推計精度がどの程度向上しているかを確認します．
まず，対象となる銘柄の期待リターンとリスクを適当に設定し，マーケットポートフォリオを算出します．リスク・リターン平面上では以下のようになります．
今回は期待リターンとリスク及び相関係数が既知であるため，ベータの値は厳密に以下のように求められます：
\beta_i = \dfrac{\sigma_{iM}}{\sigma_M^2}
ただし$\sigma_{iM}$は銘柄$i$とマーケットポートフォリオのリターンの共分散，$\sigma_M$はマーケットポートフォリオのリスクです．そのため，真のベータと3つの推計方法(線形回帰，Vasicekモデル，ベイズ修正)による値を比較できます．以下の左図は線形回帰・Vasicekモデルによる推計結果で，真の値$y = x$とは乖離が生じていることが確認できます．一方右図はベイズ修正によるモデルであり，推計精度の向上を確認できます．
真の値との乖離を平均二乗誤差(MSE)で見ると，以下のようなグラフになっています．
まずナイーブな線形回帰に比べ，ベイズ修正によるモデルの誤差が小さくなっていることが確認できます．両者の相違点は「他銘柄のリターンの情報を用いているか」というところにあります．マーケットポートフォリオが他の銘柄の組み合わせからなる以上，この差が精度向上に寄与したものと考えられます．またVasicekモデルについても若干の精度向上が見られますが，期待リターンとリスクの設定によってはむしろ誤差が大きくなる場合がありました．これは，Vasicekモデルの仮定(ベータの期待値が1となることなど)が非常に強く，設定がこの仮定と整合していない場合には推計精度が悪くなるためと考えられます．
実データによる分析
TOPIX Core30とTOPIX Large70を構成する約100社を対象に，同様に推計結果を比較しました．ただし，企業の真のベータ値を推計することはできないため，2017年11月末～2022年11月末の60カ月ベータを真の値として，2012年11月末～2017年11月末の推計値と比較するという方法を取ります(この方法については後で議論)．
数値シミュレーションと同様に，手法ごとに推計値の平均二乗誤差(MSE)を見たものが以下のグラフです．まず，数値シミュレーションと同様，ナイーブな線形回帰に比べるとBayes修正の精度が向上しています．一方で，VasicekモデルはBayes修正よりも精度が高くなっています．理由としては，
- TOPIX構成銘柄のうち時価総額が大きいものを取ってきているため，ベータが1に近い銘柄が多く，Vasicekモデルの仮定と整合していた
- 比較している2つの時点(2022年11月末と2017年11月末)では各企業の状況が大きく変化し，そもそもベータが変化している
可能性があります．
まとめ
本記事では，まずVasicekモデルをもとにMerrill Lynchの調整ベータを導出し，次により一般的なベイズ修正を議論しました．エビデンス最大化によるハイパーパラメータ最適化の方法を示し，数値シミュレーションと実データをもとに線形回帰，Vasicekモデル，Bayes修正の推計方法を比較しました．少なくとも単純な線形回帰に比べると，Bayes修正の方法はより良い結果を与える可能性がありそうです．Vasicekモデルは，その仮定に対する整合性が重要となると考えられます．ただし実データを用いたベータの比較方法については，検討の余地はありそうです．なお，関連する論文を以下に掲載しています．
https://repository.musashi.ac.jp/dspace/bitstream/11149/1817/1/musashi_2016no64_01_005.pdf
ディスクレイマー
本記事は筆者の個人的なメモであり，情報の正確性，信頼性，完全性は保証されず，本情報に基づいて被ったいかなる損害についても，一切の責任を負いかねます．