動学マクロ経済学に関する蓮見氏のレクチャーノートを拝読していて、もう少しRのソースコード部分に補足説明があったほうが分かりやすい気がしたので、ラムゼー・モデル（Ramsey–Cass–Koopmans model）に関して紹介してみたいと思います。
このモデルは連続関数を仮定して変分法で解くこともできるのですが、最近は離散問題に置き換えてコンピューターで数値解を求める方向*1で教えているようです。これは、RBC/DSGEへの拡張が容易である*2ためでしょう。
何はともあれ、こういうモデルを計算してみると、マクロ経済学を直感的に理解できるようになると思います。解き方としてはモデルの式を整理して非線形連立方程式の形に直した上で、始点と終点を整理し、ニュートン・ラフソン法などで解くだけなので、実際にやってみましょう。
1. モデル式 → 非線形連立方程式
まずはモデルを確認しておきましょう。詳しくは蓮見氏のレクチャーノートを見てください。
1.1. 家計
元々のモデルではあまり特定化された関数は置いていない気がしますが、家計の目的関数になる効用関数を以下のようにおきます。
期の消費が
になるわけですが、割引率
があるので、将来の事はあまり気にしません。
制約条件は、以下のようになります。
資本に利率
を乗じたものが利子収入、
が賃金収入です。つまり、消費
と投資
は収入以下になると言う制約です。
来期の資本の量は今期の資本と投資の合計から、減価償却を引いたものになります。
これを制約条件に代入しておきましょう。
さて、制約付最大化問題を解くために、ラグランジュアンを置きます。
家計が決定できるのは、消費量と投資量
なので、それぞれで偏微分をして一階条件を出します。
の偏微分も忘れずに。
1.2. 企業
企業の生産関数は以下のように、コブ＝ダグラス型として起きます。
人口が増えないとしてと置きます。完全競争を仮定して、利子率と賃金を求めます。利子率は資本の限界生産性、賃金は労働の限界生産物に等しくなるわけですね。
生産物は利払いと賃金に分割されるので、以下のように企業の生産物と家計の収入が等しいくなります。
2. カリブレーション
ここからが本題で、蓮見氏のレクチャーノートのソースコードにコメントが無いので、勝手に補足していきます。
数値計算する式は上の二式ですが、パラメーターや初期値を定めないといけませんし、非線形連立方程式を解いてくれるパッケージの使い方も確認する必要があります。
2.1. パラメーター
と
以外の変数はパラメーターになります。任意でいいのですが、
、
、
、
となっていました。
2.2. 初期値
マクロ経済学と言うか数値解析の問題ですが、と
に初期値が必要です。
が1から30まで計算するとすると、2×30=60変数を初期化する必要があります。レクチャーノートでは何故か定常点
が入っているのですが、0より大きな値であれば何でも問題ないです。
2.3. 計算までのソースコード
ようやくコーディングできるだけの情報が揃ったので計算してみましょう。なおレクチャーノートでは消費と投資の行列を作って崩していたのですが、ベクトルを分割するようにしています。
#
# 非線形連立方程式を解くパッケージを読み込む
#
library(nleqslv)
#
# パラメーターの設定
#
alpha <- 0.3
beta <- 0.99
delta <- 0.25
At <- 1.0
#
# 定常状態の計算
# （Kは値で、Cは作図のために関数で戻す）
#
steadystate <- function(){
sK <- ((1/beta+delta-1)/alpha/At)^(1/(alpha-1))
list(K=sK, C=function(K=sK){
At*K^alpha - delta*K
})
}
SS <- steadystate()
#
# カリブレーションを行う関数
# maxT:計算する期間の長さ, K01:資本ストック初期値
#
calibrate <- function(maxT, K01){
# レキシカルスコーピングでcalibrate内で定義された変数maxTとK01は、
# parseXとobjfunからアクセスできる一方でグローバル変数にならない
#
# 計算結果から消費と資本ストックを分ける
# X:計算結果
#
parseX <- function(X) {
# Xは(C[1], C[2], ... , C[maxT], K[2], K[3], ... , K[maxT + 1])と言うベクトル
# Xの前半部分が各期の消費
# 最後に定常状態を付け加える
C <- c(X[1:maxT], SS$C())
# Xの後半部分が各期の投資
# 最初に初期値を付け加える
K <- c(K01, X[(maxT+1):(2*maxT)])
# リスト構造で値を戻す
list(C=C, K=K)
}
#
# 右辺ゼロの非線形連立方程式を表す関数
# X:各期の消費と資本ストックを表すベクトル
#
objfun <- function(X) {
# 消費と資本ストックに分ける
X <- parseX(X)
# 戻り値のために空ベクトルを創る
rC <- numeric(maxT)
rK <- numeric(maxT)
# 線形連立方程式に値を代入していく
# 2*maxT本の連立方程式になる事に注意
for(t in 1:maxT) {
rC[t] <- X$C[t+1]/X$C[t] - beta*(alpha*At*X$K[t+1]^(alpha-1)-delta+1)
rK[t] <- X$K[t+1] - At*X$K[t]^alpha - (1-delta)*X$K[t]+X$C[t]
}
# 戻り値を返す
return(c(rC, rK))
}
# 初期値は基本的に任意の値で構わないが、連立方程式内に割り算があるため0は避ける
Xinit <- c(rep(1, maxT), rep(1, maxT))
# 非線形連立方程式を解く
r <- nleqslv(Xinit, objfun)
# 結果から消費と資本ストックを分ける
parseX(r$x)
}
# 期間の長さ
maxT <- 30
# 初期資本ストックは定常状態の半分として計算
X1 <- calibrate(maxT, 0.5*SS$K)
# 初期資本ストックは定常状態の倍として計算
X2 <- calibrate(maxT, 2*SS$K)
講義などで発散する経路が必要な場合は、以下のように計算します。
#
# 非定常経路の計算
#
inertia <- function(maxT, K01, C02){
# 戻り値のために空ベクトルを作る
C <- numeric(maxT)
K <- numeric(maxT)
K[1] <- K01
C[1] <- C02
for(t in 2:maxT){
# 資本の増加式
K[t] <- At*K[t-1]^alpha + (1 - delta)*K[t-1] - C[t-1]
# オイラー方程式に従って来期の消費を決定
C[t] <- beta*(alpha*( K[t-1]^alpha + (1 - delta)*K[t-1] - C[t-1])^(alpha-1) - delta + 1)*C[t-1]
}
list(C=C, K=K)
}# 消費過少経路を計算
NSUC <- inertia(maxT, 0.5*SS$K, X1$C[1] - 0.02)# 消費過剰経路を計算
NSOC <- inertia(maxT, 0.5*SS$K, X1$C[1] + 0.02)
2.4. 計算結果を描画する
X1$CとX1$Kを眺めてもいいのですが、見栄えが悪いので描画してみましょう。
#
# 定常状態を描画
#
curve(SS$C(x), 0, 3, ylim = c(0.5, 1.1), xlab = expression(K[t]), ylab = expression(C[t]), main="ラムゼーモデルの位相図", lty="dotted")
lines(cbind(c(SS$K, SS$K),c(0, 100)), lty="dotted")
points(SS$K, SS$C(), pch=19)
text(2.5, 0.7, expression(paste(Delta,"K = 0",sep="")), pos=4)
text(1.20, 1.0, expression(paste(Delta,"C = 0",sep="")), pos=4)
#
# 遷移を矢印で描画
# 矢印が密集していると見づらいので、間引く
#
step <- 3
for(i in seq(1, 12, step)){
arrows(X1$K[i], X1$C[i], X1$K[i+step], X1$C[i+step], length = 0.1)
arrows(X2$K[i], X2$C[i], X2$K[i+step], X2$C[i+step], length = 0.1)
}
for(i in seq(13, maxT - 1, step)){
segments(X1$K[i], X1$C[i], X1$K[i+step], X1$C[i+step])
segments(X2$K[i], X2$C[i], X2$K[i+step], X2$C[i+step])
}
描画される図形は以下のような感じになります。
発散経路を計算した場合は、以下のように続けて描画してください。
#
# 発散経路を描画
# すぐに発散するので間引かない
#
step <- 1
for(i in seq(1, maxT, step)){
arrows(NSOC$K[i], NSOC$C[i], NSOC$K[i+step], NSOC$C[i+step], length = 0.1, lty=1, col="red")
arrows(NSUC$K[i], NSUC$C[i], NSUC$K[i+step], NSUC$C[i+step], length = 0.1, lty=1, col="blue")
}
