J-REITのキャピタルリターン、インカムリターン、リスクをPlotlyの3Dバブルチャートで比較
このコードは、http://www.japan-reit.com からスクレイピングして取得したデータと株価を元に、J-REIT全銘柄のキャピタルリターン、インカムリターン、リスク(ボラティリティ・標準偏差)を比較できるインタラクティブな3Dバブルチャートを作成します。
完成したものは以下のリンクから見られます。
https://chart-studio.plotly.com/~hiroshimaeasyryo/255/#/
import pandas as pd import numpy as np from bs4 import BeautifulSoup from pandas import Series, DataFrame import requests import urllib import YahooFinanceSpider as y import datetime import plotly.offline as po import plotly.graph_objs as go import plotly.figure_factory as ff import math import codecs import io from datetime import timedelta from plotly.offline import iplot from plotly.graph_objs import Scatter, Figure, Layout, Histogram urls = [] soups = [] base_url = 'http://www.japan-reit.com/meigara/' # 銘柄コード codes = [8951,8952,8953,8954,8955,8956,8957,8958,8960, 8961,8963,8964,8966,8967,8968,8972,8975,8976, 8977,8979,8984,8985,8986,8987,3226,3227,3234, 3249,3269,3278,3279,3281,3282,3283,3287,3290, 3292,3295,3296,3298,3451,3309,3453,3455,3459, 3462,3463,3466,3468,3470,3471,3472,3473,3476, 3478,3481,3487,3488,3492,3493,2971,2972,2979] # 銘柄コードをurlに埋め込む for i in codes: url = urllib.parse.urljoin(base_url, str(i) +"/") urls.append(url) # beautifulsoupでparseする for v in urls: res = requests.get(v) soup = BeautifulSoup(res.content, 'html.parser') soups.append(soup) # titleの取得 titles = [] for t in soups: brand = t.title.string[:-40] titles.append(brand) # 分配金利回りの取得 div_int_rates = [] for d in soups: try: table_d = d.find_all("table", {"class":"meigara-table"})[0] div_int_rate = table_d.find_all("td", {"class":"r"})[2].string[:-9] div_int_rates.append(div_int_rate) except: div_int_rates.append('0') # 時価総額の取得 capitals = [] for c in soups: try: table = c.find_all("table", {"class":"meigara-table"})[0] capital = table.find_all("td", {"class":"r"})[1].string[:-3] capitals.append(capital) except: capitals.append('0') # 時価総額をint型に変換し、バブルサイズをいい感じにするために係数で割って調整 rmcs = [] for rmc in capitals: r = round(float(rmc.replace(',', ''))) / 25000 rmcs.append(r) # DataFrameに格納 df = pd.DataFrame([titles, codes, div_int_rates, capitals]).T.rename( columns={0: '銘柄名', 1: '銘柄コード', 2:'分配金利回り[%]', 3:'時価総額[百万円]'}) # j-reit.comから抽出したデータ jrc_data = df.T # 90日間のデータを使う end_time = datetime.date.today() start_time = end_time - timedelta(days = 90) # クローラーを作る c = y.Crawler() # 終値の取得 prc = [[] for i in range(len(codes))] for i in range(len(codes)): for yf in c.get_price(str(codes[i]), start_time, end_time, y.DAILY): prc[i].append(yf.close) row_df = pd.DataFrame(prc).T # 日付の取得 date = [] for yf in c.get_price(str(codes[0]), start_time, end_time, y.DAILY): date.append(yf.date) # インデックスを日付に、カラムを銘柄名に row_df.index = date row_df.columns = titles # 値をfloat型に df = row_df.astype(float) # キャピタルリターン算出の為変化率のdataframeを用意する pct_df = df.iloc[::-1].pct_change().dropna() # キャピタルリターンを'算術平均'にて算出 # 本来なら'幾何平均'を使うべき annual_return = pct_df.mean() * 3650000 / len(df.index) X = annual_return # インカムリターンの算出(スクレイピングデータより) inc_ret = pd.Series(jrc_data.iloc[2]) Y = inc_ret # リスク(標準偏差)の算出 std = pct_df.std() Z = std # チャートの描画 from plotly.graph_objs import Scatter3d from plotly.offline import init_notebook_mode, iplot init_notebook_mode() # <- Notebook出力にはこの１行が必要 scatter = Scatter3d(x=X, y=Y, z=Z, mode='markers', text=X.index, marker=dict(size=rmcs, # マーカーのサイズ color=Z, # 色分けに使う数値（任意の数値を指定可） colorscale='Magma', # 色のパターン showscale=True) # カラーバーを表示 ) # レイアウトの指定 # バブルサイズの調整 start, end = 15000, 1200000 fig = go.Figure(scatter) fig.update_layout(scene = dict( xaxis_title = 'x:年間平均リターン[%]', yaxis_title = 'y:配当利回り[%]', zaxis_title = 'z:標準偏差' ), width=750, height=650, margin=dict(r=20, b=10, l=10, t=10) ) iplot(fig) # <- Notebookに出力するにはiplot関数を使う<200b>
Plotly Chart Studioでオンライン上で表示
chart studioでオンライン上で表示する為にはchart studioの登録とapi keyの設定が必要になります。
import chart_studio.plotly as py py.plot(fig, filename='JREIT bubble chart')
ぜひ自分の手で動かしてみて欲しい。
Facebookが開発した時系列データ分析ツール「Prophet」を使って東京都新規感染者数の予測をしてみよう
データ分析の醍醐味でもあり目的ともいうべき「予測」にトライする
データ分析のアプローチ方法はいくつもあるのだが、こと「時系列データ」の分析に関しては、どちらかというと「Machine Learningアプローチ」よりもトラディショナルな「統計学アプローチ」の方が精度が高いのではないかという説が、最近気になっている。 おそらくMachine Learningの本業はClassification(つまりクラス分け)であって、ファイナンスデータを予測するには、トラディショナルな統計学アプローチの方が最短ルートなのではないだろうか。
予測するのは東京都新規感染者数
と言いつつ予測するのはファイナンスデータではなく、しかもMLアプローチで攻めていく。
データが正規分布にもとづいているのか、定常性条件を満たしているのかは検証せずに、生きた時系列データを用いて将来を予測する。
hiroshimaeasygogo.hatenadiary.com
東京都の新規感染者数のグラフを描くだけの記事は以前掲載した。
今回はこのコードを生かしていく。
そして、Facebookが開発したお手軽時系列データ分析・予測ツール「Prophet」を使用する。
当然、R派もPython派もプロもアマも使える便利仕様だ。
コードを描く前に是非ともコマンドで pip install fbprophet を叩いて準備しておいてもらいたい。
コードスニペット
import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns sns.set(style="darkgrid") import pandas as pd import numpy as np import statistics import datetime import matplotlib.pyplot as plt import seaborn as sns plt.style.use('ggplot') row_data = pd.read_csv('https://stopcovid19.metro.tokyo.lg.jp/data/130001_tokyo_covid19_patients.csv') by_day = row_data.groupby('公表_年月日') [print(str(bd[0]) + ': ' + str(len(bd[1])) + '人') for bd in by_day] dates = [] counts = [] for bd in by_day: dates.append(bd[0]) counts.append(len(bd[1])) # column名を指定する必要があるのに注意 daily_cases = pd.DataFrame({ 'ds':dates, 'y':counts }) # 日付をdatetimeに変換する datetimes = [] for d in daily_cases.ds.values: datetimes.append(datetime.datetime.strptime(d, '%Y-%m-%d')) daily_cases.ds = datetimes # Prophetを利用した予測 from fbprophet import Prophet from fbprophet.plot import add_changepoints_to_plot # MLでお馴染みパラメータチューニング次第で分析の精度が変わる # トレンドが変化するタイミングを調整することができる model = Prophet( changepoint_range=1, changepoint_prior_scale=5, daily_seasonality=False, weekly_seasonality=True, yearly_seasonality=False ) model.fit(daily_cases) # 7日間の予測をしてみる future_df = model.make_future_dataframe(7, freq='D') forecast_df = model.predict(future_df) # プロットする model.plot(forecast_df) plt.show()
# パラメータチューニングで調整した「トレンドの転換点」を表す
fig = model.plot(forecast_df)
a = add_changepoints_to_plot(fig.gca(), model, forecast_df)
# 全体のトレンドと曜日によるトレンドに分解したものを図示
model.plot_components(forecast_df)
結局明日は何人報告される予測なのか
直近10日の予測データをのぞいてみる。
forecast_df[-10:]
予測値である最終column、"yhat"の値は60.122961人とのこと。
一方 columns[['yhat_lower', 'yhat_upper']] から予測値の95%信頼区間も示されており、39.442601人〜78.499673人だそう(バッファがありすぎでは...)。
95%信頼区間については、あくまでもその95%の区間内に入る平均値から母集団を推計する為のものなので、少し数値とイメージにギャップがあるが。
いずれにせよ、この自分が作ったモデルが合っているか、答え合わせは随時行う必要がある。
また、今回はファイナンスデータではないデータを用いたが、これがファイナンスデータ(夢はもちろん株価ですね)に応用できるか、日々チェックが必要になりそうである。
先に予想しておくと、先人たちが幾度と試みて挫折してきた歴史をみると、多分株価をはじめとする「ランダムウォーク」するデータの予測は不可能と思われる。
株価の予測は何かしらの多変量分析に定性情報となるイベント情報を組み合わせて、伝統的なテクニカル分析を組み合わせるしかないと思われ、しかも都度都度人の手でパラメータチューニングを行う必要がありそうではあるが...
果たしてどうだろうか。
第2波への警戒どころか第1波が抑えられずむしろ加速している現状
COVID-19の感染拡大が止まらない
世界の新規感染者数の7日移動平均グラフをPlotlyにあげようと思ったら、ファイルサイズがいっぱいで怒られてしまった。 そのくらい全世界で感染が広がっているし、日数が長引いている。
やむなくmatplotlibで描画したらlegend(凡例)が長すぎて上図のような不格好な感じになってしまった。
恐ろしいのは、これが累計患者数のグラフではなく新規感染者数のグラフであるということだ。新しく感染した人が、昨日よりも今日、今日よりも明日の方が多い状況が続いている。
世間では「感染拡大第2波」を恐れながらいかに経済活動を再開していくかが模索されているようだが、実際には第1波を抑えることができていないどころか、増加ペースが加速してきている。
今回は先にコードスニペットを掲載する
今回は自分でデータを収集・加工する必要がないシリーズです。
import pandas as pd import numpy as np import datetime import matplotlib.pyplot as plt # 元データの取得 row_data = pd.read_csv('https://covid.ourworldindata.org/data/owid-covid-data.csv') droped_df = row_data.drop( ['iso_code', 'continent', 'total_cases', 'total_deaths','total_cases_per_million', 'new_cases_per_million', 'total_deaths_per_million', 'new_deaths_per_million', 'total_tests', 'new_tests', 'total_tests_per_thousand', 'new_tests_per_thousand', 'new_tests_smoothed', 'new_tests_smoothed_per_thousand', 'tests_units', 'stringency_index', 'population', 'population_density', 'median_age', 'aged_65_older', 'aged_70_older', 'gdp_per_capita', 'extreme_poverty', 'cvd_death_rate', 'diabetes_prevalence', 'female_smokers', 'male_smokers', 'handwashing_facilities', 'hospital_beds_per_thousand'], axis=1) # locationごとにグループ分けしたdf newdf = droped_df.groupby('location') dflist = [] for n, df in newdf: pd.DataFrame(df) dflist.append(df) # dateとnew_casesのみ抜き出す merged_df = pd.merge(dflist[10][['date', 'new_cases']], dflist[0][['date', 'new_cases']], on='date', how='left') for c in range(1, len(dflist)): merged_df = pd.merge(merged_df, dflist[c][['date', 'new_cases']], on='date', how='left') # column名の割当 countrylist = [list(set(dflist[i]['location'])) for i in range(len(dflist))] countries = [country for c in countrylist for country in c] countries[0:0] = ['date', 'dummy_Aus'] # 新規感染者数のdataframe merged_df.columns = countries cases_df = merged_df.drop('dummy_Aus', axis=1) cases_df = cases_df.drop('World', axis=1) # 欠損値補完 cases_df = cases_df.fillna(0) # 7日移動平均のdataframeを用意 day_cas_rol = cases_df.rolling(7).mean().dropna() day_cas_rol.where(day_cas_rol > 0, 0, inplace=True) # matplotlibで描画 fig, ax = plt.subplots(figsize=(4.2,3), dpi=270) day_cas_rol.set_index(cases_df.date[6:]).plot.area( ax=ax, alpha=0.75, lw=0.5, stacked=True) plt.xticks(fontsize=4) plt.yticks(fontsize=4) plt.legend(loc='upper left', bbox_to_anchor=(1, 1)) plt.rc('legend', fontsize=0.1) plt.show()
ヒートアップするUS
感染の中心地はブラジル・チリなどの南米やインド・ロシアなどの新興国に移った。
というイメージはもう古くなってしまった。
この1週間の動きは、北京市場でのクラスター発生、ドイツの食肉処理工場での1500人規模のクラスター発生、US、韓国やフランスなどでも感染再拡大の確認など、「感染拡大&抑制のパイオニア」達の経済活動再開後の再活性化である。
そして中でも特に顕著なのは、トップが「検査体制拡充の結果だ」と胸を張るUSだ。
hiroshimaeasygogo.hatenadiary.com
※前の記事に掲載したコードをもう一回実行すればこのグラフがまた出現します。
州ごとに政策が異なり違う動きをしているのだが、US全体で見ると感染拡大ペースがピーク時に近い傾きになってきている。
physical distancing・マスク着用を遵守できない国民性や、デモで人が密集していることも関係しているのでしょう。
一方、自国民からは批判の嵐なのにも関わらず、他国からはコロナ対策を高く評価されている我が国では...
緊急事態宣言解除から約1ヶ月が経過した頃、国内での新規感染者数のうち半数以上を占めているという東京都ではそのリーダーを選び直す選挙を7/5に控え、現職の小池百合子都知事も完全にそのモードに入ったと思われる。
都内の多くの事業者が経済的に致死量に到達するような壊滅的なダメージを受けた為、「自粛から自衛へ」とシフトチェンジを迫られた。
おそらく選挙が決着するまでは、一旦緩めた感染抑制モードを元に戻すことはできないだろう。
しかし目の敵にされてしまっている「夜の街」を徹底検査した結果、都内の新規感染者数は一定のペースで増加してきている。このペースが継続すれば、都知事選が行われ開票される7/5頃には、緊急事態宣言が発令された4/7の1週間前程度の水準まで新規感染者数が増加する見込みである。
ちなみにいわゆる「3密」の条件が揃っていて、なおかつ「夜の街」以上に関わる人が多いであろう「一般のオフィス」でのクラスター発生も先日確認された。
hiroshimaeasygogo.hatenadiary.com
※このグラフも、前の記事のコードをもう一回実行したら出現します
来週はイベントが盛り沢山
6/30(tue) 中国 製造業PMI(6月)
6/30(tue) US 消費者信頼感指数(6月)
7/01(wed) 日本 日銀短観
7/01(wed) US ADP雇用者数(6月)
7/01(wed) US ISM製造業指数(6月)
7/02(thu) US 失業保険申請件数
7/02(thu) US 雇用統計
それから、週末が独立記念日で米国市場が休場。
ボルトン暴露本や人種差別デモの大統領選への影響、米中貿易合意が本当に"over"ではないのか、またUSの新規感染者数の増加傾向やゴールド市場が沸騰しているのも気持ち悪いし、株高(っぽい)・債券安(っぽい)のも不穏なムードを漂わせている。
前回奇妙な結果だった雇用統計の行方も気になる。
なんとなくデンジャラスなムードの高まりを感じつつ、来週のマーケットを迎えましょう。
コロナショック後の東証REIT指数と各銘柄の相関[Systematic risk-非市場リスク-]
今回はマーケット全体とそれを構成する個別銘柄の関係性の話です
今日は、東証REIT指数とそれを構成している主な銘柄が、コロナショック後どのように相関しているかを調べる。 構成している銘柄が東証REIT指数にある程度相関するのは当たり前なのだが、一部で銘柄固有の動きをすることもあるわけで、それがどのような銘柄に起こっているのかを検証する。 時期をコロナショック後に限定しているのは、コロナヴァイラスがREIT市場にクラッシュを起こしそれまでの「US国債10年利回りとの相関」が崩壊し、Equityとして動くようになったからだ。
大変面倒で申し訳ないが、今回は情報は自分で収集してcsvファイルを用意してほしい。
日本の銘柄データを無料で日本語で取得できる数少ないweb serviceの中で自分が推奨しているのはこちら → https://jp.investing.com/
ここからデータを取得し、こういうcsvファイルを用意する。
検証する銘柄は別に何でも良いのだが、今回はREITをタイプ別に分類した時に各タイプを代表する(東証REIT指数への寄与度が高い)順に選別した。なお、ホテル型といえばジャパン・ホテル・リート投資法人が代表格だが、あえてホテル型に関しては特殊な動きをしていることを期待してインヴィンシブル投資法人をチョイスしてみた。
import pandas as pd import numpy as np import pprint import datetime import matplotlib.pyplot as plt import seaborn as sns from sklearn import linear_model import statsmodels.api as sm clf = linear_model.LinearRegression() ＃ 用意したcsvファイルを読み込む prices = pd.read_csv('reit.csv') # データ型の調整 # 日付をdatetime型に dates=[] for d in prices[prices.columns[0]].values: dates.append(datetime.datetime.strptime(d, '%Y/%m/%d')) prices.date = dates # 東証REIT指数をfloat型に Treit = [] for tr in prices[prices.columns[1]].values: Treit.append(float(tr.replace(',', ''))) prices.TREIT_index = Treit # その他個別銘柄をint型に lis = [[] for i in range(6)] for rs in range(len(prices.columns)-2): for r in prices[prices.columns[rs+2]]: lis[rs].append(int(r.replace(',', ''))) prices[prices.columns[rs+2]]= lis[rs] # indexに日付を指定 prices.index = prices.date prices = prices.drop('date', axis=1)
ここまでで元データが完成。 時系列データの相関を求めるので、元データではなく変化率のdataframeを作る。
# 変化率のdataframe pct_chg = prices.pct_change().dropna() # いったんプロットしてみる fig = plt.figure(figsize=(5,4), dpi=150) plt.plot(pct_chg, alpha=0.7) plt.show() plt.rcParams["font.size"] = 4<200b>
これは何だかさっぱりわからないぞ.....
# seabornで単回帰分析 sns.lmplot(x="TREIT_index", y="R_8951", data=pct_chg, scatter_kws={'color':'indianred'}, line_kws={'color':'indianred'}) # 決定係数だけほしいけどsummary全てを取得できるmethodを作っておく def get_summary(tick1, tick2): model = sm.OLS(tick1, tick2) r = model.fit() print(r.summary2()) # 東証REIT指数と日本ビルファンド投資法人(オフィス系)の相関サマリを表示する get_summary(pct_chg.TREIT_index,pct_chg.R_8951)
Results: Ordinary least squares
================================================================================
Model: OLS Adj. R-squared (uncentered): 0.712
Dependent Variable: TREIT_index AIC: -309.6259
Date: 2020-06-18 21:59 BIC: -307.5316
No. Observations: 60 Log-Likelihood: 155.81
Df Model: 1 F-statistic: 149.7
Df Residuals: 59 Prob (F-statistic): 7.95e-18
R-squared (uncentered): 0.717 Scale: 0.00033052
-------------------------------------------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------
R_8951 0.7855 0.0642 12.2338 0.0000 0.6570 0.9139
--------------------------------------------------------------------------------
Omnibus: 4.890 Durbin-Watson: 1.857
Prob(Omnibus): 0.087 Jarque-Bera (JB): 4.234
Skew: -0.423 Prob(JB): 0.120
Kurtosis: 3.990 Condition No.: 1
================================================================================
summaryの見るべき箇所は、左の中央あたりの
R-squared (uncentered): 0.717という部分。これがいわゆる「決定係数」あるいは「R2」と呼ばれるもので、相関係数を2乗したものになる。
マーケットモデルによれば、
総リスク = 市場リスク(unsystematic risk) + 非市場リスク(systematic risk)
これを変形すると、
1 = R^2(決定係数) + σei/σi(個別銘柄の標準偏差/マーケットの標準偏差)
となるので、
この検証に当てはめてみれば、日本ビルファンド投資法人(8951)の非市場リスク(systematic risk)、つまり「分散効果が効かない固有のリスク」は、 1 - 0.717 = 0.283 ということになる。
では他の銘柄も検証してみましょう。
# 日本プロロジスリート投資法人(物流系) sns.lmplot(x="TREIT_index", y="R_3283", data=pct_chg, scatter_kws={'color':'indianred'}, line_kws={'color':'indianred'}) get_summary(pct_chg.TREIT_index,pct_chg.R_3283) # アドバンス・レジデンス投資法人(住居系) sns.lmplot(x="TREIT_index", y="R_3269", data=pct_chg, scatter_kws={'color':'indianred'}, line_kws={'color':'indianred'}) get_summary(pct_chg.TREIT_index,pct_chg.R_3269) # 日本リテール投資法人(商業系) sns.lmplot(x="TREIT_index", y="R_8953", data=pct_chg, scatter_kws={'color':'indianred'}, line_kws={'color':'indianred'}) get_summary(pct_chg.TREIT_index,pct_chg.R_8953) # インヴィンシブル投資法人(ホテル系) sns.lmplot(x="TREIT_index", y="R_8963", data=pct_chg, scatter_kws={'color':'indianred'}, line_kws={'color':'indianred'}) get_summary(pct_chg.TREIT_index,pct_chg.R_8963) # 野村不動産マスターファンド投資法人(総合系) sns.lmplot(x="TREIT_index", y="R_3462", data=pct_chg, scatter_kws={'color':'indianred'}, line_kws={'color':'indianred'}) get_summary(pct_chg.TREIT_index,pct_chg.R_3462)
# 省略して必要なところだけ記載しています # 日本プロロジスリート投資法人(物流系) R-squared (uncentered): 0.268 # アドバンス・レジデンス投資法人(住居系) R-squared (uncentered): 0.522 # 日本リテール投資法人(商業系) R-squared (uncentered): 0.645 # インヴィンシブル投資法人(ホテル系) R-squared (uncentered): 0.389 # 野村不動産マスターファンド投資法人(総合系) R-squared (uncentered): 0.749
調査結果
・コロナの影響に強いはずの物流銘柄のプロロジが、意外にもsystematic riskが高かった。(むしろ指数を大幅に上回るパフォーマンスなのかもしれない)
・インヴィンシブルはホテル系に分類しているが実は物件の半数近くを住居が占めている。にも関わらず、「ケイマン諸島にリゾートホテルを有している」「単価が低い」「分配金が高い」「所有しているホテルがマイステイズなどリーズナブルなホテルが多い」など特徴的で個人投資家からの人気が高かった。しかしコロナの影響を受け賃料減免の決断を全REITに先駆けて発表し、分配金も今期はほぼゼロに決定した。これによりさらにボラティリティが高まった結果、systematic riskが高まったのでしょう。
・ガチガチオフィス系の筆頭、ビルファンドはコロナショックで懸念が高まるオフィス需要の中で、さすが優良物件(優良借主)を有しているだけあって、systematic riskが比較的低い方だった。最大規模のREITということもあり指数への寄与度がそもそも高いということも、影響しているかもしれない。
・野村はバランスよくかつ数多くポートフォリオを組んでいるおかげなのか、決定係数が高い(分散が図られている)。
これをもとに、好きな銘柄を検証してみると面白いでしょう。
US州ごと新規感染者数トレンド推移をPlotlyで描画する
前回記事のグラフをPlotlyで描画してみた
hiroshimaeasygogo.hatenadiary.com
前回はUSの州ごとの新規感染者数の推移を、季節調整を施したトレンドラインでmatplotlibを使って描画した。
しかしlegend(つまり凡例)が見づらく、どこが増えていてどこが減っているのかmatplotlibではわかりづらいという話をした。
ということで今回はplotlyを使ってグラフを描画してみたので、参考にしてもらえれば幸いだ。
残念ながらstaked area plot(積み上げ面グラフ)の描き方が上手くいかなかったので、stacked bar chart(積み上げ棒グラフ)になってしまったが(実数ではなくトレンドラインを示すのに棒グラフは何となく不自然)、傾向が掴めればそれでよしとしよう。
https://chart-studio.plotly.com/~hiroshimaeasyryo/159.embed
前回触れた通りブログ内に埋め込みができないのでリンクから飛んで見てきてほしい。
コードスニペット
途中まで前回と一緒ですが一応載せておきます。
import pandas as pd import numpy as np import pprint import datetime import matplotlib.pyplot as plt import seaborn as sns row_data = pd.read_csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv') states_df = row_data.groupby('Province_State') df_list = [] for dl in states_df: pd.DataFrame(dl) df_list.append(dl) daily_list = [] for i in range(len(df_list)): daily_list.append(df_list[i][1].T[11:].sum(axis=1)) daily_df = pd.DataFrame(daily_list).T states_list = [] for dl in range(len(df_list)): states_list.append(df_list[dl][0]) daily_df.columns = states_list daily_cases = daily_df.diff().dropna() daily_cases.where(daily_cases > 0, 0, inplace=True) day_cas_rol = daily_cases.rolling(7).mean().dropna() # ここから追記 # plotlyで描画 # chart studioでオンライン上で表示する為にはchart studioの登録とapi keyの設定が必要 import chart_studio.plotly as py import plotly.graph_objects as go from plotly.subplots import make_subplots # オフラインでnotebook上で見る場合はこの2行が必要 import plotly.offline as offline offline.init_notebook_mode(connected=False) traces = [] for t in range(len(day_cas_rol.columns)): traces.append('traces_' + str(t)) for s in range(len(traces)): traces[s] = go.Bar( x=day_cas_rol.index, y=day_cas_rol[day_cas_rol.columns[s]], name=day_cas_rol.columns[s] ) layout = go.Layout(barmode='stack') fig = go.Figure(data=traces, layout=layout) py.plot(fig) # オフラインで見る場合はこちら offline.iplot(fig)
(季節調整値って何やねん...)
本題の前に
もう数回目の記事なので、これまでと少しだけ違った角度から話をしようと思う。
このブログは、果たしてどの層に需要があるのかよくわからないがとりあえずは、
・金融・経済の話題について
・汎用プログラミング言語「Python」を用いて
・実際の経済データを分析しながら
取り上げていきたいと考えている。
その為、「金融や経済に興味がある人」や「Pythonというプログラミング言語に興味がある人」や「実際のデータにもとづいてデータ分析をしたい人・グラフを見たり描いたりして分析っぽいことをするのが好きな人」が読者になることを想定して独り言を綴っている。
これはいわゆる「AND条件」つまり条件全てに当てはまる人を対象にしている、というわけではなく「OR条件」つまりどれかに当てはまる人を対象にしているブログのつもりである。
娯楽のつもりで
自分はとある街の小さな金融機関に勤める「機関投資家の卵」のポジションに昨年就いたばかりであって、上記の各スペシャリストではない。
もしも金融・経済の第一線で活躍したいのなら経済誌を毎日読みこなすなりスペシャリストのドキュメントを読むなりすべきだ。
エンジニア を目指すなら非エンジニアの記事よりもエンジニアの記事を読んだ方が絶対に良い。
データ分析に興味がある・データサイエンティストになりたい人は某Gグル社で働く現役データサイエンティスト・TJO氏のブログを読むなりtwitterをフォローするなりした方が良い。
このブログの内容はそれらと比べると専門性に乏しい。
ただPythonはわかるけど金融はわからん、金融はわかるけど統計学は知らん、
そういう方々が、本来ならリンクしているはずなのに知識が乏しい分野に触れられる機会にできればいいと思ってこれから書こうと考えている。
もちろん、自分も勉強しながら。
娯楽だと思って読んで欲しい。
季節調整値って何やねん...
さて本題に入ると、自分を含め「季節調整値」が何をどう調整したものかわからないまま鵜呑みにしている金融サイドの人や、わからないので信用しない金融サイドの人が、少数かもしれないが一定程度いるような気がしている。
季節調整をしないとどういうことになるかを示したのが以下の図だ。
アメリカの州ごとのCOVID-19の新規感染者数を表した積み上げ面グラフだが、上下が激しく運動していて何だかよくわからない。
これは日本と同様、特定の曜日が他の曜日と比べて検査数が少ない為に、陽性者数の絶対値が少ない傾向がある為にギザつくのだと思われる。
そこで7日間平均を表すトレンドラインをプロットすると、
このように傾向がわかるのだ。
特にファイナンスデータなどで特に多く見られる「時系列データ」は、「特定の曜日」「特定の月日」など、分析にあたってまず真っ先に季節特性を排除する必要がある場合が多い。
本来はこういうdata setならplotlyを推薦したい
だがしかしはっきり言って、matplotlibでこの大きなグラフを描いても同じ色が何個もあってどれがどれだかわからない。
matplotlib以外にもplotlyというインタラクティブなグラフをプロットできるツールがあって、見る方がそちらの方がわかりやすいのだがはてなブログに埋め込みできない。
plotlyのグラフ描画は別に機会を設けるとしよう。
コードスニペット
import pandas as pd import numpy as np import pprint import datetime import matplotlib.pyplot as plt import seaborn as sns row_data = pd.read_csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv') states_df = row_data.groupby('Province_State') df_list = [] for dl in states_df: pd.DataFrame(dl) df_list.append(dl) daily_list = [] for i in range(len(df_list)): daily_list.append(df_list[i][1].T[11:].sum(axis=1)) daily_df = pd.DataFrame(daily_list).T states_list = [] for dl in range(len(df_list)): states_list.append(df_list[dl][0]) daily_df.columns = states_list daily_cases = daily_df.diff().dropna() daily_cases.where(daily_cases > 0, 0, inplace=True) # 何だかわからない、生データのプロット fig, ax = plt.subplots(figsize=(4.2,3), dpi=290) daily_cases.set_index(daily_cases.index).plot.area( ax=ax, alpha=0.75, lw=0.5, stacked=True) plt.xticks(fontsize=4) plt.yticks(fontsize=4) plt.legend(loc='upper left', bbox_to_anchor=(1, 1)) plt.rc('legend', fontsize=2) plt.show() # 季節調整を施したトレンドプロット day_cas_rol = daily_cases.rolling(7).mean().dropna() fig, ax = plt.subplots(figsize=(4.2,3), dpi=300) day_cas_rol.set_index(day_cas_rol.index).plot.area( ax=ax, alpha=0.75, lw=0.5, stacked=True) plt.xticks(fontsize=4) plt.yticks(fontsize=4) plt.legend(loc='upper left', bbox_to_anchor=(1, 1)) plt.rc('legend', fontsize=2) plt.show()
一気に増加してきた東京の新規感染者数
東京アラート解除宣言2日後の急増
東京都の小池百合子知事が6/12、都独自の行動ガイドラインである「東京アラート」を解除することを明らかにした。
「東京アラート」の実効性や目安となる数値の根拠、方向性が曖昧なことなどから一部で批判の声も上がっていたが、おそらく「COVID-19そのものが我々にもたらす影響」よりも、「COVID-19によって経済活動をストップすることの影響」が、遥かに大きそうだということを一旦結論づけての、経済活動再開への舵きりだったと思う。
ヴァイラスは目に見えないのでどの程度感染が拡大しているのかわからず、状況が刻一刻と変化する中で一律の基準を示すということは大変難しいことを考えると、一概にこの判断を真っ向から否定できる自信はない。
この判断が正しかったのかどうなのか、正解はまだ現時点ではわからないというべきでしょう。
そんな中、6/14の新規感染者数は47名にのぼり、久々に40名を超えてきた。どうも盛んに警戒してきた「夜の街」でクラスターが発生し、そこで18名を計上したようだ。
日本の検査体制は相変わらず「クラスター潰し」なのか
なぜ「飲食店」の営業は認められ、「接待を伴う飲食業」が自粛あるいは営業縮小を求められるのか。
そして「接待に利用される料亭」は「接待を伴う飲食業」ではないというガイダンスを後から加えざるを得なくなる、という一連の動きをみていると、一貫性に乏しく合理性に欠けるようにも見える。
そして「３密」の代表格でありライブハウスよりもホストクラブよりも最も多くの人間がそのシチュエーションを経験していると思われる「会社での勤務」がなぜ、経営者への働きかけによって強く自粛の要請(という無言の圧力)がかけられなかったのかを考えると、やはりその方針は合理性に欠けるように見える。
日本の検査数の少なさは感染拡大初期段階から、国内外からも非難轟々だったが、相変わらず現在も感染経路不明者よりもクラスター潰しを重点に検査が行われているように感じる。報道・発表の仕方が「夜の街を控えて」をゴリ押しするのだから、そういうイメージをもたれても仕方があるまい。
今日の47人という発表はブレのようなもの
全国的な緊急事態宣言解除を受けて新規感染者数が再び増加傾向に転じることは、おそらくほぼ全ての国民が予想していたはずである。
しかしグラフのトレンドラインをみてわかるように、急激な上昇傾向に転じたわけでは今のところない。
しかも半数近くが感染経路が判明しているので、いわば「イレギュラーが発生したようなもの」と捉え過度に悲観的にならない方がよいと思われる。少なくとも今の時点では。
だが今後はわからない。引き続きしばらくはしっかり手洗いうがいして、できる範囲で経済を回していく他ない。
コードスニペット
pip install statsmodelsでstatsmodelsをインストールする必要があります。
seasonalをプロットすると季節調整がどのように行われているかグラフで確認できるので、興味がある方はseasonal.plot()で確認してみましょう。
import pandas as pd import numpy as np import datetime import matplotlib.pyplot as plt import seaborn as sns import statsmodels.api as sm plt.style.use('ggplot') row_data = pd.read_csv('https://stopcovid19.metro.tokyo.lg.jp/data/130001_tokyo_covid19_patients.csv') by_day = row_data.groupby('公表_年月日') dates = [] counts = [] for bd in by_day: dates.append(bd[0]) counts.append(len(bd[1])) daily_cases = pd.DataFrame({'counts':counts} ,index=dates) # この2行は最新データを追加する為に記述したもの。普段は不要 case6_14 = pd.Series(47, index=daily_cases.columns, name="2020-06-14") daily_cases = daily_cases.append(case6_14) datetimes = [] for d in daily_cases.index: datetimes.append(datetime.datetime.strptime(d, '%Y-%m-%d')) daily_cases.index = datetimes seasonal = sm.tsa.seasonal_decompose(daily_cases.counts, period=7) daily_cases['trend'] = seasonal.trend # グラフの描画 fig = plt.figure(figsize=(6,4), dpi=200) plt.plot(daily_cases, alpha=0.7) plt.show()
