ソロー・モデル

in English or the language of your choice.

import japanize_matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import py4macro
import statsmodels.formula.api as sm

# 警告メッセージを非表示
import warnings
warnings.filterwarnings("ignore")

はじめに

前章では差分方程式について説明し，簡単な経済モデルを使いコードの書き方を説明した。本章では，差分方程式の応用となるソロー・モデルを考える。ソロー・モデルは1987年にノーベル経済学賞を受賞したRobert M. Solowによって考案された経済成長モデルであり，マクロ経済学の代表的な理論モデルの一つである。今でも盛んに研究が続く経済成長のバックボーン的な存在である。モデルの説明の後，Pythonを使い動学的な特徴を明らかにし，線形近似を使った安定性の確認もおこなう。また理論的な予測がデータと整合性があるかについてもPenn World Tableを使って検討する。本章の内容は次章で議論する所得収斂の分析の基礎となる。

モデルの説明

ここではモデルの具体的な説明については教科書に譲るとして，簡単にモデルを紹介し，重要な式をまとめることにする。

＜記号＞

• 産出量：\(Y_t\)

• 消費量：\(C_t\)

• 投資量：\(I_t\)

• 資本ストック：\(K_t\)

• 労働：\(L_t\)

• 貯蓄率（一定な外生変数）：\(0<s<1\)

• 労働人口増加率（一定な外生変数）：\(n\equiv\dfrac{L_{t+1}}{L_t}-1\geq 0\)

• 資本減耗率（一定な外生変数）：\(0<d<1\)

• 生産性（一定な外生変数）：\(A>0\)

＜一人当たりの変数＞

• 一人当たり産出量：\(y_t\equiv\dfrac{Y_t}{L_t}\)

• 一人当たり消費量：\(c_t\equiv\dfrac{C_t}{L_t}\)

• 一人当たり投資量：\(i_t\equiv\dfrac{I_t}{L_t}\)

• 一人当たり資本ストック：\(k_t\equiv\dfrac{K_t}{L_t}\)

全ての市場は完全競争である閉鎖経済を考えよう。この経済には一種類の財（ニューメレール財）しかなく，消費・貯蓄・投資に使われる。財は次の生産関数に従って生産される。

(27)\[ Y_t=AK_t^aL_t^{1-a},\quad 0<a<1 \]

両辺を\(L_t\)で割ると一人当たりの変数で表した生産関数となる。

(28)\[ y_t=Ak_t^a \]

消費者は所得の割合\(s\)を貯蓄するが，このモデルの中で消費者の役割はこれだけであり，残り全ては生産側で決定される。貯蓄は\(sY_t\)であり投資\(I_t\)と等しくなる。

(29)\[ sY_t=I_t \]

\(t\)期の投資により\(t+1\)期の資本ストックは増加するが，毎期ごと資本は\(d\)の率で減耗する。即ち，投資と資本ストックには次の関係が成立する。

(30)\[ K_{t+1}-K_{t}=I_t-dK_t \]

ここで左辺は資本ストックの変化であり，右辺は純投資である。\(I_t\)は粗投資，\(dK_t\)は減耗した資本である。式(29)，(30)，(28)を使うと資本の蓄積方程式が導出できる。

(31)\[ K_{t+1} = sAK_t^aL^{1-a} + (1-d)K_t \]

両辺を\(L_t\)で割ることで一人当たりの変数で表すことができる。右辺は単純に一人当たりの変数に直し，左辺は次のように書き換えることに注意しよう。

\[ \frac{K_{t+1}}{L_t}=\frac{K_{t+1}}{L_{t+1}}\frac{L_{t+1}}{L_t} =k_{t+1}(1+n) \]

従って，\(t+1\)期の一人当たりの資本ストックを決定する式は次式で与えられる。

(32)\[ k_{t+1} = \frac{sAk_t^a + (1-d)k_t}{1+n} \]

この式は非線形の差分方程式だが，前章でも述べたように，考え方は線形差分方程式と同じであり，数値計算のためのPythonコードに関しては大きな違いはない。また式(32)を資本ストックの成長率を示す式に書き換えることもできる。

(33)\[ \frac{k_{t+1}}{k_t}-1 = \frac{sAk_t^{-(1-a)}-(n+d)}{1+n} \]

この式から資本が蓄積され\(k_t\)が増加すると，その成長率は減少していくことが分かる。

産出量の動学は生産関数(28)を使うことによって\(y_t\)の動きを確認できる。例えば，産出量の成長率を考えてみよう。式(28)を使うと

(34)\[ \frac{y_{t+1}}{y_t}-1 =\left(\frac{k_{t+1}}{k_t}\right)^{\alpha}-1 \]

となり，一人当たり資本ストックの成長率と同じような動きをすることが分かると思う。

次に定常状態を考えよう。定常状態では資本ストックは一定なり，資本の成長率である式(33)の右辺はゼロになる。 定常値は次のように確認することができる。

(35)\[\begin{split} \begin{align*} \frac{k_*}{k_*}-1 &=\frac{sAk_*^{-(1-a)} - (n+d)}{1+n} \\ &\Downarrow \\ k_* &=\left( \frac{As}{n+d} \right)^{\frac{1}{1-a}} \end{align*} \end{split}\]

この値を生産関数(28)に代入することにより一人当たりGDPの定常値を求めることができる。

(36)\[ y_*=Ak_*^{a} \]

Fig. 7は資本ストックの動学的均衡を示している。

Fig. 7 一人当たり資本ストックの動学

動学

差分方程式(32)を使って資本ストックの変化をプロットするが，以前と同じようにDataFrameを生成する関数を定義しよう。

def solow_model(k0, A,a,s,n,d,T=20):
    """引数
            k0: 資本の初期値
            A: 生産性
            a: 資本の所得比率 (a<1)
            s: 貯蓄率 (s<1)
            n: 労働人口成長率（％）
            d: 資本減耗率 (d<1)
            T: ループによる計算回数
       戻り値
            資本と産出量からなるDataFrame"""
    
    k = k0
    y = A*k0**a
    
    k_list = [k]
    y_list = [y]

    for t in range(T):
        
        k = ( s*A*k**a+(1-d)*k )/( 1+n )
        y = A*k**a

        k_list.append(k)
        y_list.append(y)

    # DataFrameの作成
    dic = {'capital':k_list, 'output':y_list}
    df = pd.DataFrame(dic)
    
    return df 

引数に使うパラーメータには次の値を使ってプロットしてみよう。

df = solow_model(k0=200, A=10, a=0.3, s=0.3, n=0.02, d=0.05, T=100)
df.plot(subplots=True)
pass

異なる初期値を使って資本の変化をプロットしてみる。

initial_list = range(100,300,10)      # 1

fig, ax = plt.subplots()              # 2

for i in initial_list:                # 3
    df_temp = solow_model(k0=i, A=10, a=0.3, s=0.3, n=0.02, d=0.05, T=150)
    ax.plot('capital', data=df_temp)
    
ax.set(title='ソロー・モデルの移行過程',  # 4
       xlabel='ループの回数',
       ylabel='一人当たり資本ストック')
pass

コードの説明

1. range(100,300,10)は100から300-1=299までの間の数を10間隔で生成する。100，110，120，\(cdots\)，290となる。

2. .subplots()を使ってキャンバスをfigに，「軸」をaxに割り当てる。

3. initial_listに対してforループを設定する。一回のループごとに以下をおこなう。

   • iはk0の引数に使う値になる。

   • solow_model()で資本ストックの数値を計算し，そのDataFrameをdf_tempに割り当てる。

   • .plot()を使いaxに図をプロットする。これにより図が重ねて描かれる。

4. axのメソッドset()を使い，引数を使い以下を追加する。

   • title：図のタイトル

   • xlabel：横軸のラベル

   • ylabel：縦軸のラベル

   • この３行を次のように書いても同じ結果となる。

     ax.set_title('ソロー・モデルの移行過程')
     ax.set_xlabel('ループの回数')
     ax.set_ylabel('一人当たり資本ストック')
     

     この３つを.set()でまとめて書いたのが上のコードである。３つに分けて書く利点はフォントの大きさを指定できることだろう。

図から初期値に関わらず定常値に収束していることが分かる。即ち，定常状態である長期均衡は安定的である。

定常状態での変数の値

次に定常状態での変数の値を計算してみよう。

def calculate_steady_state(A,a,s,n,d):
    
    k_ss = ( s*A/(n+d) )**(1/(1-a))    
    y_ss = A*k_ss**(a/(1-a))    
    
    return k_ss, y_ss

ss = calculate_steady_state(A=10, a=0.3, s=0.3, n=0.02, d=0.05)

print(f'定常状態での資本ストック:{ss[0]:.1f}'
      f'\n定常状態での産出量: {ss[1]:.1f}')

定常状態での資本ストック:214.5
定常状態での産出量: 99.8

線形近似

説明

Fig. 7は，定常状態は安定的であることを示している。またシミュレーションの結果からも定常状態の安定性が確認できる。次に，線形近似を使って解析的に安定性を確認してみることにする。また線形近似は真の値からの乖離が発生するが，その乖離がどの程度のものかをコードを使って計算することにする。

＜テイラー展開による１次線形近似＞

• 関数\(z=f(x)\)を\(x_*\)でテイラー展開すると次式となる。

  \[ z=f(x^*)+\left.\frac{df}{dx}\right|_{x=x_*}(x-x_*) \]

ソロー・モデルの式に当てはめると次のような対応関係にある。

• \(z\;\Rightarrow\;k_{t+1}\)

• \(f(x)\;\Rightarrow\;\dfrac{Ask_{t}^{\alpha}+(1-d)k_t}{1+n}\)

• \(x^*\;\Rightarrow\;k^*\)

公式に従って計算してみよう。

(37)\[\begin{split} \begin{align*} k_{t+1} &=\frac{Ask_{*}^{\alpha}+(1-d)k_*}{1+n} +\frac{1}{1+n}\left[\frac{Asa}{k_t^{1-\alpha}}+(1-d)\right]_{k_t=k_*}(k_t-k_*) \\ &=k_{*} +\frac{1}{1+n}\left[\frac{Asa}{k_*^{1-\alpha}}+(1-d)\right](k_t-k_*) \\ &=k_*+\frac{1}{1+n}\left[\frac{Asa}{As/(n+d)}+(1-d)\right](k_t-k_*) \\ &=k_*+\frac{a(n+d)+1-d-n+n}{1+n}(k_t-k_*) \\ &=k_*+\frac{1+n-(1-a)(n+d)}{1+n}(k_t-k_*) \\ &=(1-\lambda) k_t + \lambda k_* \end{align*} \end{split}\]

ここで

(38)\[ \lambda\equiv \frac{(1-a)(n+d)}{1+n} \]

３行目は定常状態の式(32)と式(35)を使っている。式(37)は\(k_{t+1}\)と\(k_t\)の線形差分方程式になっており，\(k_t\)の係数は

(39)\[ 0<1-\lambda<1 \quad\because 0<a,d<1 \]

が成立する。Fig. 7の赤い直線が式(37)である。従って，初期値\(k_0>0\)からスタートする経済は必ず長期的均衡に収束することがわかる。

式(38)は\(k_t\)の係数だが，その裏にあるメカニズムを考えてみよう。特に，\(a\)の役割を考える。\(a\)は資本の所得比率であり，目安の値は1/3である。そしてソロー・モデルにおける重要な役割が資本の限界生産性の逓減を決定することである。この効果により資本ストックが増加する毎に産出量も増加するがその増加自体が減少する。この効果により，Fig. 7の曲線は凹関数になっており，生産関数(28)の場合は必ず45度線と交差することになる。即ち，資本の限界生産性の逓減こそが\(k_t\)が一定になる定常状態に経済が収束す理由なのである。この点がソロー・モデルの一番重要なメカニズムとなる。

ここで\(a\)が上昇したとしよう。そうなると資本の限界生産性の逓減の効果は弱くなり，資本ストックが増加しても産出量の増加自体の減少は小さくなる。また式(35)が示すように定常状態での資本ストックはより大きくなる。これはFig. 7で曲線が上方シフトしている考えると良いだろう。

内生的成長

更に\(a\)を上昇させて\(a=1\)になるとどうなるのだろう。この場合，資本の限界生産性は逓減せず一定となる。そして\(k_t\)の係数である式(39)は1になってしまい，\(k_t\)が一定になる定常状態が存在しなくなる。\(a=1\)となる極限の状態を内生的成長と呼ぶ。ここでは立ち入った議論はしないが，内生的成長の典型的な生産関数は次式となり，

\[y_t=Ak_t\]

この生産関数に基づくモデルは\(AK\)モデルと呼ばれる。資本の限界生産性は\(A\)で一定になることが分かると思う。

\(\lambda\)の解釈

上の議論から\(\lambda\)が定常状態の安定性を決定することが分かったが，\(\lambda\)の値はどのように解釈できるだろうか。例えば，\(\lambda\)が大きい場合と小さい場合では何が違うのだろうか。次式は式(37)の最後の等号を少し書き換えたものである。

(40)\[ k_*-k_{t+1}=(1-\lambda)(k_*-k_t) \]

左辺は\(t+1\)期において定常状態までの残りの「距離」であり，右辺の\(k_*-k_t\)は\(t\)においての定常状態までの残りの「距離」である。後者を次の様に定義し

\[z_t\equiv k_*-k_t\]

式(40)を整理すると次式となる。

(41)\[ \frac{z_{t+1}-z_t}{z_t}=-\lambda \]

左辺は\(t\)期と\(t+1\)期において定常状態までの「距離」が何％減少したかを示す資本ストックの収束速度である。このモデルの中での収束速度の決定要因は資本の所得比率\(\alpha\)，労働人口増加率\(n\)と資本の減耗率\(d\)ということである。

ここでは\(a\)の役割に着目し，なぜ\(a\)の上昇は収束速度の減少をもたらすのかを直感的に考えてみよう。この点を理解するために，まず\(a\)は資本の限界生産性の逓減を決定するパラーメータであることを思い出そう。\(a\)が上昇するとその効果は弱まる。即ち，資本ストックが１単位増加すると産出量は増え，その増加分が減少するのが「逓減」であるが，その減少が小さくなるのである。これにより（上で説明したように）\(k_t\)が一定となる定常状態は増加することになる。重要な点は，定常状態の増加の意味である。定常状態はマラソンのゴールの様なものである。トップランナーはゴールすると走るのを止め，後続ランナーはトップランナーとの「距離」を縮めることができる。定常状態の増加は，ゴールが遠くなることと同じである。ゴールが遠のくとトップランナーは走り続けるわけだから，それだけ距離を縮めることが難しくなり収束速度が減少することになる。極端なケースとして\(a=1\)の場合，\(k_t\)が一定になる定常状態は存在せず，ゴールがない状態が永遠に続いており，永遠に収束しないということである。言い換えると，資本の限界生産性の逓減（\(a<1\)）こそが「距離」を縮めキャッチアップを可能にするメカニズムなのだ。

労働人口増加率\(n\)と資本の減耗率\(d\)の上昇は収束速度を速くする。式(35)から分かる様に，\(n\)もしくは\(d\)の上昇は定常状態を減少させる。即ち，ゴールはより近くになるということだ。

これである程度キャッチアップのメカニズムが分かったと思うが，今までの議論で足りないものが２点あるので，それらについて簡単に言及する。第一に，ここで考えたソロー・モデルには技術進歩が抜けている（一定な\(A\)を仮定した）。この点を導入してこそソロー・モデルのフルバージョンであり，その場合の労働効率１単位当たり資本ストック（\(K_t/(A_tL_t)\)）の収束速度は次の式で与えられる。

(42)\[ \lambda\equiv \frac{(1-a)(g+n+d+ng)}{(1+n)(1+g)} \]

ここで\(g\)は技術進歩率である。ソロー・モデルでは4つの変数が収束速度の決定要因になるる。\(g=0\)の場合，式(38)と同じになることが確認できる。第二に，式(38)は資本ストックの収束速度であり一人当たりGDPの収束速度と異なるのではないかという疑問である。実は同じである。これはコブ・ダグラス生産関数(28)を仮定しているからであり，対数の近似を使えば簡単に示すことができる。式(40)を次のように書き直そう。

(43)\[ \frac{k_{t+1}}{k_*}-1=(1-\lambda)\left(\frac{k_t}{k_*}-1\right) \]

ここで\(\log(1+x-1)\approx x-1\)の近似を使い左辺を次のように書き換える。

\[ \frac{k_{t+1}}{k_*}-1 \approx\log\left(\frac{k_{t+1}}{k_*}\right) =\log\left(\frac{y_{t+1}}{y_*}\right)^{\frac{1}{\alpha}} =\frac{1}{\alpha}\left(\frac{y_{t+1}}{y_*}-1\right) \]

同様に\(\dfrac{k_{t}}{k_*}-1\)もこの形に書き換えることができる。後はこの関係を使うことにより，式(43)を整理すると次式となる。

(44)\[ y_*-y_{t+1}=(1-\lambda)(y_*-y_t) \]

式(40)と同じ形になっているので所得の収束速度も式(41)と同じである。

線形近似による誤差

次に関数solow_model()を修正して線形近似の誤差を確かめてみよう。

def solow_model_approx(k0,A,a,s,n,d,T=20):
    """引数
            k0: 資本の初期値
            A: 生産性
            a: 資本の所得比率 (a<1)
            s: 貯蓄率 (s<1)
            n: 労働人口成長率（n>=0）
            d: 資本減耗率 (d<1)
            T: ループによる計算回数
       戻り値
            線形近似モデルを使い計算した資本と産出量からなるDataFrame"""
    
    k = k0
    y = A*k0**a
    
    k_list = [k]
    y_list = [y]

    # 定常状態
    k_ss = ( s*A/(n+d) )**(1/(1-a))    
    
    for t in range(T):
        
        lamb = 1-(1-a)*(n+d)/(1+n)  # lambda
        k = lamb*k+(1-lamb)*k_ss  # 線形近似
        y = A*k**a

        k_list.append(k)
        y_list.append(y)

    # DataFrameの作成
    dic = {'capital':k_list, 'output':y_list}
    df = pd.DataFrame(dic)
    
    return df 

近似誤差を計算するために，上と同じ数値でシミュレーションをおこなう。

df_approx = solow_model_approx(k0=200, A=10, a=0.3, s=0.3, n=0.02, d=0.05, T=150)

非線形のモデル(32)と線形近似のモデル(37)で計算された資本ストックを重ねて図示してみよう。

df_approx['capital'].plot(label='線形近似モデル', legend=True)
df['capital'].plot(label='ソロー・モデル', legend=True)
pass

df_approxが図示され、その上にdfが重ねて表示されるが、殆ど同じのように見える。誤差を％で計算し図示してみよう。

( 100*( 1-df_approx['capital']/df['capital'] ) ).plot(marker='.')
pass

