餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

動学マクロ経済学と言う名の非線形連立方程式を解いてみる

動学マクロ経済学に関する蓮見氏のレクチャーノートを拝読していて、もう少しRのソースコード部分に補足説明があったほうが分かりやすい気がしたので、ラムゼー・モデル(Ramsey–Cass–Koopmans model)に関して紹介してみたいと思います。

このモデルは連続関数を仮定して変分法で解くこともできるのですが、最近は離散問題に置き換えてコンピューターで数値解を求める方向*1で教えているようです。これは、RBC/DSGEへの拡張が容易である*2ためでしょう。

何はともあれ、こういうモデルを計算してみると、マクロ経済学を直感的に理解できるようになると思います。解き方としてはモデルの式を整理して非線形連立方程式の形に直した上で、始点と終点を整理し、ニュートン・ラフソン法などで解くだけなので、実際にやってみましょう。

1. モデル式 → 非線形連立方程式

まずはモデルを確認しておきましょう。詳しくは蓮見氏のレクチャーノートを見てください。

1.1. 家計

元々のモデルではあまり特定化された関数は置いていない気がしますが、家計の目的関数になる効用関数を以下のようにおきます。

\max_{C} U_t = \sum^{\infty}_{i=t} \beta^{i-t} \ln(C_{i})

t期の消費がC_tになるわけですが、割引率\beta<1があるので、将来の事はあまり気にしません。
制約条件は、以下のようになります。

\mbox{s.t. } r_t K_i + w_i - C_i - I_{i} \geq 0

資本K_iに利率r_tを乗じたものが利子収入、w_iが賃金収入です。つまり、消費C_iと投資I_{i}は収入以下になると言う制約です。
来期の資本の量は今期の資本と投資の合計から、減価償却\delta K_iを引いたものになります。

K_{i+1}=K_i - \delta K_i + I_i

これを制約条件に代入しておきましょう。

\mbox{s.t. } r_tK_i + w_i + K_i - \delta K_i - C_i - K_{i+1} \geq 0

さて、制約付最大化問題を解くために、ラグランジュアンを置きます。

\mathit{L}=\sum^{\infty}_{i=t} \beta^{i-t} \bigg[  \ln(C_{i}) + \lambda_i\{ r_tK_i + w_i + K_i - \delta K_i - C_i - K_{i+1} \} \bigg]

家計が決定できるのは、消費量C_iと投資量K_iなので、それぞれで偏微分をして一階条件を出します。\lambda_i偏微分も忘れずに。

\frac{1}{C_i}-\lambda_i=0

\beta(r_{t+1}+1-\delta)\lambda_{i+1} - \lambda_i=0

r_tK_i + w_i + K_i - \delta K_i - C_i - K_{i+1}=0

1.2. 企業

企業の生産関数は以下のように、コブ=ダグラス型として起きます。

Y_t=A_t K^\alpha L^{1-\alpha}

人口が増えないとしてL=1と置きます。完全競争を仮定して、利子率と賃金を求めます。利子率は資本の限界生産性、賃金は労働の限界生産物に等しくなるわけですね。

r_t=\alpha A_t K_t^{\alpha-1}

w_t= A_t K_t^\alpha-r_tK_t = (1-\alpha) A_t K_t^{\alpha}

生産物は利払いと賃金に分割されるので、以下のように企業の生産物と家計の収入が等しいくなります。

Y_t=A_t K_t^\alpha=r_tK_i + w_i

1.3. 非線形連立方程式

家計の効用最大化の一階条件の三つの式を整理しましょう。\lambda_ir_tw_tを消して、方程式の形式にしています。

\frac{C_{t+1}}{C_t}-\beta(\alpha A_{t+1} K_{t+1}^{\alpha-1}-\delta+1)=0

K_{t+1}-(A_t K_t^\alpha + (1-\delta)K_t - C_t)=0

1.4. 端点の整理

実は最適経路を求める問題なので、始点と終点が無いと計算できません。始点は初期保有ストック量、終点は定常状態となります。
定常状態はt期とt+1期の消費と資本ストックが同一と言う事になります。つまり、C_{t+1}=C_tK_{t+1}=K_tとなる点です。定常点を(C^*, K^*)と置くと、非線形連立方程式の二式から以下のように整理できます。

K^*=\bigg [ \frac{\frac{1}{\beta} + \delta - 1}{\alpha A_t} \bigg ]^{\frac{1}{\alpha - 1}}
C^*= A_t \cdot  K^{*\alpha}  - \delta \cdot K^*

2. カリブレーション

ここからが本題で、蓮見氏のレクチャーノートのソースコードにコメントが無いので、勝手に補足していきます。
数値計算する式は上の二式ですが、パラメーターや初期値を定めないといけませんし、非線形連立方程式を解いてくれるパッケージの使い方も確認する必要があります。

2.1. パラメーター

C_tK_t以外の変数はパラメーターになります。任意でいいのですが、\alpha=0.3\beta=0.99\delta=0.25A_t=1となっていました。

2.2. 初期値

マクロ経済学と言うか数値解析の問題ですが、C_tK_tに初期値が必要です。tが1から30まで計算するとすると、2×30=60変数を初期化する必要があります。レクチャーノートでは何故か定常点(C^*, K^*)が入っているのですが、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")
}


3. ラムゼー・モデルから分かること

最近の動学マクロ経済学は、もっと複雑なモデルが組まれていますし、パラメーターも実データから推定する世界になっていますが、このラムゼー・モデルが一つの出発的になっているのは間違いないと思います。動学マクロ経済学が、家計が将来を見通して意思決定を行う事を前提としていることが分かります。

*1:カリブレーションと言います。

*2:学生の数学力が低くて済むと言う利点もあるとは思います。