PythonでReal Business Cycle Model その1

RBCモデルは新古典派系動学的一般均衡モデル

さて、まずRBCモデルについての説明から行う必要がありますね。RBCモデルとは、マクロ経済学の中でも景気循環分析(Business Cycle Analysis)で用いられる新古典派系動学的一般均衡モデルで、景気変動の源泉を実物要因、特に財の生産技術に対する確率的なショックに求めている点が特徴です。日本語訳は実物景気循環モデルで、その名からわかるように価格は伸縮的で貨幣中立が成り立っており、名目変数が景気循環に影響しないモデルとなっています。

Kydland and Prescott(1982)がその嚆矢

RBCモデルは以下の論文で公開されました。

Kydland, F. and E. Prescott (1982), “Time to Build and Aggregate Fluctuations,” Econometrica 50, 1345-1371.

それまでの動学的一般均衡モデルとしては、ラムゼイモデルが有名でした。RBCモデルはこのラムゼイモデルで外生として扱われていた労働供給を内生化し、また生産関数に対する技術的なショックを導入しています。のちにKydlandとPrescottはこの業績からノーベル経済学賞を受賞しています。

現在ではNew Keynesian モデルが主流だが、議論の発射台として重用されている

少なくとも実務の世界ではRBCモデルを使用している機関は存在しないと思います。RBCが仮定している価格の伸縮性が現実的な仮定ではないからです。現在は価格の硬直性を導入したNew Keynesianモデルが広く使用されています¹。では、RBCモデルは有用ではないのかというと、そういうわけではありません。New KeynesianモデルはRBCモデルに比べて、モデルが複雑になるため、モデルを拡張した際にその拡張が経済的にどのような意味を持っているかを解釈するのが大変です。そのため、モデル拡張の初期においてはプレーンなRBCモデルを用いて研究を行い、その特徴を調べ、その後実際のデータにフィッテイングする際にNew Keynesianモデルが使用されるケースがあります。よって、現在においてもRBCモデルは議論の発射台として有用です。

モデル概観

ここで、簡単にモデルを概観しておきましょう。経済主体として、無期限間²生きる代表的家計³と企業が存在します。家計は企業に労働と資本⁴を供給し、企業はそれらを用いて財を生産します。家計は企業から受け取る賃金と資本のレンタル料を、消費と貯蓄に振り分け、財の消費から効用を得ます。この経済には、資本と国債の2つの資産が存在し、家計はこれらを購入することで貯蓄を行うことができます。 ここまでの議論から、経済には以下4つの市場が存在します。

1. 財市場
2. 労働市場
3. 資本市場

これらの市場は完全競争市場であると仮定します。家計は生涯効用の最大化を、企業は利潤最大化を達成するよう最適な行動を行うと仮定します。つまり、以下のような最適化問題に議論を帰着させることができます。

Matlab(Octave) + Dynareがデファクトスタンダード

動学的確率的一般均衡モデルにはDynareと呼ばれる専用のプログラミング言語が存在します。このDynareは単体で動くのではなく、Matlab(Octave)上で動かす必要がありますが、インパルス応答関数等の数値シミュレーションに加え、状態空間モデルの推定や構造パラメータのMCMC推定など動学的一般均衡モデルの基本的な分析を行うには十分な機能を備えています。ただ、先進的な研究を行う場合、特に推定方法の高度化、Matlab単体を使用して分析が行われることが多いと感じます。Dynareは手を動かして既存のモデルを動かしたい人に向けたソフトウェアであるといってもよいかもしれません(もちろんDynareでも論文は書けますし、実際に存在もします)。

PythonでRBC Modelを分析する意義

Matlab + Dynareがデファクトスタンダードであるのに、Pythonを使用する理由は以下3点です。

1. モデル分析以外の部分の利便性が高い
   • このブログでも紹介したようにPythonにはpandasdmxやpandas_datareader等APIを使用してマクロ経済データを直接取得できる便利な関数群があります。これまでの分析では、実際にWEBサイトへ行き、手でcsvファイルをダウンロードして、それをExcelで加工し分析用のデータを作成していました。
2. より高度な分析のための外部ライブラリの統合が容易
   • 動学的一般均衡モデルの実装に際して、状態空間モデルであったりMCMCに関する実装を一から行う気は毛頭ありません。statsmodelsやPymc3といった外部ライブラリを利用することを考えています。また、粒子フィルタなどさ らに発展的な技術についてもPythonの方が実装済みコードが存在しており、かつClassという概念を用いればこれ らの特性が容易に継承可能である点からもPythonでの実装がBestであると考えました。
3. みんなPythonを使っているから
   • せっかくであれば、自分だけじゃなくいろんな人に使ってほしいので、利用者が多いPythonを選びました。

このPostでやり遂げたいこと

各モデル方程式をSympyで定義し、それを対数線形化、その後Canonical formと呼ばれる行列形式に整理するまでをクラス化したいと思います。

Canonical formとは

おそらく聞き慣れない単語であろうCanonical formについて説明しておきます。先ほど導出したRBCモデルの均衡条件式は非線形差分方程式システムとなっています。動学的確率的一般均衡モデルでは、これを定常値周りで対数線形近似して線形化するのが一般的です。線形化されたシステムを行列形式で表現すると以下のようになります。

\[ Bs_{t+1} = Cs_{t} + \Psi\epsilon_t \] ここで\(s_t\)は内生変数で\(\epsilon_t\)は外生変数、\(B,C,\Psi\)は係数行列です。このシステムをCanonical formと呼びます。今回は均衡条件式をこちらの形式へ整形することがゴールです。

利用するモジュール群のインポート

以下を使用します。

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
import sympy as sym

Sympyでの方程式の定義

先ほど導出した方程式をSympyで定義します。 まず内生変数を宣言します。

sym.init_printing()
var = sym.var('l, c, k, l, z, r, LEAD_c, LEAD_k, LEAD_l, LEAD_r, LEAD_z, LEAD_R')
endog = sym.var('l, c, k, z')
LEADs = sym.var('LEAD_l, LEAD_c, LEAD_k, LEAD_z')

次にパラメータを宣言します。

params = sym.var('beta, theta, varphi, alpha, delta, phi, Psi')

方程式を定義します。

# Eliminate Price  
LEAD_R  = (LEAD_z*alpha*(LEAD_k/LEAD_l)**(alpha-1))
w       = (1-alpha)*z*(k/l)**alpha
  
# Optimal Conditions  & state transition
labor   = l**varphi-w/(Psi*c**(theta))
euler   = c**(-theta) -(LEAD_c**(-theta))*beta*(1+LEAD_R-delta)
capital = LEAD_k - (1-delta)*k - z*(k**alpha)*(l**(1-alpha)) + c
tech    = LEAD_z - phi*z
  
optcon  = sym.Matrix([labor, euler, capital, tech])

Sympyの自動数式処理機能を用いた自動対数線形近似の実施

optconに保存されているのは非線形差分方程式システムとなっています。次に対数線形近似を行い、線形システムへ変換を行います。対数線形近似はそれ自体に特別な意味があるとは思えませんので、Sympyのjacobianを用いてヤコビ行列を計算し、自動処理します。

# Differentiation 
jopt    = optcon.jacobian(endog).subs([(LEAD,endog[i]) for i, LEAD in enumerate(LEADs)])
jopt_LEAD = optcon.jacobian(LEADs).subs([(LEAD,endog[i]) for i, LEAD in enumerate(LEADs)])

ここでは、内生変数\(s_t\)と\(s_{t+1}\)を分けて処理しています。

対数線形近似後の係数を得るためには、ヤコビ行列のそれぞれの要素に、対応する変数の定常均衡値を掛け合わせる必要があります。ヤコビ行列の各要素は\(\partial f(x_t)/\partial x_t\)ですが、対数線形近似=関数\(f(x_t)\)を定常値からの乖離\(x_t/x\)で(定常値周りで)一次近似することですので、定常値\(x\)を掛け合わせることで係数\(\partial f(x_t)/(\partial x_t/x)\)を求めることができます。
よって、定常値の計算が必要です。数値的に求めることもできますが、定常均衡はそれを求めること自体もモデルの特性を知る上で重要ですので、解析的に解いた数式をSympy expressionとして定義します。

# Steady State Equation needed to be calculated by hand
kls     = (((1/beta)+delta-1)/alpha)**(1/(alpha-1))
wstar   = (1-alpha)*(kls)**alpha
clstar  = kls**alpha - delta*kls
lstar   = ((wstar/Psi)*(clstar**(-theta)))**(1/(varphi+theta))
kstar   = kls*lstar
cstar   = clstar*lstar
zstar   = 1
Ystar   = (kstar**alpha)*(lstar**(1-alpha))
ss_eq   = sym.Matrix([lstar, cstar, kstar, zstar])

次に、この定常均衡方程式システムを構造パラメータを引数とするPython関数へ変換します。MCMCで構造パラメータの推定を行う際などは頻繁にパラメータを更新しますので、関数化しておくことは便利です。また、構造パラメータと定常値を引数として対数線形システムの係数行列を評価する関数もPython関数へ変換しておきます。

# Translating Sympy expressions into Python functions for steady state and coefficient matrix
fss     = sym.lambdify([params],ss_eq)
fcoef   = sym.lambdify([params,endog],jopt)
fcoef_LEAD = sym.lambdify([params,endog],jopt_LEAD)

これで解析的には準備が整いました。では、構造パラメータに数値を代入して、係数を計算してみます。

# Evaluate steady state and each derivative in terms of % deviations from ss
vparams = np.array([0.99, 1.5 ,2 , 0.3, 0.025, 0.8, 1])
ss = fss(vparams)
B = np.matrix(fcoef_LEAD(vparams,ss)*ss, dtype='float')
C = np.matrix(fcoef(vparams,ss)*ss, dtype='float')
A = np.dot(np.linalg.inv(np.matrix(C)),B)

print(A)

## [[-1.02747592e-02  3.21263317e-01 -5.48161413e-03 -4.68593779e-01]
##  [ 3.19823603e-02 -1.00000000e+00 -1.49198793e-03  4.01094076e-02]
##  [ 4.95211793e-02 -1.54839039e+00 -9.81949430e-01  3.57889915e+00]
##  [ 0.00000000e+00  0.00000000e+00  0.00000000e+00 -1.25000000e+00]]