ランダム変数と全要素生産性#
import japanize_matplotlib
import numpy as np
import pandas as pd
import py4macro
import statsmodels.formula.api as smf
# 警告メッセージを非表示
import warnings
warnings.filterwarnings("ignore")
はじめに#
この章の目的は２つある。 第一に，マクロ経済分析におけるランダム変数の役割と変数の持続性について考える。 第二に，景気循環を理解する上で全要素生産性の役割が強調される研究があるが，その導入として全要素生産性の特徴について検討する。またソロー・モデルを使い，全要素生産性が確率的に変動する簡単な景気循環モデルのシミュレーションを通してマクロ変数の特徴を考察する。
ランダム変数の役割#
# データの読み込み
df = py4macro.data('jpn-q')
# トレンドからの％乖離のデータの作成
df['gdp_cycle'] = np.log( df['gdp'] ) - py4macro.trend( np.log(df['gdp']) )
# プロット
ax = ( 100 * df['gdp_cycle'] ).plot(title='GDPのトレンドからの乖離（％）')
ax.axhline(0, c='red')
pass
このようなGDPの変動を説明する方法として２つある。
決定的（deterministic）な過程で生成された時系列として捉える。
ランダム（stochastic）な過程で生成された時系列として捉える。
１の考えによると，差分方程式\(x_{t+1}=f(x_t)\)のように来期のGDPは今期のGDPに依存している。マクロ変数の持続性を捉えるには必要不可欠な特徴と言えるだろう。しかし，今期のGDP（\(x_t\)）が与えられると来期のGDP（\(x_{t+1}\)）は自ずと決定さることになり，来期のGDPの予測は簡単なものとなる。もしそうであれば政策運営は非常に簡単であろうが，現実はそうではない。そういう意味では，１は現実を十分に捉えることができていない。
一方，２の考えはランダム変数の実現値の連続としてGDPが観測され，今期のGDPが与えられても来期のGDPにはランダムな要素があるため，予測が難しいという特徴がある。この特徴こそが，マクロ変数の性質を捉えるには必要な要素であり，上のプロットに現れていると考えられる。主流のマクロ経済学では，消費者や企業の行動はフォーワード・ルッキングであり（将来を見据えた最適化行動であり），且つ経済全体ではランダムな（不確実な）要素が重要な役割を果たしているという考え基づいている。そのような分析枠組みの中で，景気循環のメカニズムを解明することが目的となっている。
以下では，まずランダム変数生成過程を考え，短期分析における全要素生産性の役割について説明する。またSolowモデルを使い簡単な景気循環モデルを展開する。
ホワイト・ノイズ#
説明#
時系列のランダム変数\(\varepsilon_t\)，\(t=0,1,2,3,\cdots\)を考えよう。例えば，レストランの経営者の収入。ビジネスにはリスク（競争相手の出現やコロナ感染症問題）があるため変動すると考えるのが自然である。\(\varepsilon_t\)は\(t\)毎にある分布から抽出されると考えることができる。次の３つの性質を満たしたランダム変数をホワイト・ノイズ（White Noise）と呼ぶ。
平均は
0：\(\quad\text{E}\left[\varepsilon_t\right]=0\)分散は一定：\(\text{ E}\left[\varepsilon_t^2\right]=\sigma^2\)（\(\sigma^2\)に\(t\)の添字はない）
自己共分散は
0：\(\text{ E}\left[\varepsilon_t \varepsilon_{t-s}\right]=0\)（全ての\(s\ne 0\)に対して）。類似する概念に自己相関係数があり，値が\([-1,1]\)になるように標準化されている。（ホワイト・ノイズの自己相関係数は0となる。）\[ \text{自己相関係数}=\frac{\text{ E}\left[\varepsilon_t \varepsilon_{t-s}\right]} {\sqrt{\text{ E}\left[\varepsilon_t^2\right]\text{ E}\left[\varepsilon_{t-s}^2\right]}} \]
＜コメント＞
平均0と分散が\(\sigma^2\)のホワイト・ノイズを次のように表記する。
\[ \varepsilon_t\sim\textit{WN}(0,\sigma^2) \]ホワイト・ノイズの例：
平均
0，分散\(\sigma^2\)の正規分布から抽出され\(\varepsilon_t\)最小値
10，最大値100からから抽出され\(\varepsilon_t\)（正規分布や一様分布ではない分布でもホワイト・ノイズになり得る）
ホワイト・ノイズは独立同分布（Indepent and Identically Distributed）の１つである。
以下では実際にランダム変数を発生させ，そのプロットと統計的な特徴について考察する。その準備として次のコードを実行しよう。
rng = np.random.default_rng()
Note
NumPyを使いランダム変数を生成する方法については，「経済学のためのPython入門」のNumpy: ランダム変数を参考にしてください。
図示#
標準正規分布からのランダム変数をn個抽出しプロットしてみよう。
n = 100
vals = rng.normal(size=n)
df_wn = pd.DataFrame({'White_Noise':vals})
df_wn.plot(marker='.')
pass
２つの特徴がある。
平均
0：0を中心に周辺を動いている。分散一定：
0から外れても0に戻っている。
平均と分散を計算してみよう。
vals.mean(), vals.var()
(0.01736777770114493, 0.9804631771614523)
サンプル統計量であるため誤差が発生していることがわかる。
２つの特徴をを確認するためにPandasのplot()を使ってヒストグラムを描いてみる。
df_wn.plot(kind='hist', bins=30)
pass
概ね0を中心に左右対象に動いている。即ち，平均0を反映している。また，0から離れている観測値が少ないことが分かる。これは0方向に戻ることを示している。これは分散が一定であるためである。この性質はnを500や1000などの大きな数字に設定するより一層分かりやすいだろう。
自己共分散（自己相関）#
次に１期違いの自己共分散\(\text{E}[\varepsilon_t \varepsilon_{t-1}]\)を考えよう。
current = vals[:-1] # 最後の要素を除外
one_period_ahead = vals[1:] # 0番目の要素を除外
この２つの変数を使って散布図を描いてみる。Pandasのplot()で横軸・縦軸を指定し，引数kindで散布図を指定するだけである。
# DataFrameの作成
df_wn2 = pd.DataFrame({'current':current,'one_period_ahead':one_period_ahead})
# 散布図のプロット
df_wn2.plot(x='current', y='one_period_ahead', kind='scatter')
pass
自己共分散がゼロであれば，ランダムに散らばっているはずであり，何らかのパターンも確認できないはずである。
ここで異なるプロットの方法を紹介する。Pandasのサブパッケージplottingにあるlag_plot()関数を使うと共分散の強さを示す図を簡単に作成することが可能となる。
必須の引数：１つの列
オプションの引数：
lagは時間ラグ（デフォルトは1）
pd.plotting.lag_plot(df_wn)
pass
横軸：
t期の値縦軸：
t+1期の値
２つの図は同じであることが確認できる。
分散・自己共分散を計算してみよう。
np.cov(current, one_period_ahead)
array([[ 0.99363141, -0.04304577],
[-0.04304577, 0.99548054]])
右上と左下の値が自己共分散であり、非常に小さな値である。自己共分散は0という仮定を反映している。（左上はcurrentの分散であり、右下はone_period_aheadの分散である。）
次にメソッドautocorr()を使い自己相関係数を計算しよう。
df_wn['White_Noise'].autocorr()
-0.04328141327500157
自己共分散0の仮定の反映である。この結果からわかることは，ホワイト・ノイズだけではマクロ変数の持続性を説明できないということである。
自己回帰モデル：AR(1)#
説明#
次に確率過程として自己回帰モデルを考えよう。英語でAugoregressive Modelと呼ばれ，AR(1)と表記されるモデルは次式で表される。（AR(1)の(1)は右辺と左辺の変数は１期しか違わないことを表している。）
ここで
\(y_t\)：\(t\)期の\(y\)の値
\(-1<\rho<1\)
\(\varepsilon_t\sim WN(0,\sigma^2)\)（ホワイト・ノイズ）
\(\varepsilon_t\)を所与とすると\(y_t\)の差分方程式となっている。\(-1<\rho<1\)となっているため，安定的な過程であることがわかる。即ち，\(\varepsilon_t=0\)であれば、\(y_t=\)は定常状態である\(0\)に近づいて行くことになる。しかしホワイト・ノイズである\(\varepsilon_t\)により、毎期確率的なショックが発生し\(y_t\)を定常状態から乖離させることになる。また式(148)からわかるように，今期の値\(y_t\)は前期の値\(y_{t-1}\)に依存しており，その依存度はパラメータ\(\rho\)によって決定される。この\(\rho\)こそが\(y_t\)の持続性の強さを決定することになり，自己回帰モデル(148)の自己相関係数は\(\rho\)と等しくなる。
３つの例：持続性の違い#
直感的にマクロ変数の特徴である持続性とは，前期の値が今期の値にどれだけ影響を持っているかを示す。\(\rho\)の値を変えて持続性の違いを視覚的に確認するための関数を作成しよう。
def ar1_model(rho, T=100):
"""引数：
rho: AR(1)の持続性を捉えるパラメータ
T: シミュレーションの回数（デフォルト：0）
戻り値：
matplotlibの図を示す
自己相関係数の値を表示する"""
y0 = 0 # 1
y_list = [y0] # 2
y = y0 # 3
for t in range(1,T): # 4
e = rng.normal() # 5
y = rho*y + e # 6
y_list.append(y) # 7
df_ar1 = pd.Series(y_list) # 8
ac = df_ar1.autocorr() # 9
ax_ = df_ar1.plot(marker='.') # 10
ax_.set_title(fr'$\rho$={rho} 自己相関係数：{ac:.3f}',size=20) # 11
ax_.axhline(0, c='red') # 12
コードの説明
yの初期値を設定。yの値を格納するリスト。初期値を入れてある。左辺の
yはforループで使うアップデート用のyであり，初期値を割り当てている。T期間のforループの開始。変数tは(4)と(5)に入っていないので単にループの回数を数えている。range(1,T)の1はループの計算が\(t=1\)期のyの計算から始まるため。\(t=0\)期のyは(1)で与えられている。
t期のホワイト・ノイズの生成。右辺の
eはt期のホワイトノイズであり，y前期のyとなっている。y_listにyを追加。y_listを使ってSeriesを作成しdf_ar1に割り当てる。.autocorr()を使い己相関係数を計算しacに割り当てる。yのプロットmarker='.'はデータのマーカーを点に設定。図の「軸」を
ax_に割り当てる。
図のタイトルの設定
frのfはf-stringを表し，{rho}と{ac:.3f}に数値を代入する。:.3fは小数点第三位までの表示を設定frのrは$\rho$をギリシャ文字に変換するためののも。$\rho$はLaTeXのコード。
.axhline()は横線をひくメソッド。
ar1_model(0.1)
この例では\(\rho\)の値が低いため，\(y\)は定常状態の0に直ぐに戻ろうとする力が強い。従って、前期の値の今期の値に対する影響力が小さいため、持続性が低い。
ar1_model(0.5)
\(\rho\)の値が高くなると、定常状態の\(0\)に戻ろうとする作用が弱くなり、持続性が強くなることがわかる。
ar1_model(0.9)
\(\rho\)が非常に高いため、持続性も非常に強くなっている。即ち、今期の値は前期の値に対する依存度が大きい。
全要素生産性#
説明#
マクロ変数の特徴を捉えるにはホワイト・ノイズよりもAR(1)の方がより適しているということが分かった。では，どの様な形でAR(1)をマクロ・モデルに導入すれば良いのだろうか。現在盛んに行われている景気循環研究はDSGEモデル（Dynamic Stochastic General Equilibrium Model）と呼ばれるモデルに基づいており、DSGEは名前が示すように次のような特徴がある。
時間軸に沿ってマクロ変数の変化を考察する動学的な特徴
確率的な要因の重要な役割
合理的期待に基づく将来を見据えた消費者と企業の最適行動
一般均衡のモデル
そして，DSGE研究の先駆けとなったモデルが実物的景気循環モデル（Real Business Cycles Model; RBCモデル）と呼ばれる理論である。詳細については後の章に譲るが，基本的なアイデアは簡単で全要素生産性（TFP）がAR(1)に従って確率的に変動し，それが消費者と企業の経済行動に影響することにより景気循環が発生するという理論である。ここではデータを使い全要素生産性の性質を考察する。RBCモデルで需要な役割を果たす労働時間についてもデータを使い特徴を探ることにする。
データ#
次の生産関数を仮定する。
\(Y_t\)：GDP
\(K_t\)：資本ストック
\(H_t\)：総労働時間（労働者数\(\times\)労働者平均労働時間）
\(A_t\)：TFP
TFPについて解くと次式を得る。
景気循環を考えているため、TFPを計算する際には次の点を考慮することが望ましい。
年次データではなく四半期データを使う（年次データでは短期的な変動を捉えることができない）。
\(H_t\)は労働者数ではなく，総労働時間とする（短期的な労働供給の変化を捉えるため）。
\(K_t\)は資本ストックの稼働率を考慮したデータとする（短期的に実際に生産に使われた資本ストックの変化を捉えるため）。
ここでは１と２のみを考慮したデータを使いTFPの変動を考えることにする。これにより資本ストックの稼働率はTFPの一部となってしまうが，RBCモデルの基本的なアイデアは変わらない。
上で読み込んだデータdfには次の変数が含まれている（詳細はGDP：水準・トレンド・変動を参照）。
期間：1980年Q1〜2021年Q4
gdp：国内総生産（支出側）capital：資本ストックemployed：就業者hours：労働者1人あたり月平均就業時間（以下では「平均労働時間」と呼ぶ）total_hours：総労働時間（hours\(\times\)employed）
労働時間と雇用の特徴#
平均労働時間hours，就業者数employed，総労働時間total_hoursの特徴を考えてみる。まずサイクルを計算し，それぞれの変数を図示しよう。（employedについては対数化した就業者数のトレンドを計算しても良いが，結果は大きく変わらないため，ここでは対数を使わずに計算する。）
df['hours_trend'] = py4macro.trend(df['hours'])
df['hours_cycle'] = np.log( df['hours']/df['hours_trend'] )
df['employed_trend'] = py4macro.trend(df['employed'])
df['employed_cycle'] = np.log( df['employed']/df['employed_trend'] )
df['total_hours_trend'] = py4macro.trend(df['total_hours'])
df['total_hours_cycle'] = np.log( df['total_hours']/df['total_hours_trend'] )
最初に，平均労働時間を考えよう。
py4macro.data('jpn-q',description=True)
| `gdp`: 国内総生産（GDP）
| `consumption`: 消費
| `investment`: 投資
| `government`: 政府支出
| `exports`: 輸出
| `imports`: 輸入
| `capital`: 資本ストック
| `employed`: 就業者数
| `unemployed`: 失業者数
| `unemployment_rate`: 失業率
| `hours`: 労働者一人当たり月平均労働時間
| `total_hours`: 月平均総労働時間（`employed`X`hours`）
| `inflation`: インフレ率
| `price`: 消費者物価指数
| `deflator`: GDPデフレーター
|
| * 四半期データ
|
| ＜出典＞
| GDPとその構成要素
| * 1994年Q1~2021年Q4
| * 実額・四半期・実質季節調整系列（年換算）
| * 2015暦年（平成27年）連鎖価格
| * 単位：10億円
| * 国民経済計算（GDP統計）
| * 1980年Q1~1993年Q4
| * 実額・四半期・実質季節調整系列（年換算）
| * 平成27年基準支出側GDP系列簡易遡及（参考系列であり上のデータと接続可能）
| * 単位：10億円
| * 国民経済計算（GDP統計）
|
| 実質資本ストック
| * 1994年Q1~2019年Q4
| * 平成25年基準
| * 単位：10億円
| * 国民経済計算（GDP統計）
| * 1980年Q1~1993年Q4
| * 平成25年基準遡及系列
| * 単位：10億円
| * 国民経済計算（GDP統計）
|
| 就業者数，失業者数，失業率
| * 総務省「労働力調査」
| * 単位：万人，％
|
| 労働者一人当たり月平均労働時間
| * 厚生労働省「毎月勤労統計調査」
| * 30 人以上(一般・パート)、月間実労働時間(総実労働時間)
| * 2020年の平均を100に基準化
|
| インフレ率
| * 景気動向指数（速報，改訂値，月次，原数値の前年同月比）から四半期平均として計算
|
| 消費者物価指数
| * 2020年基準
| * 「中分類指数（全国）＜時系列表＞【月次】」の四半期平均として計算
| * 簡便的に移動平均を使い季節調整を施している
|
| GDPデフレーター
| * 1994年Q1~2019年Q4
| * 2015年（平成27年）基準
| * 季節調整系列
| * 1980年Q1~1993年Q4
| * 2015年（平成27年）基準遡及系列
| * 季節調整系列
df[['hours','hours_trend']].plot()
pass
平均労働時間は長期的に減少傾向にある。2020年の平均が100に標準化しているので，1980年代には約25％長かったことを示している。次に，就業者数を図示する。
df[['employed','employed_trend']].plot()
pass
一方，雇用者数は増加傾向にある。女性の労働市場参加率の増加や，雇用形態の変化（非正規など）の影響と考えられる。上の２つの変化を反映したのが総労働時間の変化である。
df[['total_hours','total_hours_trend']].plot()
pass
バブル崩壊後は減少傾向にあるが，過去１０年間に持ち直して来ている。これは雇用の拡大の要因が大きい。
雇用と労働時間の変動（トレンドからの乖離率）を比べてみよう。
ax_ = df[['hours_cycle','employed_cycle']].plot()
ax_.axhline(0, c='red', lw=0.5)
pass
次の２つの点を指摘できる。
労働時間の方が変動が大きい。即ち，総労働時間の調整は就業者数よりも労働時間を通して行われていることが伺える。
就業者数の変化は、労働時間の変化を後追いする傾向がある。これは就業者の調整費用がより高いためであろう。
自己相関係数を計算してみる。
cycle_list = ['hours_cycle','employed_cycle','total_hours_cycle']
var_list = ['平均労働時間','就業者数','総労働時間']
print(' 自己相関係数\n----------------------')
for c, v in zip(cycle_list,var_list):
ac = df[c].autocorr()
print(f'{v:<5}\t{ac:.3f}') # 1
自己相関係数
----------------------
平均労働時間 0.454
就業者数 0.833
総労働時間 0.574
コードの説明
:<5はf-stringで代入する変数vのスペースを5に設定し左寄せ<にすることを指定している。また:.3fは小数点第三位までの表示を指定している。
就業者数がより大きな自己相関係数の値となっており，就業者での調整により時間がかかるためである。即ち，今期の就業者数は前期の就業者数に依存するところが大きいという意味である。
次に，GDPのトレンドからの乖離との相関係数を計算してみる。
print('GDPサイクルとの相関係数')
print('-------------------------')
for c, v in zip(cycle_list,var_list):
corr = df[['gdp_cycle',c]].corr().iloc[0,1]
print(f'{v:<9}\t{corr:.3f}')
GDPサイクルとの相関係数
-------------------------
平均労働時間 0.723
就業者数 0.482
総労働時間 0.814
やはり非常に強い正の相関があり，平均労働時間の影響が強いようである。
TFPの特徴#
TFPの水準，トレンド，変動（サイクル）を計算するが，資本の所得分配率を次のように仮定する。
a = 0.36
# 全要素生産性（対数）の計算
df['tfp_log'] = np.log( df['gdp']/( df['capital']**a * df['total_hours']**(1-a) ) )
# 全要素生産性のトレンド（対数）の計算
df['tfp_log_trend'] = py4macro.trend( df['tfp_log'] )
# 全要素生産性のトレンドからの乖離率の計算
df['tfp_cycle'] = df['tfp_log'] - df['tfp_log_trend']
# 全要素生産性とトレンドのプロット
ax = df[['tfp_log','tfp_log_trend']].plot()
ax.set_title('全要素生産性とトレンド（対数）', size=20)
pass
全要素生産性の変動をプロットしよう。
ax = (100 * df['tfp_cycle']).plot()
ax.axhline(0, c='red')
ax.set_title('全要素生産性の変動（％）', size=20)
pass
この図から全要素生産性の変動は大きいことが分かる。変動の最小値と最大値を確認してみよう。
df['tfp_cycle'].min(), df['tfp_cycle'].max()
(-0.04527476816105558, 0.01918027625853458)
正の乖離は2％近くあり，負の乖離は4％以上となる。リーマンショック時を除くと正も負も絶対値で概ね2％以内に収まっている。
GDPと一緒にプロットしてみよう。
ax_ = (100 * df[['gdp_cycle','tfp_cycle']]).plot()
ax_.set_title('GDPとTFPの変動（％）',size=20)
pass
同じ方向に動いているのが確認できる。GDPの変動との相関係数を計算してみると，非常に高いことがわかる。
df[['gdp_cycle','tfp_cycle']].corr().iloc[0,1]
0.8968819031634699
以前の計算結果を使うと，この数字はGDPのどの構成要素よりも高いことが確認できる。次に，自己相関係数を計算してみよう。
df['tfp_cycle'].autocorr()
0.5097428451922277
gdp_cycle程ではないが、強い持続性があることが確認できる。
df['gdp_cycle'].autocorr()
0.6808888095459286
＜これらの結果の含意＞
全要素生産性の変動は大きく持続性も高い。またGDPとの相関性も大きい。これらから全要素生産性が景気循環を引き起こす要因ではないかということが示唆される。全要素生産性は資本や労働時間で説明できないGDPの要素であり，そのような「その他」の部分に景気循環を引き起こすランダムな要素が隠されているかも知れない，という考え方である。
AR(1)としてのTFP#
全要素生産性を景気循環の要因として捉えるために，全要素生産性をAR(1)としてモデル化してみようというのがこの節の目的である。そのためにstatsmodelsを使い回帰分析をしてみよう。まず１期ずらした変数を作成する。
df['tfp_cycle_lag'] = df['tfp_cycle'].shift()
コードの説明
shift()はSeriesを１期ずらすメソッドである。
df[['tfp_cycle','tfp_cycle_lag']].head()
| tfp_cycle | tfp_cycle_lag | |
|---|---|---|
| 1980-03-31 | 0.003481 | NaN |
| 1980-06-30 | -0.011712 | 0.003481 |
| 1980-09-30 | 0.000090 | -0.011712 |
| 1980-12-31 | 0.011619 | 0.000090 |
| 1981-03-31 | 0.010427 | 0.011619 |
列tfp_cycle_lagは列tfp_cycleを１期ずらしていることが確認できる。次に回帰分析をおこなう。
res_tfp = smf.ols('tfp_cycle ~ -1 + tfp_cycle_lag', data=df).fit()
rho = res_tfp.params.iloc[0]
rho
0.5109092142476828
コードの説明
回帰式の中に-1とあるが、定数項を省く場合に指定する。
自己相関係数と殆ど変わらない値である。２つの変数の散布図を表示してみよう。
df.plot(x='tfp_cycle_lag', y='tfp_cycle', kind='scatter')
pass
正の相関が図でも確認できる。次に\(\epsilon_t\)について考える。この変数はホワイト・ノイズであり、TFPショックの「源泉」である。平均0の正規分布としてモデル化するが，上の回帰式の残差の標準偏差をTFPの標準偏差の推定値として考えよう。
sigma = res_tfp.resid.std()
sigma
0.00845643508172743
コードの説明
res_tfpの属性.residは残差を返している。.residはSeriesであり，その属性std()を使い標本標準偏差を計算している（デフォルトではddof=1)。
初期値は0としてrhoとsigmaを使ってAR(1)の値を計算してプロットしてみよう。
T = 100 # 標本の大きさ
A = 0 # 初期値
A_list = [A] # 値を格納するリスト
for i in range(T):
e = np.random.normal(loc=0,scale=sigma)
A = rho*A+e
A_list.append(A)
df_A = pd.Series(A_list) # Seriesの作成
ax_ = df_A.plot() # 図示
ax_.axhline(0, c='red') # 0の平行線
ac_tfp = df_A.autocorr() # 自己相関係数の計算
print(f'自己相関係数：{ac_tfp:.3f}')
自己相関係数：0.512
シミュレーションを行う度に図と自己相関係数は異なることになる。
確率的ソロー・モデル#
説明#
上で日本のデータを使い全要素生産性のrhoとsigmaを計算しシミュレーションを行なったが，そのような全要素生産性の変動が景気循環の主な要因と考えるのがRBCモデルと呼ばれるものである。RBCモデルの詳細については後の章で検討するが，ここではその導入としてソロー・モデルにAR(1)としてのTFPを組み込んでシミュレーションをおこない，計算結果をデータと比べてどの程度整合性があるかを検討する。
モデルの均衡式として以下を使う。
生産関数：\(Y_t=A_tB_tK_t^aH^{1-a}\)
\(A_tB_t\)が全要素生産性
\(A_t\)は変動を捉える項
\[\begin{split} \begin{align*} \log(A_{t})&=\rho\log(A_{t-1})+\varepsilon_t \\ \varepsilon_t&\sim \textit{WH}(0,\sigma^2) \end{align*} \end{split}\]\(B_t\)はトレンドを捉える項
\[ B_t=B_0(1+g)^t \]\(H\)は一定と仮定する
資本の蓄積方程式：\(K_{t+1}=sY_t+(1-d)K_t\)
消費：\(C_t=(1-s)Y_t\)
投資（貯蓄）：\(I_t=sY_t\)
これらの式から次の2式で\(K_t\)と\(A_t\)の値が計算できる。
変数の初期値を次のように仮定する。
\(A_0=1\)
\(K_0=\) 1980年の資本ストックの平均
\(H\)は次の値として一定とする。
\(H=\) 1980年の総労働時間の平均
パラメータの値を次を仮定する。
\(\rho\)と\(\sigma\)は上で計算した
rhoとsigmaの値\(g=\) 1980年から最後の年までの平均四半期成長率
\(s=\) 1980年から最後の年までのGDPに占める投資の割合の平均
\(d=0.025\)：年率約10％を想定
\(a=0.36\)：資本の所得割合
シミュレーションの結果を次の変数を含むDataFrameにまとめる。
\(Y\)，\(C\)，\(I\)，\(K\)，\(A\)の水準の系列
\(Y\)，\(C\)，\(I\)，\(K\)，\(A\)のトレンドの系列
\(Y\)，\(C\)，\(I\)，\(K\)，\(A\)の変動（サイクル）の系列
初期値の計算#
まずGDP，資本，総労働時間の1980年と2019年の平均を計算する。
var_list = ['gdp','capital','total_hours'] #1
year_list = [df.index[0].year,df.index[-1].year] #2
gdp_dict = {} #3
K_dict = {} #4
H_dict = {} #5
for yr in year_list: #6
cond = ( df.index.year == yr) #7
mean = df.loc[cond, var_list].mean() #8
gdp_dict[yr] = mean.iloc[0] #9
K_dict[yr] = mean.iloc[1] #10
H_dict[yr] = mean.iloc[2] #11
コードの説明
#1：３つの変数のリスト#2：平均を計算する年のリストであり，df.index[0].yearとdf.index[-1].yearはインデックスから最初と最後の年を抽出する#3：それぞれの年のGDPの平均を格納する辞書#4：それぞれの年の資本ストックの平均を格納する辞書#5：それぞれの年の総労働時間の平均を格納する辞書#6：year_listに対してのforループ#7：df.indexの年がyrと同じ場合はTrueを返し，そうでない場合はFalseを返すSeriesをcondに割り当てる。#8：それぞれの年の行，var_listにある列を抽出し平均を計算する。.mean()は列の平均を計算する。計算した平均は
Seriesとして返されるので，それを左辺にある変数meanに割り当てる。
#9：meanの0番目の要素をgdp_dictのキーyr（年）に対応する値をして設定する。#10：meanの1番目の要素をK_dictのキーyr（年）に対応する値をして設定する。#11：meanの2番目の要素をH_dictのキーyr（年）に対応する値をして設定する。
このコードの結果を確認してみよう。例えば，
gdp_dict
{1980: 273013.55, 2021: 536811.775}
キー・値の２つの対がある。キー1980の値には1980年のGDPの平均が設定されており，同様に，キー2021の値には2021年のGDPの平均が設定されている。従って，次のコードで1980年の値にアクセスできる。
gdp_dict[1980]
他の辞書も同様である。
なぜforループを使うのか。
上のコードでforループを使っているが，もし使わなければ次のようなコードになる。
var_list = ['gdp','capital', 'total_hours']
mean1980 = df.loc['1980-03-31':'1980-12-31',var_list].mean()
mean2019 = df.loc['2019-03-31':'2019-12-31',var_list].mean()
gdp1980 = mean1980[0]
K1980 = mean1980[1]
H1980 = mean1980[2]
gdp2019 = mean2019[0]
K2019 = mean2019[1]
H2019 = mean2019[2]
同じようなコードが並んでいる。py4macroのデータがアップデートされて2020年までのデータが含まれることになったとしよう。forループを使うコードであれば2019を2020に一ヵ所だけ変更すれば良い（変更する必要もないコードを書くことも可能）。forループを使わないコードで何ヵ所修正する必要があるか数えてみよう。面倒と感じないだろうか。それに修正箇所が増えればbugが紛れ込む可能性が増えていく。
次に全要素生産性の平均四半期成長率を計算するためにdfの列名を確認しよう。
df.columns
Index(['gdp', 'consumption', 'investment', 'government', 'exports', 'imports',
'capital', 'employed', 'unemployed', 'unemployment_rate', 'hours',
'total_hours', 'inflation', 'price', 'deflator', 'gdp_cycle',
'hours_trend', 'hours_cycle', 'employed_trend', 'employed_cycle',
'total_hours_trend', 'total_hours_cycle', 'tfp_log', 'tfp_log_trend',
'tfp_cycle', 'tfp_cycle_lag'],
dtype='object')
gdpは0番目，capitalは6番目，total_hoursは11番目の列にある。この情報をもとに次のように計算しよう。
var_idx = [0,6,11] # 1
df0 = df.iloc[0,var_idx] # 2
df1 = df.iloc[-1,var_idx] # 3
# 4
B0 = df0['gdp']/( df0['capital']**a * df0['total_hours']**(1-a))
# 5
B1 = df1['gdp']/( df1['capital']**a * df1['total_hours']**(1-a))
no_quarters = len(df) # 6
g = (B1/B0)**(1/no_quarters)-1 # 7
g
0.0024043695054836167
コードの説明
3つの変数の列インデックのリスト
最初の行で3つの変数だけを抽出し
df0に割り当てる。最後の行で3つの変数だけを抽出し
df1に割り当てる。最初の四半期の全要素生産性の計算
最後の四半期の全要素生産性の計算
四半期の数
平均四半期成長率の計算
次にGDPに対する投資の割合の平均を計算する。
s = (df['investment']/df['gdp']).mean()
s
0.21777291680724095
約22％が投資に回されている。
シミュレーション#
関数#
シミュレーションのコードを関数としてまとめよう。
def stochastic_solow(T=160, # 160=40*4 40年間の四半期の数
rho=rho,
sigma=sigma,
g=g,
s=s,
d=0.025,
a=a,
H=H_dict[1980], # 1980年平均に固定
K0=K_dict[1980],
B0=B0):
"""引数：
T: シミュレーションの回数
rho: AR(1)のパラメータ
sigma: ホワイト・ノイズ分散
g: TFPの平均四半期成長率
s: 貯蓄率
d: 資本減耗率
a: 生産関数のパラメータ（資本の所得割合）
H: 総労働時間
返り値
次の変数を含むDataFrame
GDP，消費，投資，資本ストック，全要素生産性の水準，トレンド，変動"""
# ========== A,B,Kの計算 ==========
# 計算結果を格納するリスト
A_list = [1] # 1
B_list = [B0] # 2
K_list = [K0] # 3
# 初期値でありアップデート用変数
A = 1 # 4
B = B0 # 5
K = K0 # 6
# A,B,Kの時系列の計算
for t in range(1,T):
K = s * A*B * K**a *H**(1-a) + (1-d)*K # 7
e = np.random.normal(0,scale=sigma) # 8
A = A**rho * np.exp(e) # 9
B = B0*(1+g)**t # 10
A_list.append(A) # 11
B_list.append(B)
K_list.append(K)
# ========== DataFrameの作成 ==========
df_sim = pd.DataFrame({'K':K_list, # 12
'A':A_list,
'B':B_list})
# ========== Y,C,Iの計算 ========== # 13
df_sim['Y'] = df_sim['A']*df_sim['B']*df_sim['K']**a * H**(1-a)
df_sim['C'] = (1-s)*df_sim['Y']
df_sim['I'] = s*df_sim['Y']
# ========== トレンドとサイクルの計算 ==========
for v in ['K','A','Y','C','I']:
# ---------- レンドの計算 ---------- # 14
df_sim[f'{v}_log_trend'] = py4macro.trend( np.log(df_sim[v]) )
# ---------- サイクルの計算 ---------- # 15
df_sim[f'{v}の変動'] = np.log( df_sim[v] ) - df_sim[f'{v}_log_trend']
return df_sim
コードの説明
Aの値を格納する空のリスト。初期値が入っている。Bの値を格納する空のリスト。初期値が入っている。Kの値を格納する空のリスト。初期値が入っている。forループでアップデートされるAの値を一次的に格納する。最初は初期値が入っている。forループでアップデートされるBの値を一次的に格納する。最初は初期値が入っている。forループでアップデートされるKの値を一次的に格納する。最初は初期値が入っている。右辺では(4)~(6)の値（
A，B，K）を使い来期のKを計算し，その値を左辺のKに割り当てる。この瞬間に(6)のKがアップデートされる。来期のホワイト・ノイズの値を生成し
eに割り当てる。右辺では(4)の値（
A）を使い来期のAの値を計算し，その値を左辺のAに割り当てる。この瞬間に(4)のAがアップデートされる。右辺では来期の
Bの値を計算し，その値を左辺のBに割り当てる。この瞬間に(5)のBがアップデートされる。リスト(1)~(3)に来期値を追加する。
リスト(1)~(3)を使いDataFrameの作成
Y，C，Iを計算し新たな列としてDataFrameに追加トレンドを計算し新たな列としてDataFrameに追加
f-stringを使って新しく作る列名をK_trendのように_trendを追加する。
トレンドからの乖離を計算し新たな列としてDataFrameに追加
f-stringを使って新しく作る列名をK_cycleのように_cycleを追加する。
実行#
デフォルトの値でシミュレーションを行い、最初の5行を表示する。
df_sim = stochastic_solow()
df_sim.head()
| K | A | B | Y | C | I | K_log_trend | Kの変動 | A_log_trend | Aの変動 | Y_log_trend | Yの変動 | C_log_trend | Cの変動 | I_log_trend | Iの変動 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 8.500536e+05 | 1.000000 | 0.366275 | 272036.578508 | 212794.379328 | 59242.199180 | 13.692619 | -0.039565 | 0.003195 | -0.003195 | 12.531130 | -0.017438 | 12.285520 | -0.017438 | 11.006828 | -0.017438 |
| 1 | 8.880445e+05 | 1.010306 | 0.367155 | 279871.577227 | 218923.127523 | 60948.449704 | 13.725641 | -0.028864 | 0.002517 | 0.007736 | 12.544741 | -0.002655 | 12.299131 | -0.002655 | 11.020439 | -0.002655 |
| 2 | 9.267918e+05 | 1.011750 | 0.368038 | 285298.349628 | 223168.095869 | 62130.253759 | 13.758638 | -0.019154 | 0.001837 | 0.009845 | 12.558341 | 0.002949 | 12.312731 | 0.002949 | 11.034039 | 0.002949 |
| 3 | 9.657523e+05 | 1.005766 | 0.368923 | 288538.617366 | 225702.721051 | 62835.896315 | 13.791568 | -0.010905 | 0.001157 | 0.004592 | 12.571918 | 0.000666 | 12.326308 | 0.000666 | 11.047616 | 0.000666 |
| 4 | 1.004444e+06 | 0.997704 | 0.369810 | 291000.251128 | 227628.277648 | 63371.973480 | 13.824375 | -0.004430 | 0.000488 | -0.002786 | 12.585460 | -0.004381 | 12.339850 | -0.004381 | 11.061158 | -0.004381 |
\(Y\)の水準とトレンドの図示（対数）。
np.log(df_sim['Y']).plot(title='シミュレーション：Yとトレンド（対数）')
df_sim['Y_log_trend'].plot()
pass
トレンドの傾きは成長率を表しているが，徐々に減少していることが分かる。最初は成長率に対する資本蓄積の効果が大きいが，定常状態に近づくにつれて資本蓄積の効果が減少しているためである。
\(Y\)の変動（サイクル）の図示。
ax_ = ( 100 * df_sim['Yの変動'] ).plot(title='シミュレーション：Yの変動（％）')
ax_.axhline(0, c='red')
pass
絶対値で約2％内に収まった変動となっている。データでは，リーマンショック時はマイナス約4％程乖離したが，それ以外は絶対値で2％以内に収まっており，シミュレーションは概ねデータと同じような特徴を持っていると言える。
結果#
それぞれの変数の自己相関係数を計算してみる。
var_list = ['Y','C','I','K','A']
print('\n--- シミュレーション：変動の自己相関係数 ---\n')
for v in var_list:
ac = df_sim[f'{v}の変動'].autocorr()
print(f'{v}の変動： {ac:.3f}')
--- シミュレーション：変動の自己相関係数 ---
Yの変動： 0.448
Cの変動： 0.448
Iの変動： 0.448
Kの変動： 0.969
Aの変動： 0.458
まずシミュレーションを行う度にこれらの値は変化することに留意し結果を考えよう。ある程度の自己相関が発生しており，\(A\)の影響が反映されている。\(Y\)，\(C\)，\(I\)の値は同じなのは，\(C\)，\(I\)は\(Y\)の線形関数であるためであり，これらの３つの変数は\(A\)の変動と似た動きになっている。また\(K\)の値が大きいのは、ストック変数であるため変化に時間が掛かるためある。データdfに含まれる変数の自己相関係数と比べるために，まず消費，投資，資本の変動を計算する。
df.columns
Index(['gdp', 'consumption', 'investment', 'government', 'exports', 'imports',
'capital', 'employed', 'unemployed', 'unemployment_rate', 'hours',
'total_hours', 'inflation', 'price', 'deflator', 'gdp_cycle',
'hours_trend', 'hours_cycle', 'employed_trend', 'employed_cycle',
'total_hours_trend', 'total_hours_cycle', 'tfp_log', 'tfp_log_trend',
'tfp_cycle', 'tfp_cycle_lag'],
dtype='object')
data_var_list = ['consumption','investment','capital']
for v in data_var_list:
df[f'{v}_cycle'] = np.log( df[v]/py4macro.trend(df[v]) )
結果を示すために，data_var_listの最初と最後にgdpとtfpを追加する。
data_var_list = ['gdp']+data_var_list+['tfp']
print('\n--- データ：変動の自己相関係数 ---\n')
for v in data_var_list:
ac = df[f'{v}_cycle'].autocorr()
print(f'{v:>11}の変動： {ac:>5.3f}') # 1
--- データ：変動の自己相関係数 ---
gdpの変動： 0.681
consumptionの変動： 0.460
investmentの変動： 0.840
capitalの変動： 0.884
tfpの変動： 0.510
コードの説明
f-stringを使って変数vとacを代入している。:>11はvを表示する幅を11に設定し右寄せ（>）にしている。:>5はacを表示する幅を5に設定して右寄せ（>）にしている。.3fはacを小数点第三位までを表示するように指定している。
実際の資本ストックは非常に大きな持続性があり、シミュレーションの結果はそれを捉えている。しかし他の変数に関しては持続性が足りないようである。
しかしシミュレーションの結果はランダム変数の実現値に依存しているので，シミュレーションを行う度に異なる結果が出てくる。また，このシミュレーションではT=160として40年間を考えたが，Tの値によっても結果は大きく左右される。Tを大きくするとより安定的な結果になっていくだろう。例えば，T=10_000を試してみよう。
次に、\(Y\)との相関係数を計算する。
print('\n--- シミュレーション：GDPとの相関係数 ---\n')
for v in var_list[1:]:
cov = df_sim[['Yの変動',f'{v}の変動']].corr().iloc[0,1]
print(f'{v}の変動： {cov:.3f}')
--- シミュレーション：GDPとの相関係数 ---
Cの変動： 1.000
Iの変動： 1.000
Kの変動： 0.074
Aの変動： 0.979
\(C\)と\(I\)の値は1.0となっている理由は\(Y\)と線形になっているためであり，全く同じように変動するためである。
\(Y\)と\(A\)は非常に強く相関している。生産関数を見ると，両変数も線形の関係があり，\(K\)の変化によって変化のズレが生じるメカニズムとなっているためである。
相関度をデータと比べてみよう。
print('\n--- データ：GDPとの相関係数 ---\n')
for v in data_var_list[1:]:
cor = df[['gdp_cycle',f'{v}_cycle']].corr().iloc[0,1]
print(f'{v:>11}の変動：{cor:>5.3f}')
--- データ：GDPとの相関係数 ---
consumptionの変動：0.766
investmentの変動：0.769
capitalの変動：0.154
tfpの変動：0.897
\(C\)と\(I\)に関してはシミュレーションの相関係数は大きすぎる結果となっている。
確率的ソロー・モデルによるシミュレーション結果は，マクロ変数の特徴を捉えている部分もあるが景気循環モデルとして考えるには難しいだろう。消費者の効用最大化問題が欠落しているため，\(Y\)，\(C\)，\(I\)が全く同じように動くのがネックになり，消費や投資の特徴を捉えることができていない。後の章で議論するが，消費者の最適化問題や将来の予測などを導入したRBCモデルは，よりデータに近い結果となる。しかし，いくつか問題が残るのも事実である。
TFPショックとは何なのか？
生産関数を見る限りTFPは生産性を表している。技術水準，政府の規制，政治システムや経済制度の変化を捉えると考えられ，例えば，ビジネスに対する規制の強化によって負のショックが発生する。しかし，それだけで実際に観測されるGDPの大きな変動（例えば、リーマン・ショックやコロナ・ショック）を説明するのには疑問である。
ソロー・モデルも含めて価格調整が瞬時に行われるが，次の問題が残る。
古典派の二分法が成立すると仮定するが、短期的には無理がある仮定（例えば、賃金）。
景気循環に対して金融政策は無力なのか？
このような問題意識を持って，RBCモデルとニューケインジアン・モデルを議論していくことにする。