Corresponds to section 6.3 in EDTC
Presentation here is terse--see that section for details
Market for a commodity whose price at t is pt
Commodity is purchased by consumers and speculators (investors)
Consumers demand quantity D(p) corresponding to price p
The inverse demand function P is assumed to be continuous, decreasing, etc. (see EDTC)
Let It denote the quantity purchased by speculators at t
This quantity depreciates at rate α between periods
The quantity supplied to the market at time t is denoted Xt
The "harvest" of the commodity at time t is Wt
We assume that
Total supply is the sum of the harvest plus supply carry-over from speculation
As shown in the text, in equilibrium we have
and
The initial condition X0 is treated as given
Equilibrium prices and quantities are found by computing a pricing function
Analysis in the text shows that such a function is given by the
where
To compute
Let p be a guess of
Consider the new function on S constructed by associating to each x in S the real number r satisfying
and
We denote the new function by Tp, where Tp(x) is the r that solves this equation
It can be shown that this r exists, is unique, and, moreover,
and
The next listing implements this logic to compute Tp(x) from p and x
## Filename: cpdynam.py
## Author: John Stachurski
from scipy import mean
from scipy.stats import beta
from scipy.optimize import brentq
alpha, a, c = 0.8, 5.0, 2.0
W = beta(5, 5).rvs(1000) * c + a # Shock observations
D = P = lambda x: 1.0 / x
def fix_point(h, lower, upper):
"""Computes the fixed point of h on [upper, lower] using SciPy's
brentq routine, which finds the zeros (roots) of a univariate function.
Parameters: h is a function and lower and upper are numbers."""
return brentq(lambda x: x - h(x), lower, upper)
def T(p, x):
"""Computes Tp(x), where T is the pricing functional operator.
Parameters: p is a vectorized function (i.e., acts pointwise on
arrays) and x is a number."""
y = alpha * mean(p(W))
if y <= P(x):
return P(x)
h = lambda r: alpha * mean(p(alpha*(x - D(r)) + W))
return fix_point(h, P(x), y)
Here we are assuming that
The listing corresponds to listing 6.7 in the text
Theorem: The operator T is a contraction mapping on the set
with the supremum metric. Moreover,
Problem:
Using the code in cpdynam.py, compute an approximation to
Plot the sequence of iterates that results
Your plot should look like this
The smallest (lowest) function is the starting point P
Solution:
import numpy as np
from lininterp import LinInterp
from cpdynam import *
import matplotlib.pyplot as plt
gridsize = 150
grid = np.linspace(a, 35, gridsize)
tol = 0.0005
vals = P(grid)
while 1:
plt.plot(grid, vals, 'k-')
p = LinInterp(grid, vals)
f = np.vectorize(lambda x: T(p, x))
new_vals = f(grid)
if max(np.abs(new_vals - vals)) < tol:
break
vals = new_vals
plt.show()