マクロ経済学の動的計画法説明用コードのMatlabからRへの移植

f:id:uncorrelated:20161224234007p:plain

一橋大学の阿部教授が院生向けのレクチャー・ノートをアップロードしていたのを見つけ、拝読して勉強させて頂いていたのですが、説明用のコード（pp.16—20）がMatlabだったのでRに移植してみました。プロプライエタリなMatlab嫌いな人に役立つかも知れないので、公開します。

単に移植するだけではなく、ループ周りを中心に無駄を省いています。理由は、行列にしておいた方が変数の中身を説明しやすいのと、コード全体の見通しが良くなるのと、速度面で有利になるからです。資本量の区切り幅inK=0.4ぐらいにして、行列の中身を表示させると、何を演算しているのか理解しやすいと思います。

alpha <- 0.25 # production parameter
beta <- 0.9 # subjective discount factor
gamma <- -2 # preference parameter
delta <- 1 # 1 -depreciation rate
minK <- 0.2 # minimum value of the capital grid
maxK <- 1.8 # maximum value of the capital grid
inK <- 0.001 # size of capital grid increments
nK <- ceiling((maxK-minK)/inK + 1) # number of grid points
# 横（列）の当期資本と、縦（行）の来期資本で定まる効用を表す行列
U <- matrix(rep(-Inf, nK^2), nK, nK)
K <- t(replicate(nK, ((1:nK)-1)*inK + minK)) # 当期資本。nK×nKに拡張。当期資本は列ごとに変えているのに注意。
K_prime <- replicate(nK, ((1:nK)-1)*inK + minK) # 来期資本k_{t+1}。nK×nKに拡張。来期資本は行ごとに変えているのに注意。
I <- K_prime - delta*K # 投資量。横（列）の当期資本と、縦（行）の来期資本で定まる。
C <- ((1-beta)/(alpha*beta))*(K^(alpha)) - I # (79)式を(80)式に代入した式から消費量を計算
U[C>0] <- ((C^(gamma+1))/(gamma+1))[C>0] # 効用を計算して行列に格納。消費がゼロ以下の場合は、初期値-Infを保存している
#
# whileループのための変数を初期化
#
V <- rep(0, nK) # value functionの代わりになるベクトル
Decis <- rep(0, nK) # Vに対応する来期資本を表すグリッド位置
iter <- 0 # ループ回数
metric <- 10 # 更新前後のvalue functionの差分の最大要素の値が入る。初期値は終了条件より大きければ、何でも良い。
while(metric > 0.001){
  # (82)式の最大化問題に対応する来期資本を探す
  # beta*Vは1列分のベクターだが、自動的にnK列分に複製拡大して足されるのに注意
  r <- U + beta*V
  # それぞれの列の最大値を選んで、value functionとする。
  # 当期資本に対応する列の中の最大値は、最大化問題の解となる来期資本
  tmpV <- apply(r, 2, max) # value functionの値
  tmpDecis <- apply(r, 2, which.max) # value functionに対応する来期資本を表すグリッド位置
  metric <- max(abs(V - tmpV)) # 終了条件と比較する更新によって生じた差分を計算
  V <- tmpV;
  Decis <- tmpDecis;
  iter <- iter + 1 # ループ回数をカウント
# 計算状況を表示
  vfor8 <- V[1] # value of V for k=kmin=0.2
  ufor8 <- (Decis[1] - 1)*inK + minK # value of control for k=0.2
  print(sprintf("ループ回数:%d, VF変化量:%.4f, vfor8:%.4f, ufor8:%.4f", iter, metric, vfor8, ufor8))
}
# policy functionと対応する各期の効用を計算
policy <- (Decis - 1)*inK + minK # policy function
U_op <- rep(-Inf, nK) # Utility under the optimal policy
K <- ((1:nK)-1)*inK + minK
I <- policy - delta*K;
C <- (1-beta)/(alpha*beta)*(K^alpha) - I
U_op[C>0] <- ((C^(gamma+1))/(gamma+1))[C>0]
# 各期の効用の割引価値を計算
betam <- beta*rep(1, nK)
value <- U_op/(rep(1, nK)-betam)
# 結果をプロット
x_at <- seq(1, length(value), length.out=17)
x_labels <- sprintf("%.2f", seq(minK, maxK, length.out=length(x_at)))
y_unit <- 5
ylim <- c(floor(min(value)/y_unit)*y_unit, ceiling(max(value)/y_unit)*y_unit)
y_at <- sprintf("%.2f", seq(ylim[1], ylim[2], length.out=5))
plot(value, main="DETERMINISTIC GROWTH MODEL: VALUE FUNCTION", xlab="Capital", ylab="Present Value", axes=FALSE, ylim=ylim)
axis(1, at=x_at, labels=x_labels, las=2)
axis(2, at=y_at)