%matplotlib inline
1部門モデルを考える.
生産関数と効用関数を特定化して、オイラー条件
\begin{align} u_2(k_{t-1}, k_t) + \rho u_1(k_t, k_{t+1}) = 0 \end{align}によって定まる力学系を可視化したい
3つのプログラミング・パラダイム
## 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):
return A * k ** α
from math import log
def U(c):
return log(c)
def u(x, y):
return U(f(x) - y)
より
\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}def F(x, y):
return ((1 + ρ * α) * A * y ** α
- ρ * α * (A ** 2) * (x ** α) * (y ** (α - 1)))
def G(x):
return [
x[1],
F(x[0], x[1])
]
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]
以上をまとめたものが listing01.py
%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
で呼び出された場合には実行されます.
# listing02.py
import numpy as np
import matplotlib.pyplot as plt # プロッティング・ライブラリを読み込む
from listing01 import f
fig = plt.figure()
ax = fig.add_subplot(111) # 作図領域の作成
grids = np.linspace(0.0, 1.2, 120)
ax.plot(grids, f(grids))
plt.show() # 図の表示
# listing03.py
import matplotlib.pyplot as plt
from listing01 import f, G
duration = 10
x0 = [0.8, 0.43]
fig = plt.figure()
ax = fig.add_subplot(111)
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()
これでは動きが分からない ....
# 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()
ax = fig.add_subplot(111, aspect='equal')
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')
ax.arrow(x[0], x[1], dx[0], dx[1], length_includes_head=True)
x = x1
plt.show()
対象と動作をまとめてやる. たとえば次のような記法でシミュレーションができるようにしたい.
>>> 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]]
# 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] ** α)
]
# 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)
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]]
# 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')
ax.arrow(x0[0], x0[1], dx[0], dx[1], length_includes_head=True)
x0 = x[:]
初期値を変えてシミュレーションしたものを重ねることも簡単にできる. Simulation
クラスを書き換える.
# 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')
ax.arrow(x0[0], x0[1], dx[0], dx[1], length_includes_head=True)
x0 = x[:]
%run listing06
初期値を決めたときに、軌道を計算・プロットする方法の説明は終わり.
残る問題は安定多様体上の軌道をどのように計算するか?
長期均衡に収束するような初期値をなんとか探せばよい
あるいは、このモデルなら長期均衡の近傍から出発する逆向きの軌道を計算する方がスマートにできる. 生産関数と効用関数の特定化次第では逆向きの時間発展が計算できる. ヤコビアンの安定固有ベクトル方向に摂動を加えて、逆向きに時間発展させれば安定多様体上の軌道が計算できるはず