初期の資本ストック\(k_0\)は同じなので誤差はゼロであり，\(t=1\)から線形近似の誤差が現れることになる。誤差は単調ではない。Fig. 7の図が示しているように，階段のような形で増加していくためであり，その階段お大きさや進み具合が異なるためである。線形近似の値は大き過ぎるため負の値になっているが、誤差は大きくても約0.02％であり、定常状態に近づくにつれて誤差はゼロに近づいている。もちろん，誤差の値は初期値が定常値から離れればそれだけ大きくなっていく。

長期均衡の予測

説明

この節では長期均衡（定常状態）に焦点を当て，理論的な予測のデータとの整合性をチェックする。まず定常状態の特徴をまとめよう。式(35)と(36)を使いうと定常状態での一人当たり資本ストックとGDPは次式で与えれる。

\[ k_*= \left( \frac{sA}{n+d} \right)^{\frac{1}{1-a}}, \quad y_*=Ak_*^a \]

この２つをそれぞれ試すこともできるが，同時に捉えるために２つの式の比率を考える。

(45)\[ \frac{k_*}{y_*}=\frac{s}{n+d} =\left.\frac{K_t/L_t}{Y_t/L_t}\right|_{\text{定常状態}} =\left.\frac{K_t}{Y_t}\right|_{\text{定常状態}} \]

この値は資本ストック対GDP比と等しいく，次のことが分かる。

• 貯蓄率\(s\)の上昇は資本ストック対GDP比を増加させる。。

• 労働人口成長率\(n\)の上昇は資本ストック対GDP比を減少させる。

• 資本減耗率\(d\)の上昇は資本ストック対GDP比を減少させる。

この３つの予測が成立するか確かめるためにpy4macroモジュールに含まれるPenn World Dataの次の変数を使う。

• cgdpo：GDP（2019年;生産側）

• cn：物的資本ストック（2019年）

• csh_i：対GDP比資本形成の比率

  • 投資の対GDP比である。

  • 貯蓄率\(s\)の代わりに使う。

  • 1960年〜2019年の平均を使う。

• emp：雇用者数

  • 労働人口の代わりに使う。

  • 1960年〜2019年の平均成長率\(n\)の計算に使う。

• delta：資本ストックの年平均減耗率

  • 1960年〜2019年の平均を使う。

データ

1960年以降のデータをpwtに割り当てる。

pwt = py4macro.data('pwt').query('year >= 1960')

貯蓄率

1960年から2019年までの国別の貯蓄率の平均を計算するが必要がある。ここで説明しているDataFrameのメソッド.groupby()を使うのが最も簡単な計算方法だろう。ここでは異なる方法としてDataFrameのメソッド.pivot()を紹介する。.pivot()はデータを整形する上で非常に便利なメソッドなので知って損はないだろう。

.pivot()は，元のDataFrameから列を選び，その列から新たなDataFrameを作成する便利なメソッドである。実際にコードを実行して説明しよう。

saving = pwt.pivot(index='year', columns='country', values='csh_i')
saving.head()

country	Albania	Algeria	Angola	Anguilla	Antigua and Barbuda	Argentina	Armenia	Aruba	Australia	Austria	...	United Arab Emirates	United Kingdom	United States	Uruguay	Uzbekistan	Venezuela (Bolivarian Republic of)	Viet Nam	Yemen	Zambia	Zimbabwe
year
1960	NaN	0.303861	NaN	NaN	NaN	0.148988	NaN	NaN	0.356183	0.215909	...	NaN	0.250929	0.243929	0.154078	NaN	0.377831	NaN	NaN	0.053204	0.201060
1961	NaN	0.337310	NaN	NaN	NaN	0.151162	NaN	NaN	0.302833	0.205103	...	NaN	0.242859	0.235292	0.157036	NaN	0.344871	NaN	NaN	0.047308	0.184336
1962	NaN	0.372718	NaN	NaN	NaN	0.131594	NaN	NaN	0.320442	0.184987	...	NaN	0.224224	0.237818	0.133648	NaN	0.351696	NaN	NaN	0.042713	0.129202
1963	NaN	0.359685	NaN	NaN	NaN	0.113241	NaN	NaN	0.315644	0.184251	...	NaN	0.226413	0.243220	0.123745	NaN	0.324748	NaN	NaN	0.031628	0.113312
1964	NaN	0.272713	NaN	NaN	NaN	0.139779	NaN	NaN	0.341544	0.208624	...	NaN	0.270936	0.246120	0.102930	NaN	0.397737	NaN	NaN	0.020669	0.117670

5 rows × 183 columns

コードの説明

pwtの３つの列year，country，csh_iを使って新たなDataFrameを作成しsavingに割り当てている。引数は次の役割をする。

• index：新たなDataFrameの行ラベルを指定する。

  • コードではyearが指定され行ラベルになっている。

• columns：新たなDataFrameの列ラベルを指定する。

  • コードではcountryが指定され列ラベルになっている。

• values：新たなDataFrameの値を指定する。

  • コードではcsh_iが指定され，その値でDataFrameが埋め尽くされている。

1960年以降欠損値がない国だけを使うことにしよう。NaNがある列を削除する必要があるので.dropna()を使う。

saving = saving.dropna(axis='columns')
saving.head()

country	Algeria	Argentina	Australia	Austria	Bangladesh	Barbados	Belgium	Benin	Bolivia (Plurinational State of)	Botswana	...	Tunisia	Turkey	U.R. of Tanzania: Mainland	Uganda	United Kingdom	United States	Uruguay	Venezuela (Bolivarian Republic of)	Zambia	Zimbabwe
year
1960	0.303861	0.148988	0.356183	0.215909	0.015934	0.097599	0.290086	0.041995	0.155338	0.109739	...	0.119855	0.206398	0.110611	0.105968	0.250929	0.243929	0.154078	0.377831	0.053204	0.201060
1961	0.337310	0.151162	0.302833	0.205103	0.017991	0.080713	0.299679	0.038844	0.116987	0.126386	...	0.148971	0.200566	0.119069	0.116082	0.242859	0.235292	0.157036	0.344871	0.047308	0.184336
1962	0.372718	0.131594	0.320442	0.184987	0.026607	0.073272	0.281366	0.038015	0.194461	0.143890	...	0.187644	0.196202	0.099676	0.116963	0.224224	0.237818	0.133648	0.351696	0.042713	0.129202
1963	0.359685	0.113241	0.315644	0.184251	0.019315	0.079214	0.272968	0.039580	0.165575	0.160907	...	0.194986	0.201216	0.089344	0.146264	0.226413	0.243220	0.123745	0.324748	0.031628	0.113312
1964	0.272713	0.139779	0.341544	0.208624	0.024646	0.078267	0.309144	0.033927	0.171541	0.200124	...	0.210442	0.199917	0.119029	0.187897	0.270936	0.246120	0.102930	0.397737	0.020669	0.117670

5 rows × 111 columns

コードの説明

.dropna()はNaNがある行か列を削除する。行と列のどちらを削除するかは引数axisで指定するが，デフォルトは'rows'である。即ち，引数なしで.dropna()を実行するとNaNがある行が削除される。ここでは列を削除したいので，引数に'columns'を指定している。

AlbaniaやAngolaなどが削除されていることが確認できる。何ヵ国残っているか確認してみよう。

saving.shape

(60, 111)

コードの説明

属性.shapeはDataFrameの行の数（左の数字）と列（右の数字）の数を返す。

111ヵ国含まれていることが確認できた。次に，それぞれの列の平均を計算する。

saving = saving.mean().to_frame('saving_rate')   # 1
saving.head()

	saving_rate
country
Algeria	0.349476
Argentina	0.147927
Australia	0.285702
Austria	0.271485
Bangladesh	0.135590

コードの説明

1. .mean()はそれぞれの列の平均を計算し，Seriesを返す。後でDataFrameを結合するメソッド.merge()を使うために.to_frame()を使ってSeriesをDataFrameに変換しており，引数saving_rateは列ラベルを指定している。もちろん引数を使わずに２行に分けることも可能である。

saving = saving.mean().to_frame()
saving.columns = ['saving_rate']

Tip

上の計算では1960年以降に欠損値が一つでもあればその国は排除されたが，全ての年でデータが揃っている経済だけを扱いたい場合に便利に使える方法である。一方で，.groupby()を使うと欠損値があっても平均は計算されるので，単純に.groupby()を使うと1960年以降に欠損値がある経済も含まれることになる。それを避けるためには一捻り必要だが，それについては貯蓄率・資本減耗率・労働人口増加率の平均が参考になるだろう。

資本減耗率

savingと同じ方法で資本減耗率の平均からなるDataFrameを作成する。

depreciation = pwt.pivot(index='year', columns='country', values='delta')
depreciation = depreciation.dropna(axis='columns')

depreciation.shape

(60, 110)

110ヵ国含まれている。

depreciation = depreciation.mean().to_frame('depreciation')

労働人口成長率

平均成長率を計算するには1960年と2019年の労働人口だけで計算できるが，上と同じ方法で計算してみる。

emp = pwt.pivot(index='year', columns='country', values='emp')
emp = emp.dropna(axis='columns')

emp.shape

(60, 91)

91ヵ国しか含まれていない。

emp_growth = ( ( emp.loc[2019,:]/emp.loc[1960,:] )**(1/(len(emp)-1))-1 
             ).to_frame('employment_growth')
emp_growth.head()

	employment_growth
country
Algeria	0.030113
Argentina	0.016983
Australia	0.019336
Austria	0.004750
Bangladesh	0.022279

資本ストック対GDP比

2019年のcgdpoとcnを使って資本ストック対GDP比を含むDataFrameを作成する。

ky_ratio = pwt.query('year == 2019') \
              .loc[:,['country','cgdpo','cn']] \
              .set_index('country') \
              .dropna()
ky_ratio.head()

	cgdpo	cn
country
Aruba	3466.241943	1.742797e+04
Angola	223289.312500	1.299232e+06
Anguilla	241.384537	2.266816e+03
Albania	36288.328125	2.239714e+05
United Arab Emirates	635332.812500	4.453746e+06

コードの説明

.set_index(''country')を使ってcountryを行ラベルに指定し，.dropna()によって欠損値がある行は削除する。

資本ストック対GDP比の列の作成しよう。

ky_ratio['ky_ratio'] = np.log( ky_ratio['cn']/ky_ratio['cgdpo'] )

含まれる国数を確認する。

ky_ratio.shape

(180, 3)

180ヵ国含まれており，savingやdepreciationの国数よりも多くの国が含まれている。

データの結合

上で作成したDataFrameを結合する必要があり，そのためのPandasの関数.merge()の使い方を説明する。df_leftとdf_rightの２つのDataFrameがあるとしよう。df_leftを左のDataFrame，df_rightを右のDataFrameと呼ぶことにする。２つを結合する場合，次のコードとなる。

pd.merge(df_left, df_right) 

しかし注意が必要な点が２つある。

1. 行数が同じでもdf_leftとdf_rightでは行の並びが異なる可能性がある。

2. 行数が異なる可能性がある。

これらの問題に対応するためのに引数が用意されている。

まず１つ目の問題は，行ラベルを基準に，もしくはある列に合わせて行を並び替えることにより対応できる。例えば，上で作成したDataFrameであれば，行ラベルがcountryになっているので，それに合わせて結合すれば良い。その場合の引数を含めたコードは次の様になる。

pd.merge(df_left, df_right, left_index=True, right_index=True) 

ここでのleft_index=Trueとright_index=Trueは行ラベルを基準に結合することを指定しており，デフォルトは両方ともFalseである。行ラベルではなく，ある列を基準に結合したい場合もあるだろう。その場合は次の引数を使う。

pd.merge(df_left, df_right,
         left_on=＜`df_left`の基準列のラベル（文字列）＞,
         right_on=＜`df_right`の基準列のラベル（文字列）＞) 

left_onは基準列に使うdf_leftにある列ラベルを文字列で指定する。同様にright_onは基準列に使うdf_rightにある列ラベルを文字列で指定する。デフォルトは両方ともNoneとなっている。

２つ目の問題はhowという引数を使うことにより対処できる。使える値は次の４つであり，いずれも文字列で指定する。

• 'inner'：df_leftとdf_rightの両方の基準列にある共通の行だけを残す（デフォルト）。

• 'left'：df_leftの行は全て残し，df_rightからはマッチする行だけが残り，対応する行がない場合はNaNが入る。

• 'right'：df_rightの行は全て残し，df_leftからはマッチする行だけが残り，対応する行がない場合はNaNが入る。

• 'outer'：df_leftとdf_rightの両方の行を全て残し，マッチする行がない場合はNaNを入れる。

では実際に上で作成したDataFrameを結合しよう。

for df_right in [saving, depreciation, emp_growth]:  # 1
    ky_ratio = pd.merge(ky_ratio, df_right,          # 2
                        left_index=True,             # 3
                        right_index=True,            # 4
                        how='outer')                 # 5

コードの説明

1. df_rightが上の説明のdf_rightに対応している。[saving, depreciation, emp_growth]は上で作成したDataFrameのリスト。

2. ky_ratioが上の説明のdf_leftに対応している。右辺で結合したDataFrameを左辺にあるky_ratioに割り当てている（際割り当て）。

3. ky_ratioの行ラベルを基準とすることを指定する。

4. df_rightの行ラベルを基準とすることを指定する。

5. ここでの'outer'は左右のDataFrameのそれぞれの行ラベルを残し，値がない箇所にはNaNを入れることを指定する。

結合の結果を表示してみよう。

ky_ratio.head()

	cgdpo	cn	ky_ratio	saving_rate	depreciation	employment_growth
country
Albania	36288.328125	2.239714e+05	1.820022	NaN	NaN	NaN
Algeria	503848.312500	2.290438e+06	1.514223	0.349476	0.041502	0.030113
Angola	223289.312500	1.299232e+06	1.761060	NaN	NaN	NaN
Anguilla	241.384537	2.266816e+03	2.239740	NaN	NaN	NaN
Antigua and Barbuda	1631.650513	1.022517e+04	1.835260	NaN	NaN	NaN

列AlbaniaやAngolaはsaving，depreciation，emp_growthのDataFramには無いためNaNが入っている。

結合後のky_ratioの情報を表示してみよう。

ky_ratio.info()

<class 'pandas.core.frame.DataFrame'>
Index: 180 entries, Albania to Zimbabwe
Data columns (total 6 columns):
 #   Column             Non-Null Count  Dtype  
---  ------             --------------  -----  
 0   cgdpo              180 non-null    float64
 1   cn                 180 non-null    float64
 2   ky_ratio           180 non-null    float64
 3   saving_rate        111 non-null    float64
 4   depreciation       110 non-null    float64
 5   employment_growth  91 non-null     float64
dtypes: float64(6)
memory usage: 9.8+ KB

180ヵ国が含まれるがNaNがある国も多いことが分かる。

トレンド線と散布図

ここでは次の３つをおこなう。

• 資本ストック対GDP比と次の３つの変数の散布図の表示

  • 貯蓄率

  • 資本減耗率

  • 労働人口増加率

• 回帰分析に基づいて計算したトレンド線の表示

• トレンド線の傾きの統計的優位性の表示

forループを使ってこれらを同時に計算・表示する。まずky_ratioの列ラベルをみると，回帰分析の説明変数に使う変数が最後の３つに並んでいる。

ky_ratio.columns[-3:]

Index(['saving_rate', 'depreciation', 'employment_growth'], dtype='object')

これを使いforループを組んでみよう。

for var in ky_ratio.columns[-3:]:              # 1
    
    df_temp = ky_ratio.copy()                  # 2
    res = sm.ols(f'ky_ratio ~ {var}',          # 3
                 data=df_temp).fit()           # 4
    bhat = res.params.iloc[1]                  # 5
    pval = res.pvalues.iloc[1]                 # 6
    
    df_temp['Trend'] = res.fittedvalues        # 7

    fig, ax = plt.subplots(tight_layout=True)  # 8
    ax.scatter(var, 'ky_ratio', data=df_temp)  # 9
    ax.plot(var, 'Trend',                      # 10
            data=df_temp.sort_values('Trend'), # 11
            c='r')                             # 12
    ax.set_title(f'トレンドの傾き：{bhat:.2f}\n'  # 13
                 f'p値：{pval:.3f}', size=20)   # 14
    ax.set_ylabel('資本ストック対GDP比（対数）',   # 15
                  size=15)                     # 16
    ax.set_xlabel(f'{var}', size=20)           # 17

コードの説明

1. ky_ratioの最後の３列のラベルに対してforループを組んで，varはその列ラベルを指す。

2. ky_ratioのコピーを作りdf_tempに割り当てる。.copy()は実態としては別物のコピーを作成する。詳細は割愛するが.copy()がないと実態は同じで参照記号のみが異なることになり，予期しない結果につながることを防ぐために.copy()を使っている。

3. 回帰分析をおこなっているが，f-stringを使い回帰式の説明変数を指定している。また回帰分析の結果をresに割り当てている。

4. dataで回帰分析のデータの指定をおこない，.fit()で自動計算！

5. resの属性.paramsは推定値を２つ返すが，1番目に傾きの推定値が格納されているため[1]で抽出し，bhatに割り当てている。

6. resの属性.pvaluesは\(p\)値を２つ返すが，1番目に傾きの\(p\)値が格納されているため[1]で抽出し，pvalに割り当てている。

7. resの属性fittedvaluesは予測値を返すが，それを新たな列（ラベルはTrend）としてdf_tempに追加している。

8. .subplots()によるキャンバスをfigに，「軸」をaxに割り当てている。引数のtight_layout=Trueはスペースに余裕がある様に表示にするためであり，なくても構わない。

9. .scatter()は散布図を描くメソッド。横軸はvar，縦軸はky_ratioとし，df_tempをデータとして使う。

10. .plot()を使いトレンド線を描く。横軸はvar，縦軸はTrend。

11. データにはdf_tempを使うが，メソッド.sort_values('Trend')で列Trendに従って昇順に並び替える。これはトレンド線が正確に表示されるために必要。

12. トレンド線の色を赤に指定する。

13. 図のタイトルを設定する。

    • f-stringを使いbhatの値を代入する。

    • :.2fは小数点第２位までの表示を指定。

    • \nは改行の意味。

    • 行を変えるので\nの後に'が必要となる。

14. f-stringを使いpvalの値を代入する。:.3fは小数点第３位までの表示を指定し，size=20はフォントの大きさの指定。

15. 縦軸のラベルの設定。

16. フォントの大きさの指定。

17. 横軸のラベルの設定であり，f-stringを使いvarの値を代入している。size=20はフォントの大きさの指定。

３つの図からソロー・モデルの理論的予測はデータと整合性があることが確認できる。ここで注意する点が一つある。式(45)は因果関係を予測している。例えば，貯蓄率が高くなることにより長期的な一人 当たりGDPは増加する。一方，トレンド線は因果関係を示しているのではなく単なる相関関係を表している。ソロー・モデルの因果関係を計量経済学的に検討するにはさまざまな要因の検討が必要になり，本章の域を超える事になる。