# ラムゼーモデルの局所動学¶

### 経済動学 講義資料¶

2015-5-25

https://bitbucket.org/kenjisato/economic_dynamics

In [1]:
%matplotlib inline


## 目標¶

1部門モデルを考える.

\begin{align} u_2(k_{t-1}, k_t) + \rho u_1(k_t, k_{t+1}) = 0 \end{align}

によって定まる力学系を可視化したい

3つのプログラミング・パラダイム

1. 手続き型
2. オブジェクト指向
3. [Optional] 関数型

## 手続き型 (素朴な方法)¶

### パラメータ¶

In [2]:
## Unicode identifiers are only valid in Python3
## In Python2, use less readable notation such as
# alpha = 0.4
# rho = 0.9

A = 1.1
α = 0.4
ρ = 0.9


### 生産関数¶

\begin{align} f(k) = Ak^\alpha \end{align}
In [3]:
def f(k):
return A * k ** α


### 消費の効用関数¶

\begin{align} U(c) = \log c \end{align}
In [4]:
from math import log

def U(c):
return log(c)


### 規約型効用関数¶

\begin{align} u(x, y) = U(f(x) - y) \end{align}
In [5]:
def u(x, y):
return U(f(x) - y)


### オイラー条件¶

\begin{align} u_2(k_{t-1}, k_t) = \frac{-1}{Ak_{t-1}^\alpha - k_t}\\ u_1(k_t, k_{t+1}) = \frac{\alpha A k_t^{\alpha - 1}}{Ak_t^\alpha - k_{t+1}} \end{align}

より

\begin{align} k_{t+1} = F(k_{t-1}, k_t) = (1+\rho \alpha) A k_t^\alpha - \rho \alpha A^2 k_{t-1}^\alpha k_t^{\alpha-1} \end{align}
In [6]:
def F(x, y):
return ((1 + ρ * α) * A * y ** α
- ρ * α * (A ** 2) * (x ** α) * (y ** (α - 1)))


### 相空間上の力学系¶

In [7]:
def G(x):
return [
x[1],
F(x[0], x[1])
]


### シミュレーション¶

In [8]:
duration = 4
x0 = [0.8, 0.43]

x = x0[:]
for t in range(duration):
print(x)
x = G(x)

[0.8, 0.43]
[0.43, 0.4063168783383537]
[0.4063168783383537, 0.509937247085211]
[0.509937247085211, 0.6875981937381642]


In [9]:
%cat listing01.py

#!/usr/bin/env python3

from math import log

## Unicode identifiers are only valid in Python3
## In Python2, use less readable notation such as
# alpha = 0.4
# rho = 0.9

A = 1.1
α = 0.4
ρ = 0.9

def f(k):
"""Production function"""
return A * k ** α

def U(c):
"""Utility function"""
return log(c)

def u(x, y):
"""Reduced form utility function"""
return U(f(x) - y)

def F(x, y):
"""Solution of Euler equation"""
return ((1 + ρ * α) * A * y ** α
- ρ * α * (A ** 2) * (x ** α) * (y ** (α - 1)))

def G(x):
"""Dynamical system"""
return [
x[1],
F(x[0], x[1])
]

if __name__ == "__main__":

duration = 4
x0 = [0.8, 0.43]

x = x0[:]
for t in range(duration):
print(x)
x = G(x)


### 補足¶

if __name__ == "_main__":
[do something]



これはPythonのイディオムで、 import文で呼び出された場合には実行されずに、run で呼び出された場合には実行されます.

### 可視化の基本¶

In [10]:
# listing02.py

import numpy as np
import matplotlib.pyplot as plt   # プロッティング・ライブラリを読み込む
from listing01 import f

fig = plt.figure()

grids = np.linspace(0.0, 1.2, 120)
ax.plot(grids, f(grids))

plt.show()   # 図の表示


### 軌道の可視化¶

In [11]:
# listing03.py

import matplotlib.pyplot as plt
from listing01 import f, G

duration = 10
x0 = [0.8, 0.43]

fig = plt.figure()

grids = np.linspace(0.0, 1.2, 120)
ax.plot(grids, f(grids))

x = x0[:]
for t in range(duration):
ax.plot(x[0], x[1], marker='o', linestyle='')
x = G(x)

plt.show()


これでは動きが分からない ....

In [12]:
# listing04.py

import numpy as np
import matplotlib.pyplot as plt
from listing01 import f, G

duration = 10
x0 = [0.8, 0.43]

fig = plt.figure()

grids = np.linspace(0.0, 1.2, 120)
ax.plot(grids, f(grids))

x = x0[:]
for t in range(duration):
x1 = G(x)
dx = [x1[0] - x[0], x1[1] - x[1]]
ax.plot(x[0], x[1], marker='', linestyle='', color='black')
x = x1
plt.show()


## OOP¶

>>> ramsey = Ramsey(A=1.1, α=0.4, ρ=0.9)
>>> sim = Simulation(ramsey, x0=[0.8, 0.43], duration=5)
>>> path = [x for x in sim]
>>> path
[[0.8, 0.43],
[0.43, 0.4063168783383537],
[0.4063168783383537, 0.5099372470852112],
[0.5099372470852112, 0.6875981937381646],
[0.6875981937381646, 0.8712750948974344]]
In [13]:
# listing05.py の一部

class Ramsey:
"""One-sector Ramsey model"""

def __init__(self, A, α, ρ):
self.A = A
self.α = α
self.ρ = ρ

def f(self, x):
return self.A * x ** self.α

def U(self, x):
return log(x)

def u(self, x, y):
return self.U(self.f(x) - y)

def forward(self, x):
"""1ステップの時間発展"""
A, α, ρ = self.A, self.α, self.ρ
return [
x[1],
(1 + ρ * α) * A * x[1] ** α - ρ * α * (A ** 2) * (x[1] ** (α - 1)) * (x[0] ** α)
]

In [14]:
# listing05.py の一部

class Simulation:
"""Simulation of a dynamical system"""

def __init__(self, sys, x0, duration):
self.sys = sys
self.x0 = x0
self.duration = duration

def __iter__(self):
x = self.x0[:]
for _ in range(self.duration):
yield x
x = self.sys.forward(x)

In [15]:
ramsey = Ramsey(A=1.1, α=0.4, ρ=0.9)
sim = Simulation(ramsey, x0=[0.8, 0.43], duration=5)

path = [x for x in sim]
path

Out[15]:
[[0.8, 0.43],
[0.43, 0.4063168783383537],
[0.4063168783383537, 0.5099372470852112],
[0.5099372470852112, 0.6875981937381646],
[0.6875981937381646, 0.8712750948974344]]
In [16]:
# listing05.py の一部

ramsey = Ramsey(A=1.1, α=0.4, ρ=0.9)
sim = Simulation(ramsey, x0=[0.8, 0.43], duration=10)

fig, ax = plt.subplots(subplot_kw={'aspect':'equal'})

grids = np.linspace(0.0, 1.2, 120)
ax.plot(grids, ramsey.f(grids))

for i, x in enumerate(sim):
if i == 0:
ax.plot(x[0], x[1], marker='', linestyle='', color='black')
else:
dx = [x[0] - x0[0], x[1] - x0[1]]
ax.plot(x[0], x[1], marker='', linestyle='', color='black')
x0 = x[:]


In [17]:
# listing06.py の一部

class Simulation:
"""__init__ と __iter__ は変更なし"""

def reset(self, *, x0=None, duration=None):
if x0 is not None:
self.x0 = x0[:]
if duration is not None:
self.duration = duration

def plot(self, ax):
for i, x in enumerate(self):
if i == 0:
ax.plot(x[0], x[1], marker='', linestyle='', color='black')
else:
dx = [x[0] - x0[0], x[1] - x0[1]]
ax.plot(x[0], x[1], marker='', linestyle='', color='black')
x0 = x[:]

In [18]:
%run listing06


あるいは、このモデルなら長期均衡の近傍から出発する逆向きの軌道を計算する方がスマートにできる. 生産関数と効用関数の特定化次第では逆向きの時間発展が計算できる. ヤコビアンの安定固有ベクトル方向に摂動を加えて、逆向きに時間発展させれば安定多様体上の軌道が計算できるはず

In [ ]: