おめさんの備忘録

特定の金融商品/銘柄の売買の勧誘/推奨等を目的とするものではないです。市場全般の推奨や証券市場等の動向の上昇/下落を示唆するものでもないです。掲載の金融商品/資金運用等による投資によって生じたいかなる損失について、一切の責任を負いません。

コロンビア大 周辺のジム巡り

Dodge Fitness Center

大学のジム。キャンパス内にある。学生証があれば追加の会費なし利用できる。

下は年末のガラガラのときに撮った。普段はかなり混んでる...。
深夜でもキャンパス内とその周辺は安全なので、閉館ギリギリの時間帯か、朝早くいくとよい。

ラック・ダンベル・マシン・カーディオ系は一通り揃ってるけど、常に乱雑に散らかっている。

そしてなぜか色々とぶっ壊れている。

たまにイベントがあって閉館する。
今月(2024年3月)はアイビーリーグ関連の運動会で4日間しまってた。


② New York Sport Club

ハーレム方面にある民間のジム。

普段どれほど混んでいるかはわからないけど、2023年末に1 Day Free Passで行ってみたときはかなり空いてた。
下はその時の写真。

Dodgeと同程度かそれ以上に器具が充実してる。

1 Day Free Passを利用すると、その後1ヶ月間くらいまで営業電話が毎日かかってくる。

NYSCの直ぐ側にある、普通の民間ジム。普通すぎて写真が↓しかない。

こちらも1 Day Free Passで一日限り利用した。ラック・ダンベルあたりのエリアが混んでた。

④ Planet Fitness 125th St

ここも1 Day Free Pass。
一番整ってる。ただしラック・バーベルがない。
写真を撮ると怒られそうな雰囲気が漂っていたので人が写ってない写真しかない。




ちなみに、こっちの民間のジムは退会に関して面倒な手続きとかポリシーがあることが多い。
短期間だけ利用したい場合は、1ヶ月のみ入会するのではなく、1(~5) Day Free Passで転々と利用するか、もしくはClasspassというアプリの1 month trialでCrunch gymという民間のジムに通うなどして、特定のジムにcommitしないことを勧めたい。

JuliaでSparse行列のCholesky Decomposition

