「PythonのMCMCライブラリPyMC」章

テキスト「PythonのMCMCライブラリPyMC」の執筆者
渡辺祥則 先生

この記事は、テキストの「PythonのMCMCライブラリ PyMC」章のサンプルコード code4.py に掲載された例題４「無題」の実践を取り扱います。
株価のボラティリティを推論している模様です。
コードを解読してPyMCモデルに迫っていきましょう！

はじめに

岩波データサイエンスVol.1の紹介

この記事は書籍「岩波データサイエンス vol.1」（岩波書店、以下「テキスト」と呼びます）の特集記事「ベイズ推論とMCMCのフリーソフト」のベイズモデルを用いて、PyMC Ver.5で「実験的」に実装する様子を描いた統計ドキュメンタリーです。

テキストは、2015年10月に発売され、ベイズモデリングの様々なソフトウェアを用いたモデリング事例を多数掲載し、ベイズモデリングの楽しさを紹介する素晴らしい書籍です。
入門的なモデルから2次階差を取り扱う空間モデルまで、幅広い難易度のモデルを満喫できます！

このシリーズは、テキストのベイズモデルをPyMC Ver.5に書き換えて実践します。

引用表記

この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「岩波データサイエンス vol.1」
第９刷、編者 岩波データサイエンス刊行委員会 　岩波書店

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

PyMC環境の準備

Anacondaを用いる環境構築とGoogle ColaboratoryでPyMCを動かす方法について、次の記事にまとめています。
「PyMCを動かすまでの準備」章をご覧ください。

PythonとPyMCの最新化
テキストのPyMCのバージョンはVer.3です。
また、テキスト発行の2015年から現在2024年までの約10年で、Python及び各種パッケージの利用方法が変わっています。
そこでこの記事では、テキスト掲載のPython＆PyMC3コードを自己流Python＆PyMC5コードに変換することを、実践事項の１つに取り上げようと思います。

例題３ 離散変数のサンプリング

1. 前説

テキストには一切文章が掲載されていないのですが、サンプルコードにひっそり佇む「code4.py」を発見したので、PyMC5で実践したいと思います！
テキストから解説・補助を得ることができないため、コードを解析して「やりたいこと・分析したいこと」を想像しながらの作業となりました。
ミス・手落ちなどはどうぞご容赦くださいね🍀

サンプルコードでは Webサイトから日経225インデックスの日別CSVデータをダウンロードするようになっています。
ただ、現在、このWebサイトにアクセスすることができません。
そこで、Yahoo! ファイナンス から日経225インデックスの株価データをダウンロードして、テキストの時系列モデルを再現してみます。

2. 準備

「PythonのMCMCライブラリPyMC」章で利用するパッケージをインポートします。

### インポート

# ユーティリティ
import pickle

# 数値・確率計算
import pandas as pd
import numpy as np

# PyMC
import pymc as pm
import pytensor.tensor as pt
import arviz as az

# 回帰分析
import scipy.stats as stats
import statsmodels.formula.api as smf
from sklearn.linear_model import LinearRegression

# 分類タスク
from sklearn.metrics import classification_report

# WEBアクセス
import requests

# yfinance株価データ取得
from pandas_datareader import data as pdr
import yfinance as yf

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'

# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')

yfinance パッケージと pandas_datareader パッケージを利用します。
必要に応じてインストールしてください。
以下はAnaconda環境のインストール例です。

### yfinance、pandas-datareaderのインストール
# !conda install -c conda-forge yfinance
# !pip install pandas-datareader

3. データの読み込みと前処理

Yahoo! ファイナンス から日経225の日足データを取得して、pandas_datareader でpandas のデータフレームに格納します。
「# 検索条件の設定」の３つの変数に、銘柄コード、取得開始年月日、取得終了年月日を指定します。

### Yahoo!Financeから日経225株価データを取得

# 検索条件の設定
ticker = '^N225'        # 銘柄コード
start = '2021-01-01'    # 開始日
end = '2024-01-20'      # 終了日

# yfinanceで株価検索してデータフレームに格納
yf.pdr_override()
data4 = pdr.get_data_yahoo(tickers=ticker, start=start, end=end)
display(data4)

# データフレームをcsvファイルに出力
data4.to_csv('./data/N225.csv')

【実行結果】
746日分の日足データを取得しました。

終値（Close）データのみ利用します。
テキストの処理どおり、当日と前日の終値の差＝リターン（株価変動）の対数を算出します。

### データの前処理
close_price = data4['Close']
log_returns = np.log(close_price).diff().dropna()
display(log_returns.to_frame())

【実行結果】
データが１件減りました。

4. データの外観の確認

終値の時系列プロットと株価変動の時系列プロットを描画します。
サンプルコードの「title="Close Price of NIKKEI index"」と「title='return of NIKKEI index close price'」の描画に相当します。

### 日経225の株価と株価変動の推移
# 株価推移
close_price.plot(figsize=(15, 7), xlabel='日付', ylabel='円',
                 title='日経225終値推移')
plt.grid(lw=0.5)
plt.show()
# 株価変動推移
log_returns.plot(figsize=(15, 7), lw=0.7, xlabel='日付', ylabel='株価変動の対数', 
                 title='日経225株価変動 終値ベース・対数')
plt.axhline(0, color='red', lw=0.5)
plt.grid(lw=0.5)
plt.show()

【実行結果】
2023年5月頃から日経225インデックスは上昇基調に入った模様です。
コロナの規制が解除されたころでしょうか。

一方で、株価変動は０付近を上下に行き来しています。
定常過程のように見えます。

階層ベイズモデル

1. 前説

■ モデル数式
テキストのPyMCモデルから数式を逆引きして、若干変更を加えました。
事前分布から尤度関数の順で書きます。

サンプルコードが株価分析用の規定のモデルを適用しているのか、独自のモデルを作成しているのかはよく分かりません。
概観すると次のようになっています。

• 目的変数である株価変動（対数）はスチューデントの$${t}$$分布に従います。

• 注目するパラメータは$${s}$$です。
  $${s}$$はスチューデントの$${t}$$分布の平均パラメータに関わり、指数関数$${\exp(-2s)}$$とされています。
  $${s}$$は対数のさらに対数なのでしょうか・・・・？
  $${\boldsymbol{s}}$$はランダムウォークに従うとしています。
  時系列分析らしさが出ています。

サンプルコードに変更を加えた箇所は、主に、MCMCで分布が収束するように微調整した部分になります。

2. モデルの定義

数式モデルどおり実直にPyMCのモデルを定義しましょう。
また、サンプルコードに合わせて観測データから直近600件を抽出して、モデルに与えることとしています。

### モデルの定義

## データの前処理 直前600日分に絞り込み
log_returns600 =log_returns[-600:] 

## モデルの定義
with pm.Model() as model4:
    
    ### データ関連定義
    ## coordの定義
    # 観測データのインデックス
    model4.add_coord('data', values=log_returns600.index, mutable=True)
    
    ## dataの定義
    # 観測値Y :株価変動(対数)
    Y = pm.ConstantData('Y', value=log_returns600, dims='data')

    ### 事前分布
    # ガウシアンランダムウォーク
    sigmaS = pm.Exponential('sigmaS', lam=1/0.02,
                               transform=pm.distributions.transforms.log)
    init_dist = pm.Normal.dist(mu=0, sigma=0.001)
    s = pm.GaussianRandomWalk('s', mu=0, sigma=sigmaS, init_dist=init_dist,
                              dims='data')  # sigmaを変更 元：sigmaS**-2
    # 尤度のパラメータ
    nu = pm.Exponential('nu', lam=1/10)
    sigmaY = pm.Uniform('sigmaY', lower=0, upper=0.01) # 追加
    
    ### 尤度:スチューデントのt分布
    y = pm.StudentT('y', nu=nu, mu=pt.exp(-2*s), sigma=sigmaY, observed=Y,
                     dims='data')

モデルの内容を表示・可視化してみましょう。
モデル名 model4 とタイプするとモデルの数式の概要が表示されます。

### モデルの表示
model4

【実行結果】

続いてモデルを可視化します。
１行のコードでグラフ形式の可視化ツール「Graphviz」に自動連携して、変数の関係を図示してくれます！

### モデルの可視化
pm.model_to_graphviz(model4)

【実行結果】
データやパラメータの繋がりが示されます。
丸い形状の変数には「確率分布」が併記されます。
データ値を示す変数はグレー色が付いています。

3. MCMCの実行と収束の確認

■ MCMC
いよいよベイズプログラミングの真骨頂 MCMC の実行です。
マルコフ連鎖=chainsを4本、バーンイン期間=tuneを1000、利用するサンプル=drawを１chainあたり1000に設定して、合計4000個の事後分布からのサンプリングデータを生成します。テキストの指定数とは異なっています。

サンプリングデータは idata4 に格納されます。

サンプル方法に numpyro を指定しています。処理速度が早くなります。
numpyro を使わない場合は「nuts_sampler='numpyro',」を削除します。

### 事後分布からのサンプリング 1分50秒
with model4:
    idata4 = pm.sample(draws=1000, tune=1000, chains=4, random_seed=1234,
                       nuts_sampler='numpyro', target_accept=0.9)

【実行結果】
処理時間は約1分50秒です。

■ $${\hat{R}}$$で収束確認
収束の確認をします。
指標$${\hat{R}}$$（あーるはっと）の値が$${1.1}$$以下のとき、収束したこととします。
次のコードではしきい値を$${1.06}$$に設定して、$${\hat{R}>1.06}$$のパラメータの個数をカウントします。
個数が０ならば、$${\hat{R} \leq1.06}$$なので収束したとみなします。

### r_hat>1.1の確認
# 設定
idata_in = idata4        # idata名
threshold = 1.06         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
$${\hat{R}>1.06}$$のパラメータは０個でした。
$${\hat{R} \leq1.1}$$なので収束したとみなします。

■ $${\hat{R}}$$の値把握と事後統計量の表示

### 事後統計量の表示
var_names = ['nu', 'sigmaS','sigmaY', 's']
pm.summary(idata4, hdi_prob=0.95, var_names=var_names)

【実行結果】
MCMCでサンプリングデータを生成したすべてのパラメータについて、統計量や診断情報が一覧表示されました。

■ トレースプロットで収束確認
トレースプロットを描画します。
テキストの図7に該当します。
こちらの図でも収束の確認を行えます。

### トレースプロットの描画
pm.plot_trace(idata4, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】
左の曲線は、４本のマルコフ連鎖の分布です。
４本の曲線がほぼ同じ分布を描いているので、サンプリングデータが同一の分布から生成されたことが推察できます。
右のノイズみたいな図は、均等に乱雑に描画されています。
この状態は収束していると推察できます。
均等・乱雑ではなくて、何らかの傾向を示すような描画の場合は収束を疑います。

$${\hat{R}}$$とトレースプロットから、分布の収束が確認できましたので、MCMCのサンプリングデータを用いて、推論を継続できます！
やったね！

4. 分析

サンプルコードの２つのプロットを描画してみます。
モデルやパラメータの意味合いをあまり理解できていないので、プロットすることが限界、結果を分析するには至っていません（汗）
詳細な分析はやりません（やれません）。

■ パラメータ$${s}$$の推論値のプロット
テキストの「'log volatility'」に相当します。
$${s}$$は「対数ボラティリティ」だったのでしょうか？？？

### ボラティリティ？の描画 ※sの事後分布サンプリングデータのプロット

## データ準備
# sの事後分布サンプリングデータの取り出し
s_post_sel = idata4.posterior.s.stack(sample=('chain', 'draw')).T[::10].T

## 描画
# 描画領域の指定
plt.figure(figsize=(10, 3))
ax = plt.subplot()
# プロット
ax.plot(log_returns600.index, s_post_sel, color='blue', alpha=0.01)
# 修飾
ax.set(title='log volatility?', xlabel='日付')
plt.grid(lw=0.5);

【実行結果】
左下から右へ煙のように青い靄が広がっています。
時間が経過するほどボラティリティが大きくなることを示しているのかもしれません。

■ 株価変動とボラティリティの描画
冒頭で示した株価変動データと上のボラティリティデータを同時に描画します。
サンプルコードのボラティリティ「exp(s)」のプロットを変形し、さらに株価変動を並べてみました。

### 株価変動とボラティリティの描画
# ボラティリティsを株価変動と同じ単位に変換
exps = np.exp(-2*s_post_sel)
# 描画領域の指定
plt.figure(figsize=(10, 3))
ax = plt.subplot()
# 株価変動の描画
ax.plot(log_returns600.index[2:], log_returns600[2:].values, lw=0.5,
        label='株価変動（対数）')
# ボラティリティの描画
ax.plot(log_returns600.index[2:], exps[2:], 'r', alpha=0.05)
ax.plot(log_returns600.index[2:], -exps[2:], 'r', alpha=0.05)
# ボラティリティの凡例作成用のダミー
ax.plot([], [], 'r', lw=0.4, label='ボラティリティ（補正後）')
ax.set(xlabel='日付', ylabel='株価変動（対数）',
       title='日経225株価変動とボラティリティ')
ax.legend()
ax.grid(lw=0.5);

【実行結果】
赤い線がボラティリティです。上下に広がっている時期がボラティリティの高い時期だったのかもしれません。

うまく纏められませんですが、サンプルコードのような株価データの変動の分析にもベイズが使えるということは分かったような気がします。

5. 推論データの保存

pickle で MCMCサンプリングデータ idata4 を保存できます。

### 推論データの保存
file = r'idata4_ch1.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata4, f)

次のコードで保存したファイルを読み込みできます。

### 推論データの読み込み
file = r'idata4_ch1.pkl'
with open(file, 'rb') as f:
    idata4_ch1_load = pickle.load(f)

ベイズ記事は以上です。

次回予告

「階層ベイズ最初の一歩」章
　モデル１「直線的あてはめ」

結び

株価分析のベイズモデリングに興味を持ちました。
ただテキストにベイズモデル自体の説明や株価変動の分析方法が掲載されておらず、悶々としています。
ベイズの株価分析の基本的な型があったら知りたいです（緩募）。

おわり

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。