[R] ふだんと少し違うソローモデル

概要

普段とは趣向を変えて, R でソローモデルのシミュレーションをする.
そのままだとつまらないので, ソローモデルの人口成長の前提条件をロジスティック法則としてみる.
分量はPDF換算で 6 ページ.

ロジスティック法則について

ソローモデルや, 他の多くの経済モデルでは, 人口は一定に成長する「マルサス法則」に従っていると仮定することが多い. マルサス法則のもとでは, 一定の人口成長率 $r$ *1に対して, $t$ 時点での人口 $N(t)$ を

$\begin{align} \dot{N}(t)= & rN(t)\end{align}$

とすることで記述される. 上の点は $t$ での微分を表す: $\dot{N}(t):=\partial N(t)/dt$ . これば単純な線型微分方程式で, $N(t)=N(0)\exp(nt)$ が導ける. つまり, 人口成長率が一定 (かつ1より大きい) なら, 人口は指数関数的に増加するため, 未来の食料不足はこのままでは不可避であるということを主張したのがマルサスの人口論である.

しかし, 現在の状況を確認すると, 人口成長は一定ではない. 日本や欧米の先進国は少子化のため将来の人口減少が確実となっているし, 途上国であっても多くの国で出生率の低下が確認されている.

ここで, Verhulst*2 によって提案されたロジスティック法則を紹介したい. Verhurst は, 生物の種の個体数の増加率が一定でなく, 個体数の増加にともなって増加率が鈍ることに着目し, 次のような個体数 (人口) 成長モデル, ロジスティック方程式を考案した. あらたにパラメータ $s>0$ を用いて,

$\begin{align} \dot{N}(t)= & rN(t)-sN(t)^{2}\end{align}$

と表す.

つまり, マルサスの方程式の右辺は人口の増加要素だけだが, ロジスティック方程式では $-sN(t)^{2}$ という, 人口増加に対するブレーキとなる項が追加されている. この微分方程式も容易に解くことができ*3て,

$\begin{align} N(t)= & \frac{rN(0)\mathrm{e}^{rt}}{r+sN(0)(\mathrm{e}^{rt}-1)}\end{align}$

となる. 既視感を感じる方もいるかもしれないが, これがいわゆるロジスティック関数で, 当初の用法とは外れてしまったが, ロジスティック回帰やロジットリンク関数などとして, 統計学・機械学習で広範に応用されている. ロジスティック回帰の「ロジスティック」もここに由来する. そもそもなぜロジスティックと呼ぶかというと, Verhulst が兵站 (logistics) を専門にしていたから, と言われている.

ロジスティック関数はゆるやかな S字を描くような曲線となるため, 人口がロジスティック法則にしたがうなら, ある水準に達すると急速に成長が鈍化し, 頭打ちになる. 別の表現として,

$\begin{align} \dot{N}(t)= & rN(t)(1-\frac{s}{r}N(t))\\ = & rN(t)\left(\frac{r/s-N(t)}{r/s}\right)\\ = & rN(t)\left(\frac{C-N(t)}{C}\right)\end{align}$

というものがある. $N$ が $C=s/r$ の値に近づくほど, カッコ内の数字はゼロに近づくから, 一定水準に近づくほど増加率が鈍るというのがより直観的にわかる. $C$ を用いて陽に表すなら,

$\begin{align} N(t)= & \frac{N(0)\mathrm{e}^{rt}}{C+N(0)\mathrm{e}^{rt}}\end{align}$

となる. この $C$ は「環境収容力」と呼ばれ, ある環境がどれだけの数の生物を養う「収容力」があるかを表すパラメータと解釈できる. 一応, 適当にパラメータを決めたものを, マルサス法則の場合と比較して図示しておく.

f:id:ill-identified:20150712225301p:plain

f:id:ill-identified:20150712225307p:plain

f:id:ill-identified:20150712225309p:plain

マルサス法則は見ての通り指数関数的な増加をするのに対し, ロジスティックはパラメータによって形状が変わる. 「停滞期」は現在の人口100に対して収容力が小さい場合で, 「成長期」は収容力が十分大きい場合である. 以下の図のように, ロジスティック法則では人口が収容力を超えることがない. 始めはマルサス法則のように指数関数的に増加するが, 収容の限界に近づくほど増加が鈍る.

ソローモデルのシミュレーション

では, このロジスティック法則をソローモデルにあてはめるとどうなるか. Romer (2012, Advanced Macroeconomics) での表現に沿ってソローモデルを記述する*4と, 1人あたりの資本蓄積 $k(t)=K(t)/A(t)L(t)$ は

$\begin{align} \dot{k}(t)= & sf(k(t))-(r+g+\delta)k(t)\end{align}$

と表せる. $r$ が人口成長率*5, $g$ は技術進歩率で, 技術進歩もマルサス法則のように一定とする. $\delta$ は減価償却率で, $0<\delta<1$ とする. $f$ は労働効率あたりの生産関数だが, とりあえずコブ・ダグラス型とする. 第1項が $k$ の直線で, 第2項が曲線だから, $\dot{k}=0$ となるとき, つまり $k$ の不動点はこの2つの線の交点で表される.

しかし, 今回はロジスティック法則を考えたいので, $r$ は定数ではなく, $t$ の関数として表される人口 (労働人口) の変化率で, いわば $\dot{N}(t)$ である. そのため, 明らかに第1項は直線ではなく, より複雑な曲線になる. 実は, ロジスティック法則を適用したソローモデルについて, Ritelli, Mingari Scarpello, & Brida (2008, The Solow model with logistic manpower: a stability analysis) という発表がすでに存在する. この論文では, ちゃんと解析的に解いているが, 今回は R の練習も兼ねて, 変数 $k$ や $y$ がどう変化するかを見るために, R でグラフを作成する. プログラミングでは連続的な変化を扱えないため, 微分を差分で近似して表現する: 微小時間 $\Delta$ が経過した前後の $k$ の変化量は,

$\begin{align} k(t+\Delta)-k(t)= & \left\{ sf(k(t))-(r+g+\delta)k(t)\right\} \Delta\end{align}$

だから, 1期ごとの変化を

$\begin{align} k(t+1)= & k(t)+sf(k(t)-(r+g+\delta)k(t))\end{align}$

と表せる. あとは初期値を与え, 繰り返し $r$ の毎期の人口変化率を計算して代入すれば $k(t)$ の推移を計算できる. なお, 二神・堀 (2009, マクロ経済学) では最初から差分方程式を用いてソローモデルを説明しており, その結果と若干式が異なるが, 結果に大きな違いはない.

結果

プログラムは最後に掲載した. 結果は ts オブジェクトで出力されるが, 見栄えが悪いので ggplot2 で加工した. 人口成長の曲線はすでに示した通りで, 技術進歩も一定であるのであえて図示する必要はない. あとは $k$ の動きがどうなるかが焦点となるので, $k$ と $y$ だけを図示する.

f:id:ill-identified:20150712224831p:plain

この通り, 人口成長期であっても, 停滞期であっても, $k$ が収束することがわかる. あまり結果が変わらないからか, ロジスティック法則に言及している経済学の研究は少ないようだ.

############################
#
# Solow-Swan model with 
# logistic population growth
#
# Solow-Swan model
# f(k)=k^alpha
# s*f(k) -delta*k = (n+g)*k+k'
# k' = s*f(k) - (L'/L +g+delta)*k
# A'=gA
# Malthusian law
# L' = n*L
# Logistic law
# L'=r*L*(C-L)/K
#   =r*L-(r*C)*L^2
#
############################
library(ggplot2)  # ver. 1.0.1
library(ggthemes) # ver. 2.2.1
library(tidyr)    # ver. 0.2.0
library(dplyr)    # ver. 0.4.2

d.logistic <- function(x,r,C){
  x+ r*x +r*((C-x)/C)*x
}
f.prod <- function(x,a){
  x^a
}
# calculate Solow growth path
solow.growth.path <- function(span, AKLY0, type.L="Logistic", alpha, delta, growth, saving, 
                              n.growth=NULL, reproduction=NULL, Capacity=NULL){
  if(toupper(substr(type.L,1,1))=="M" ){
    if(is.null(n.growth)) stop("n.growth parameter needed if you select 'Malthusian' growth")
    is.Malthusian <- T
  }
  else if (toupper(substr(type.L,1,1))=="L"){
    if(is.null(Capacity) || is.null(reproduction) ) stop("Both 'reproduction' and 'Capacity' parameter needed if you select 'Logistic' growth ")
    is.Malthusian <- F
  }
  else stop("you must select type.L whether 'M' or 'L' ")
  X <- ts(data.frame(A   = ts(rep(NA,span)),
                     L   = ts(rep(NA,span)),
                     y   = ts(rep(NA,span)),
                     L.g = ts(rep(NA,span)),
                     k   = ts(rep(NA,span))
  ))
  # substitute initial values
  X[1,] <- c(AKLY0[-2], NA, AKLY0[2]/(AKLY0[1]*AKLY0[3]) )
  # Malthusian population growh law
  
  for(t in 2:span) {
    # Logistic population growth law
    if(is.Malthusian){
      X[t,"L"]   <- (1+n.growth)*X[t-1,"L"]  
    }
    else{
      X[t,"L"]   <- d.logistic(X[t-1,"L"],r = reproduction, C=Capacity)
    }
    X[t,"L.g"] <- (X[t,"L"]-X[t-1,"L"])/ X[t-1,"L"]
    X[t,"k"]   <- X[t-1,"k"] + saving*f.prod(X[t-1,"k"],a = alpha) - (growth + delta + X[t,"L.g"] ) * X[t-1,"k"]
    X[t,"A"]     <- X[t-1,"A"] + growth*X[t-1,"A"]
  }
  # state variables
  X[,"y"] <- f.prod(X[,"k"],a=alpha)
  return(X)
}

L0 <- 100
K0 <- 100
A0 <- 100

span <- 50
alpha <- .4
delta <- .3
growth <- .3
saving <- .4
# logistic parameter
reproduction <- .2
Capacity <- 500
# Malthusian parameter
n.growth <- reproduction


# Logistic
x.l <- solow.growth.path(span = span, AKLY0 =  c(A0, K0, L0, f.prod(x = K0/(A0*L0),alpha) ),
                         type.L = "L", alpha = alpha, delta = delta,
                         growth = growth, saving = saving, reproduction=reproduction, Capacity=Capacity
                        )
plot(x.l, main="人口成長の余地があるときの経済成長")
# over capacity case
Capacity <- 70
x.l.over <- solow.growth.path(span = span, AKLY0 =  c(A0, K0, L0, f.prod(x = K0/(A0*L0),alpha) ),
                              type.L = "L", alpha = alpha, delta = delta,
                              growth = growth, saving = saving, reproduction=reproduction, Capacity=Capacity
                             )
plot(x.l.over, main="人口成長が停滞するときの経済成長")

# Malthusian
Capacity <- 1000
x.m <- solow.growth.path(span = span, AKLY0 =  c(A0, K0, L0, f.prod(K0/(A0*L0),alpha) ),
                         type="M", alpha = alpha, delta = delta,
                         growth = growth, saving = saving, n.growth=n.growth
                        )
plot(x.m, main="マルサス法則を想定した経済成長")

# 比較
df1 <- as.data.frame(x.l); df1$t <- 1:nrow(df1);      df1$situation <- "成長期"
df2 <- as.data.frame(x.l.over); df2$t <- 1:nrow(df2); df2$situation <- "停滞期"
df3 <- as.data.frame(x.m); df3$t <- 1:nrow(df3); df3$situation <- "Malthusian"
df <- rbind(df1,df2, df3) %>% select(-A) %>% gather(key=series, value=val,-t, -situation) 
                          %>% mutate(col=paste0(situation, series))
ggplot(filter(df, !series %in%c("L","L.g") )) + geom_line(aes(x=t, y=val, group=col, color=series, linetype=situation))  +
theme_economist() + scale_color_colorblind() + theme(legend.title=element_blank()) + labs(y='')



# 参考: 人口成長のプロット
g.stasis <- ggplot(filter(df, situation=='停滞期', series=='L')) + geom_line(aes(x=t, y=val, group=col))
g.growth <- ggplot(filter(df, situation=='成長期', series=='L')) + geom_line(aes(x=t, y=val, group=col))
g.m <- ggplot(filter(df, situation=='Malthusian', series=='L')) + geom_line(aes(x=t, y=val, group=col))

g.style <- theme_classic() + theme(axis.title.y=element_text(angle=0, vjust=1))
g.stasis + labs(x="t", y="人口", title='ロジスティック (停滞期)')  + g.style
g.growth + labs(x="t", y="人口", title='ロジスティック (成長期)')  + g.style
g.m      + labs(x="t", y="人口", title='マルサス法則') + g.style

参考文献

Ritelli, D., Mingari Scarpello, G., & Brida, J. G. (2008). The Solow model with logistic manpower: a stability analysis. Journal of World Economics Review, 3(2), 161–166.