let
    N = 1000000

    a = rand(N)
    S = spdiagm(0 => ones(N+1), -1 => a)
    S = S' * S

    C_1 = cholesky(S)
    L_1 = sparse(C_1.L)

    C_2 = cholesky(S, perm=1:N+1)
    L_2 = sparse(C_2.L)
    
    @assert !(L_1 * L_1' ≈ S)
    @assert L_2 * L_2' ≈ S
end

例:精度行列から多変量正規に従うRVを乱数生成

let
    μ + (sparse(cholesky(Σ_inv, perm=1:D).L)' \ rand(MvNormal(I(D))))
end

Dirichlet Process Mixture ModelのSlice Sampling

これは何

・DPMMのslice samplingによるMCMCが(おそらく)上手くいったよという記録。
・けんきうでつかう。参考文献を記録しておきたい。

  • Casarin, R., Costantini, M., and Osuntuyi, A. (2023), "Bayesian Nonparametric Panel Markov-Switching GARCH Models", Journal of Business and Economic Statistics; https://doi.org/10.1080/07350015.2023.2166049
  • Kalli, M., Griffin, J. E., and Walker, S. G. (2011), "Slice Sampling Mixture Models", Statistics and Computing, 21, 93–105.
  • Liu, J. (2021), "A Bayesian Semiparametric Realized Stochastic Volatility Model", Journal of Risk and Financial Management, 14(12); https://doi.org/10.3390/jrfm14120617
  • Walker, S. G. (2007), "Sampling the Dirichlet Mixture Model with Slices", Communications in Statistics - Simulation and Computation, 36, 45–54.

メリカのデカい食べもの

アイスクリーム。41オンス(※約1.16233キロ)。コレストロール0mg!(?)


Rosa's Pizzaというお店のチキンなんとかピザ、1スライス。5.99ドルだったと思う。皿からはみ出とる。ルームメイトがノリで2スライス頼んでて食べきるのに四苦八苦してた。
www.yelp.com


バナナ(?)。怖い。


Smith & WollenskyというステーキハウスのFilet mignon。500gはありそうだった。美味しかった。
www.yelp.com

Bayesian線形回帰のGibbs Sampler

これは何

Bayesian線形回帰のGibbs sampler.
よくわすれちゃうから


\def\abs#1{\lvert{#1}\rvert}
\def\set#1{\left\{{#1}\right\}}


\def\Bracket#1#2{\left({#1}|{#2}\right)}
\def\bracket#1{\left({#1}\right)}
\def\sqbracket#1{\left[{#1}\right}
]

Probabilistic Model

 \boldsymbol{y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}, ~~~ \boldsymbol{\epsilon} \sim \mathcal{N}(\boldsymbol{0}_N, \boldsymbol{\Sigma}_{\epsilon}); ~~~ i=1,\dots,N.

なお、
 \boldsymbol{y}=(y_1,\dots,y_N)^{\prime},
 \boldsymbol{X}:\text{design matrix},
 \boldsymbol{\beta}=(\beta_0, \beta_1, \dots, \beta_k)^{\prime},
 \boldsymbol{\Sigma}_{\epsilon}=\sigma^2_{\epsilon}\boldsymbol{I}_Nだけど必ずしもそうでなくてよい。

Priorは
 \begin{align*}
        p(\boldsymbol{\beta})&=(2\pi)^{-N/2}\abs{\boldsymbol{V}_{\beta}}^{-1/2}\exp{\set{-\frac{1}{2}(\boldsymbol{\beta}-\boldsymbol{M}_{\beta})^{\prime}\boldsymbol{V}_{\beta}^{-1}(\boldsymbol{\beta}-\boldsymbol{M}_{\beta})}}=_d\mathcal{N}(\boldsymbol{M}_{\beta},\boldsymbol{V}_{\beta}), \\
        p(\sigma^2_{\epsilon})&=\frac{(\beta_{\sigma})^{\alpha_{\sigma}}}{\Gamma(\alpha_{\sigma})}\bracket{\sigma^2_{\epsilon}}^{-(\alpha_{\sigma}+1)}\exp{\set{-\frac{\beta_{\sigma}}{\sigma^2_{\epsilon}}}}=_d\text{InvGamma}(\alpha_{\sigma},\beta_{\sigma}).
    \end{align*}

Likelihoodは、 \boldsymbol{y} \sim \mathcal{N}(\boldsymbol{X\beta},\sigma^2_{\epsilon}\boldsymbol{I}_N)より、

 \begin{align*}
        p\Bracket{\boldsymbol{y}}{\boldsymbol{\beta},\sigma^2_{\epsilon}}
        &= \abs{2\pi\boldsymbol{\Sigma}_{\epsilon}}^{-1/2}\exp{\set{-\frac{1}{2}(\boldsymbol{y}-\boldsymbol{X\beta})^{\prime}\boldsymbol{\Sigma}_{\epsilon}^{-1}(\boldsymbol{y}-\boldsymbol{X\beta})}} \\
        &= (2\pi\sigma^2_{\epsilon})^{-N/2}\exp{\set{-\frac{1}{2\sigma^2_{\epsilon}}(\boldsymbol{y}-\boldsymbol{X\beta})^{\prime}(\boldsymbol{y}-\boldsymbol{X\beta})}}.
    \end{align*}

Full Conditionals

 \begin{align*}
        p\Bracket{\sigma^2_{\epsilon}}{\boldsymbol{y},\boldsymbol{\beta}}
        & \propto p(\boldsymbol{\beta})P(\sigma^2_{\epsilon})p\Bracket{\boldsymbol{y}}{\boldsymbol{\beta},\sigma^2_{\epsilon}} \\
        & \propto (\sigma^2_{\epsilon})^{-(N/2+\alpha_{\sigma}+1)}\exp{\set{-\frac{1}{\sigma^2_{\epsilon}}\bracket{\frac{1}{2}(\boldsymbol{y}-\boldsymbol{X\beta})^{\prime}(\boldsymbol{y}-\boldsymbol{X\beta})+\beta_{\sigma}}}} \\
        &= \text{InvGamma}\bracket{\frac{N}{2}+\alpha_{\sigma},\frac{1}{2}(\boldsymbol{Y}-\boldsymbol{X\beta})^{\prime}(\boldsymbol{Y}-\boldsymbol{X\beta})+\beta_{\sigma}}.
    \end{align*}


\begin{align*}
        p\Bracket{\boldsymbol{\beta}}{\boldsymbol{y},\sigma^2_{\epsilon}}
        &\propto
        \exp{\set{
        -\frac{1}{2}\bracket{
        \boldsymbol{\beta}^{\prime}\boldsymbol{\Sigma}_{\beta}^{-1}\boldsymbol{\beta}
        -2\boldsymbol{\beta}^{\prime}\boldsymbol{\Sigma}_{\beta}^{-1}\boldsymbol{\mu}_{\beta}
        }}}.
    \end{align*}

ここで、
 \boldsymbol{\Sigma}_{\beta}=\bracket{\boldsymbol{V}_{\beta}^{-1}+\boldsymbol{X}^{\prime}\boldsymbol{X}/\sigma^2_{\epsilon}}^{-1}, ~~~ \boldsymbol{\mu}_{\beta}=\boldsymbol{\Sigma}_{\beta}\bracket{\boldsymbol{V}_{\beta}^{-1}\boldsymbol{M}_{\beta} + \boldsymbol{X}^{\prime }\boldsymbol{y}/\sigma^2_{\epsilon}}.

よって、 p\Bracket{\boldsymbol{\beta}}{\boldsymbol{y},\sigma^2_{\epsilon}}\sim\mathcal{N}(\boldsymbol{\mu}_{\beta},\boldsymbol{\Sigma}_{\beta}).

上のFC  p\Bracket{\boldsymbol{\beta}}{\text{他}} は下の特別な場合。
 \boldsymbol{\Sigma}_{\beta}=\bracket{\boldsymbol{V}_{\beta}^{-1}+\boldsymbol{X}^{\prime}\boldsymbol{\Sigma}_{\epsilon}^{-1}\boldsymbol{X}}^{-1}, ~~~ \boldsymbol{\mu}_{\beta}=\boldsymbol{\Sigma}_{\beta}\bracket{\boldsymbol{V}_{\beta}^{-1}\boldsymbol{M}_{\beta} + \boldsymbol{X}^{\prime }\boldsymbol{\Sigma}_{\epsilon}^{-1}\boldsymbol{y}}.

(Exponential familyのcanonical formにおける自然パラメータを彷彿とさせる形。)

IGもIWと同様の関係にある。

やってみよう

学びたてのころにPyで書いたコードをおいてみる。

期初 r=0において、初期値  \boldsymbol{\beta}_{(r)}\in\mathbb{R}^{k+1},(\sigma^2_{\epsilon})_{(r)}\in\mathbb{R}_{++} を設定する。 r=1,\dots, Rで、

  •  \boldsymbol{\beta}_{(r)} \sim p\Bracket{\boldsymbol{\beta}}{\boldsymbol{y},(\sigma^2_{\epsilon})_{(r-1)}}
  •  (\sigma^2_{\epsilon})_{(r)} \sim p\Bracket{\sigma^2_{\epsilon}}{\boldsymbol{y},\boldsymbol{\beta}_{(r)}}

とサンプリングを行う。

ライブラリなど
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from tqdm import tqdm
R = 3_000
BURN_IN = 1_000
N = 500

trueBeta = np.array([[2, 4]]).transpose()
K = len(trueBeta)

trueSigma2 = 2
"""Artificial Data"""
X = [np.ones((N, 1)), np.sort(np.random.uniform(size=(N, 1)), axis=0)]
X += [np.random.uniform(size=(N, 1)) for _i in range(K - 2)]
X = np.hstack(X)

U = np.random.normal(loc=0, scale=trueSigma2, size=(N, 1))

Y = X.dot(trueBeta) + U
"""r = 0"""
Beta_r0 = np.linalg.solve(X.transpose().dot(X), X.transpose().dot(Y)) # ols
sigma2_r0 = ((Y - X.dot(Beta_r0)) ** 2).mean() # mle

Theta = dict()
Theta['Beta'] = np.zeros((BURN_IN + R + 1, K))
Theta['Beta'][0] = Beta_r0.flatten()
Theta['sigma2'] = np.zeros((BURN_IN + R + 1, 1))
Theta['sigma2'][0] = sigma2_r0

実行

for r in tqdm(range(1, R + BURN_IN + 1)):
    Sigma_beta = np.linalg.solve(np.linalg.inv(V_Beta_prior()) + X.transpose().dot(X) / Theta['sigma2'][r-1], np.eye(K))
    Mu_beta = Sigma_beta.dot(np.linalg.inv(V_Beta_prior()).dot(M_Beta_prior()) + X.transpose().dot(Y) / Theta['sigma2'][r-1])
    
    C = np.linalg.cholesky(Sigma_beta)
    Beta_r = Mu_beta + C.dot(np.random.randn(K, 1))
    Theta['Beta'][r] = Beta_r.flatten()
    
    e = Y - X.dot(Beta_r)
    sigma2_r = 1 / np.random.gamma(N/2 + alpha_sigma2_prior, 1 / (e.transpose().dot(e) / 2 + beta_sigma2_prior))
    Theta['sigma2'][r] = sigma2_r
col = 'Beta'
ax = pd.DataFrame(Theta[col]).add_prefix(f'{col}_').plot(figsize=(14,6), alpha=0.6)
ax.vlines(BURN_IN, Theta[col].min(), Theta[col].max())
ax.hlines(trueBeta, 0, BURN_IN + R, linestyles='--', label='True')
ax.legend()


col = 'Beta'
ax = pd.DataFrame(Theta[col]).add_prefix(f'{col}_').plot.hist(bins=R // 10, figsize=(14,6), alpha=0.6)
ax.vlines(Theta[col][BURN_IN:].mean(axis=0), 0, 60)
ax.vlines(trueBeta, 0, 60, linestyles='--', label='True')
ax.legend()


Theta['Beta'][BURN_IN:][::100].mean(axis=0)
# array([2.03793947, 4.02184397])
fig = plt.figure(figsize=(9, 6))
ax = fig.add_subplot(
    xlabel='X',
    ylabel='Y',
    # title='Reg Lines over Draws',
)

ax.scatter(X[:, 1], Y, label='Obs')#, c=colors, cmap=plt.get_cmap("viridis"), lw=0)

_X = [np.ones((N, 1)), np.sort(np.random.uniform(size=(N, 1)), axis=0)]
_X += [np.linspace(0, 1, N).reshape(-1,1) for _i in range(K - 2)]
_X = np.hstack(_X)

Betas = Theta['Beta'][BURN_IN:][::30]#.transpose()
for _r in range(len(Betas)):
    ax.plot(_X[:,1], _X.dot(Betas[[_r]].transpose()), linestyle='solid', alpha=0.4, color='r', label=_r*'_' + 'Reg')

ax.plot(_X[:,1], _X.dot(trueBeta), color='black', label="'True'")
ax.legend()

Chan, J., Koop, G., Poirier, D., & Tobias, J. (2019). Bayesian Econometric Methods (2nd ed., Econometric Exercises). Cambridge: Cambridge University Press.

レジーム転換モデル+ベイジアンニューラルネットで景気後退確率を予測する方法

これは何

↓を読んでいたら面白そうなことを思いついたので勢いで記す。

Kim, Y. and Kang, K. (2022) "Bayesian Inference of Multivariate Regression Models with Endogenous Markov Regime-Switching Parameters", Journal of Financial Econometrics, 20(3), 391–436.
https://academic.oup.com/jfec/article/20/3/391/5909218

内生的レジーム転換モデル

観測可能な確率ベクトル  (\textbf{y}_t, \textbf{x}_t) について次のDGPを仮定。

 \underset{j \times 1}{\textbf{y}_t} = \underset{j \times 1}{\textbf{m}(s_t)} + \underset{j \times k}{\textbf{B}(s_t)} ~ \underset{k \times 1}{\textbf{x}_t} + \underset{j \times j}{\textbf{L}(s_t)} ~ \underset{j \times 1}{\textbf{e}_t}

ここで、
 \textbf{L}(s_t)\textbf{L}(s_t)':\text{正定値}, \quad s_t \in S:=\{1,...,N\}, \quad t \in \mathbb{Z}.

追加的に、次の仮定が置くことで、レジーム転換を endogenous とする。

 
\underset{(N-1 + j) \times 1}{\begin{bmatrix} \textbf{z}_t \\ \textbf{e}_t \end{bmatrix}}
\stackrel{\text{iid}}{\sim}
\mathcal{N} \left(
\begin{bmatrix} \textbf{z} \\ \textbf{e} \end{bmatrix} \mid
\begin{bmatrix} \textbf{g}_{s_{t-1}} \\ \textbf{0}_j \end{bmatrix},
\begin{bmatrix} \textbf{S}_z & \textbf{R} \\ \textbf{R}' & \textbf{I}_j \end{bmatrix}
\right),
\qquad
s_t = f_{\textbf{z}}(\textbf{z}_t)

ここで、

  •  \textbf{S}_z \in \mathbb{R}^{(N-1) \times (N-1)} は正定値;
  •  \textbf{R} = \textbf{0}_{(n-1) \times j} \Leftrightarrow (\textbf{z}_t \bot \textbf{e}_t) (i.e. レジーム転換が exogenous).

適当に事前分布を設定したのちに、下のパラメータを推定したい。

  •  (\textbf{m}(s), \textbf{B}(s), \text{vech}(\textbf{L}(s)))_{s \in S};
  •  (s_t)_{t\in\mathbb{Z}};
  •  (\textbf{S}_z, \textbf{R})

拡張

元論文では、 \textbf{g}_{s_{t-1}} \in \mathbb{R}^{N-1} は定数であった。
多分、こうすることで、 (s_t)_{t \in \mathbb{Z}} の Markovian property を明確に保持したかったんだと推測。

一方、状態の予測(フォアキャスト)に生かしたいのであれば、 \textbf{z}_t の marginal な mean  \mathbb{E}\textbf{z} は予測に使えそうな time-varying な量でもあってほしい。そこで、

 f_s = f_1 \circ \dots \circ f_{L-1} \circ f_L  : \mathbb{R}^{N-1} \rightarrow \mathbb{R}^{N-1}; \quad \mathbb{E}\textbf{z} = f_s(\textbf{x}_{t-1})

という (Bayesian) feedforward neural network を組み込む。FFNNでなければならない理由はない。とにかく、↑は、ある状態への遷移確率を NN でモデル化したことに相当する。特に、それは  \textbf{x}_{t-1} に依る。

従って、Filardo (1994) と異なる点は、

  • 時変的遷移確率を多項ロジット的な単純なモデルでパラメータ化せず、agnostic に Bayesian Neural Net を取り入れた点;
  • Kim and Kang (2022) の枠組みのもとで、endogeneity を含む点;
  • (かつ、 1次元の潜在変数による endogeneity でなく、 N-1次元の潜在変数による endogeneity においてBayesian Neural Net を取り入れた点)

など。

簡単な設定で試す

データ

McCracken and Ng (2015) の FRED-MD の

の月次データ(2000-10 ~ 2023-01)をフィーリングで選択。追加でFF先物も入れている。

観測方程式(レジームスイッチング・2変量VAR)


\textbf{y}_t = \textbf{m}(s_t) + \textbf{B}(s_t) ~ \textbf{y}_{t-1} + \textbf{L}(s_t) ~ \textbf{e}_t;
\qquad
\textbf{y}_t = \begin{bmatrix} \text{米 鉱工業指数}_t \\ \text{米 失業率}_t \end{bmatrix}.

NNのアーキテクチャ

計算リソース的な観点から、むやみにパラメータ数を増やせないため、
(1) 入力次元は5次元( \textbf{x}_t \in \mathbb{R}^5 \text{ for } t=1,...,T-1 だから);
(2) 隠れ層は2層、10次元;
(3) 出力次元は N-1=2-1=1次元
と、割と conservative(もといテキトー)に設定した。

結果

  • 景気後退の定義

任意の s_t' \in Sに対し、 \hat{m}_1(s_t) \leq \hat{m}_1(s_t') を満たす  s_t \in S を景気後退期と定義する。

トレース

この場合、景気後退期とは  s_t=1。なお、収束診断をしていない暫定的なものであることに注意。


遷移確率の事後分布の例:


もう少しオサレなアニメーション

ココ
github.com

気になる点

  • 遷移確率の変動が少ない
  • MCMCに時間がかかる
  • VIで試したけど、出力結果を評価しづらい
  • model selectionは難題か

参考

  • Filardo, A. J. (1994). "Business-cycle phases and their transitional dynamics", Journal of Business & Economic Statistics, 12(3), 299–308.
  • Kim, Y. and Kang, K. (2022). "Bayesian Inference of Multivariate Regression Models with Endogenous Markov Regime-Switching Parameters", Journal of Financial Econometrics, 20(3), 391–436.
  • McCracken, M. W. and Ng, S. (2015). "FRED-MD: A Monthly Database for Macroeconomic Research", Working Papers 2015-12, Federal Reserve Bank of St. Louis.