Commodity Pricing Model

Corresponds to section 6.3 in EDTC –see for details

Outline of the Model

Market for a commodity whose price at time t is \(p_t\)

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 \(I_t\) denote the quantity purchased by speculators at time t

This quantity depreciates at rate \(\alpha\) between periods

Supply

The quantity supplied to the market at time t is denoted \(X_t\)

The “harvest” of the commodity at time t is \(W_t\)

We assume that

(1)\[\begin{split}(W_t)_{t \geq 1} \sim \phi \text{ on } S := [a,\infty), \; \text{ where } a > 0\end{split}\]

Total supply is the sum of the harvest plus supply carry-over from speculation

(2)\[X_t = \alpha I_{t-1} + W_t\]

Equilibrium

As shown in the text, in equilibrium we have

(3)\[\alpha \mathbb{E}_t p_{t+1} - p_t \leq 0\]

(4)\[\begin{split}\alpha \mathbb{E}_t p_{t+1} - p_t < 0 \text{ implies } I_t = 0\end{split}\]

and

(5)\[X_t = D(p_t) + I_t\]

The initial condition \(X_0\) is treated as given

Equilibrium prices and quantities are found by computing a pricing function

• A function \(p\) of the state variable \(X_t\)

Analysis in the text shows that such a function is given by the \(p^*\) that solves

(6)\[p^*(x) = \max \left\{ \alpha \int p^*(\alpha I(x) + z) \phi(z)dz, P(x) \right\} \qquad (x \in S)\]

where

(7)\[I(x) := x - D(p^*(x)) \qquad (x \in S)\]

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

(8)\[r = \max \left\{ \alpha \int p(\alpha (x-D(r)) + z) \phi(z)dz, \; P(x) \right\}\]

and

(9)\[P(x) \leq r \leq \max \left\{\alpha \int p(z) \phi(z)dz, P(x) \right\}\]

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,

(10)\[r = P(x) \text{ whenever } \alpha \int p(z)\phi(z)dz \leq P(x)\]

and

(11)\[r = \alpha \int p(\alpha (x-D(r)) + z) \phi(z)dz \; \text{ otherwise}\]

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

(12)\[W = c U + a, \; U \sim \text{Beta}(5, 5), \quad D(x) = 1/x\]

The listing corresponds to listing 6.7 in the text

Theorem: The operator \(T\) is a contraction mapping on the set

(13)\[\mathcal{C} := \{q \colon S \to \mathbb{R} \,|\, q \text{ continuous, decreasing, } q \geq P\}\]

with the supremum metric. Moreover, \(p^*\) is the unique fixed point.

Problem:

Using the code in cpdynam.py, compute an approximation to \(p^*\)

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