Commodity Pricing Model
Corresponds to section 6.3 in EDTC
Presentation here is terse--see that section for details
Outline of the Model
Market for a commodity whose price at t is pt
Demand
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
Supply
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
Equilibrium
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
- A function p of the state variable Xt
Analysis in the text shows that such a function is given by the p* that solves
where
Computing the Solution
To compute p we introduce the pricing function operator T*:
Let p be a guess of p*
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, p* is the unique fixed point.
Problem:
Using the code in p*
cpdynam.py, compute an approximation to - Bear in mind that the functions are defined on S (see above)
- Start your iteration at the inverse demand function P
- Now iterate with T
- Approximate at each step using linear interpolation
- Iterate until supremum distance between successive iterates is < 0.005
- See the previous lecture if you don't know how to do this
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()
0 comments:
Post a Comment