2015 経済動学 講義資料 (2015/05/25)

2015-05-25.ipynb

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)

listing02.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*

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

#!/usr/bin/env python3

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

#!/usr/bin/env python3

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()

listing05.py

#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt

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] ** α)
        ]


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)


if __name__ == "__main__":

    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[:]

listing06.py

import numpy as np
import matplotlib.pyplot as plt
from collections import namedtuple
from listing05 import Ramsey


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

    def __init__(self, sys, x0=None, duration=None):
        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)

    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[:]


if __name__ == "__main__":

    ramsey = Ramsey(A=1.1, α=0.4, ρ=0.9)
    sim = Simulation(ramsey)

    SimParam = namedtuple('SimParam', ('x0', 'duration'))
    params = [
        SimParam((1.2, 0.460), 10),
        SimParam((1.2, 0.430), 10),
        SimParam((1.2, 0.425), 6),
        SimParam((0.05, 0.15), 10)
    ]

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

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

    for param in params:
        sim.reset(x0=param.x0, duration=param.duration)
        sim.plot(ax)