Daisuke Oyama
Faculty of Economics, University of Tokyo
In [1]:
%matplotlib inline
In [2]:
import time
import numpy as np
import matplotlib.pyplot as plt
import quantecon as qe
from quantecon.compute_fp import _print_after_skip
from quantecon.game_theory import Player, NormalFormGame, lemke_howson
In [3]:
def compute_fixed_point_ig(f, x_init, error_tol=1e-3, max_iter=50, verbose=1,
*args, **kwargs):
_skip = 1
if verbose:
start_time = time.time()
_print_after_skip(_skip, it=None)
x_new = x_init
iterate = 0
y_new = f(x_new, *args, **kwargs)
error = np.max(np.abs(y_new - x_new))
if error <= error_tol or iterate >= max_iter:
if verbose:
etime = time.time() - start_time
_print_after_skip(_skip, iterate, error, etime)
return x_new
X = np.array([x_new])
Y = np.array([y_new])
x_new = Y[0]
iterate += 1
while True:
y_new = f(x_new, *args, **kwargs)
error = np.max(np.abs(y_new - x_new))
if error <= error_tol or iterate >= max_iter:
break
X = np.append(X, np.expand_dims(x_new, axis=0), axis=0)
Y = np.append(Y, np.expand_dims(y_new, axis=0), axis=0)
m = len(X)
D = np.expand_dims(X, axis=1) - Y
D *= D
A = np.add.reduce(np.atleast_3d(D), axis=-1) * (-1)
B = np.identity(m)
g = NormalFormGame((Player(A), Player(B)))
_, rho = lemke_howson(g, init_pivot=m-1)
x_new = rho.dot(Y)
iterate += 1
if verbose:
etime = time.time() - start_time
_print_after_skip(_skip, iterate, error, etime)
return x_new
In [4]:
# Just a warmup
compute_fixed_point_ig(lambda x: 0.5*x, 1)
Out[4]:
Let us try the logistic function which is well known to generate chaotic behavior.
In [5]:
def logistic(x, r):
return r * x * (1 - x)
In [6]:
x = np.linspace(0, 1, 100)
y = logistic(x, r=4)
fig, ax = plt.subplots()
ax.plot(x, y)
ax.plot([0, 1], [0, 1], ':', color='k')
ax.set_aspect(1)
plt.show()
In [7]:
tol = 1e-5
x_init = 0.99
compute_fixed_point_ig(logistic, x_init, error_tol=tol, r=4)
Out[7]:
Comare compute_fixed_point from quantecon:
In [8]:
qe.compute_fixed_point(logistic, x_init, error_tol=tol, r=4)
Out[8]:
In [9]:
def f(x, M, c):
return -np.arctan(np.dot(M, (x - c)**3)) + c
In [10]:
x_min, x_max = -np.pi/2, np.pi/2
x = np.linspace(x_min, x_max, 100)
M = np.abs(np.random.randn())
c = 0
y = f(x, M, c)
fig, ax = plt.subplots()
ax.plot(x, y)
ax.set_xlim(x_min, x_max)
ax.set_ylim(x_min, x_max)
ax.plot([x_min, x_max], [x_min, x_max], ':', color='k')
ax.set_aspect(1)
plt.show()
In [11]:
n = 500
tol = 1e-5
max_iter = 200
In [12]:
num_trials = 3
for i in range(num_trials):
print("===== Experiment {} =====\n".format(i))
c = np.random.standard_normal(n)
M = np.abs(np.random.standard_normal(size=(n, n)))
x_init = (np.random.rand(n)-1/2)*np.pi + c
print("***Imitation game***")
x0 = compute_fixed_point_ig(f, x_init, tol, max_iter, M=M, c=c)
print("")
print("***Function iteration***")
x1 = qe.compute_fixed_point(f, x_init, tol, max_iter, verbose=1, print_skip=200, M=M, c=c)
print("")
In [13]:
num_trials = 3
for i in range(num_trials):
print("===== Experiment {} =====\n".format(i))
c = np.random.standard_normal(n)
M = np.random.normal(0, 1/13, size=(n, n))
x_init = (np.random.rand(n)-1/2)*np.pi + c
print("***Imitation game***")
x0 = compute_fixed_point_ig(f, x_init, tol, max_iter, M=M, c=c)
print("")
print("***Function iteration***")
x1 = qe.compute_fixed_point(f, x_init, tol, max_iter, verbose=1, print_skip=200, M=M, c=c)
print("")
In [ ]: