ill-identified diary

所属組織の業務や見解とは一切無関係なアフィリエイト付きメモ帳。所属とは関係ないけどここを見て所属先にも興味を持っていただけると喜びます。

[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.

マクロ経済学

マクロ経済学

Romer は第3版の邦訳も存在する:

Romer, D. (2012). Advanced Macroeconomics (4th ed.). The McGraw-Hill Companies.

Advanced Macroeconomics (The Mcgraw-Hill Series in Economics)

Advanced Macroeconomics (The Mcgraw-Hill Series in Economics)

Romer のマクロ経済学のテキストは、第3版の邦訳も存在する.

上級マクロ経済学

上級マクロ経済学

二神孝一・堀敬一. (2009). マクロ経済学. 有斐閣.

マクロ経済学

マクロ経済学

*1:この r は再生産 (reproduction) パラメータとも呼ばれる.

*2:Pierre-François Verhulst (1804-1849)

*3:変数分離形で解くことができる. 読者の練習問題とする.

*4:ここではソローモデル自体の解説は行わない. ネット上にも多くの有用な資料があるためそちらを参照されたい.

*5:多くのマクロ経済学の教科書では, n で表されていることが多いので注意.