28. Commodity Prices#

28.1. Outline#

For more than half of all countries around the globe, commodities account for the majority of total exports.

Examples of commodities include copper, diamonds, iron ore, lithium, cotton and coffee beans.

In this lecture we give an introduction to the theory of commodity prices.

The lecture is quite advanced relative to other lectures in this series.

We need to compute an equilibrium, and that equilibrium is described by a price function.

We will solve an equation where the price function is the unknown.

This is harder than solving an equation for an unknown number, or vector.

The lecture will discuss one way to solve a functional equation (an equation where the unknown object is a function).

For this lecture we need the yfinance library.

!pip install yfinance
Hide code cell output
Collecting yfinance
  Downloading yfinance-0.2.63-py2.py3-none-any.whl.metadata (5.8 kB)
Requirement already satisfied: pandas>=1.3.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from yfinance) (2.2.2)
Requirement already satisfied: numpy>=1.16.5 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from yfinance) (1.26.4)
Requirement already satisfied: requests>=2.31 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from yfinance) (2.32.3)
Collecting multitasking>=0.0.7 (from yfinance)
  Downloading multitasking-0.0.11-py3-none-any.whl.metadata (5.5 kB)
Requirement already satisfied: platformdirs>=2.0.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from yfinance) (3.10.0)
Requirement already satisfied: pytz>=2022.5 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from yfinance) (2024.1)
Collecting frozendict>=2.3.4 (from yfinance)
  Downloading frozendict-2.4.6-py312-none-any.whl.metadata (23 kB)
Collecting peewee>=3.16.2 (from yfinance)
  Downloading peewee-3.18.1.tar.gz (3.0 MB)
?25l     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/3.0 MB ? eta -:--:--
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 3.0/3.0 MB 45.3 MB/s eta 0:00:00
?25h
  Installing build dependencies ... ?25l-
 \
 |
 /
 done
?25h  Getting requirements to build wheel ... ?25l-
 done
?25h  Preparing metadata (pyproject.toml) ... ?25l-
 done
?25hRequirement already satisfied: beautifulsoup4>=4.11.1 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from yfinance) (4.12.3)
Collecting curl_cffi>=0.7 (from yfinance)
  Downloading curl_cffi-0.11.3-cp39-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (14 kB)
Requirement already satisfied: protobuf>=3.19.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from yfinance) (4.25.3)
Collecting websockets>=13.0 (from yfinance)
  Downloading websockets-15.0.1-cp312-cp312-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (6.8 kB)
Requirement already satisfied: soupsieve>1.2 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from beautifulsoup4>=4.11.1->yfinance) (2.5)
Requirement already satisfied: cffi>=1.12.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from curl_cffi>=0.7->yfinance) (1.17.1)
Requirement already satisfied: certifi>=2024.2.2 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from curl_cffi>=0.7->yfinance) (2024.8.30)
Requirement already satisfied: python-dateutil>=2.8.2 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from pandas>=1.3.0->yfinance) (2.9.0.post0)
Requirement already satisfied: tzdata>=2022.7 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from pandas>=1.3.0->yfinance) (2023.3)
Requirement already satisfied: charset-normalizer<4,>=2 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from requests>=2.31->yfinance) (3.3.2)
Requirement already satisfied: idna<4,>=2.5 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from requests>=2.31->yfinance) (3.7)
Requirement already satisfied: urllib3<3,>=1.21.1 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from requests>=2.31->yfinance) (2.2.3)
Requirement already satisfied: pycparser in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from cffi>=1.12.0->curl_cffi>=0.7->yfinance) (2.21)
Requirement already satisfied: six>=1.5 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from python-dateutil>=2.8.2->pandas>=1.3.0->yfinance) (1.16.0)
Downloading yfinance-0.2.63-py2.py3-none-any.whl (118 kB)
Downloading curl_cffi-0.11.3-cp39-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (8.5 MB)
?25l   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/8.5 MB ? eta -:--:--
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 8.5/8.5 MB 123.4 MB/s eta 0:00:00
?25hDownloading frozendict-2.4.6-py312-none-any.whl (16 kB)
Downloading multitasking-0.0.11-py3-none-any.whl (8.5 kB)
Downloading websockets-15.0.1-cp312-cp312-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl (182 kB)
Building wheels for collected packages: peewee
  Building wheel for peewee (pyproject.toml) ... ?25l-
 \
 |
 done
?25h  Created wheel for peewee: filename=peewee-3.18.1-cp312-cp312-linux_x86_64.whl size=303841 sha256=baf58bd4cdf44f4ab84d3ffc8e5379f2524e28509400b725aa2812e1f8a01684
  Stored in directory: /home/runner/.cache/pip/wheels/1a/57/6a/bb71346381d0d911cd4ce3026f1fa720da76707e4f01cf27dd
Successfully built peewee
Installing collected packages: peewee, multitasking, websockets, frozendict, curl_cffi, yfinance
Successfully installed curl_cffi-0.11.3 frozendict-2.4.6 multitasking-0.0.11 peewee-3.18.1 websockets-15.0.1 yfinance-0.2.63

We will use the following imports

import numpy as np
import yfinance as yf
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.optimize import brentq
from scipy.stats import beta

28.2. Data#

The figure below shows the price of cotton in USD since the start of 2016.

Hide code cell source
s = yf.download('CT=F', '2016-1-1', '2023-4-1')['Close']
Hide code cell output
/tmp/ipykernel_7477/3037438160.py:1: FutureWarning: YF.download() has changed argument auto_adjust default to True
  s = yf.download('CT=F', '2016-1-1', '2023-4-1')['Close']
[*********************100%***********************]  1 of 1 completed

Hide code cell source
fig, ax = plt.subplots()

ax.plot(s, marker='o', alpha=0.5, ms=1)
ax.set_ylabel('cotton price in USD', fontsize=12)
ax.set_xlabel('date', fontsize=12)

plt.show()
_images/f17bfe7250554804c8970fa6d45f251f20fff56138e0c9292515587192cd1471.png

The figure shows surprisingly large movements in the price of cotton.

What causes these movements?

In general, prices depend on the choices and actions of

  1. suppliers,

  2. consumers, and

  3. speculators.

Our focus will be on the interaction between these parties.

We will connect them together in a dynamic model of supply and demand, called the competitive storage model.

This model was developed by [Samuelson, 1971], [Wright and Williams, 1982], [Scheinkman and Schechtman, 1983], [Deaton and Laroque, 1992], [Deaton and Laroque, 1996], and [Chambers and Bailey, 1996].

28.3. The competitive storage model#

In the competitive storage model, commodities are assets that

  1. can be traded by speculators and

  2. have intrinsic value to consumers.

Total demand is the sum of consumer demand and demand by speculators.

Supply is exogenous, depending on “harvests”.

Note

These days, goods such as basic computer chips and integrated circuits are often treated as commodities in financial markets, being highly standardized, and, for these kinds of commodities, the word “harvest” is not appropriate.

Nonetheless, we maintain it for simplicity.

The equilibrium price is determined competitively.

It is a function of the current state (which determines current harvests and predicts future harvests).

28.4. The model#

Consider a market for a single commodity, whose price is given at tt by ptp_t.

The harvest of the commodity at time tt is ZtZ_t.

We assume that the sequence {Zt}t1\{ Z_t \}_{t \geq 1} is IID with common density function ϕ\phi, where ϕ\phi is nonnegative.

Speculators can store the commodity between periods, with ItI_t units purchased in the current period yielding αIt\alpha I_t units in the next.

Here the parameter α(0,1)\alpha \in (0,1) is a depreciation rate for the commodity.

For simplicity, the risk free interest rate is taken to be zero, so expected profit on purchasing ItI_t units is

Etpt+1αItptIt=(αEtpt+1pt)It \mathbb{E}_t \, p_{t+1} \cdot \alpha I_t - p_t I_t = (\alpha \mathbb{E}_t \, p_{t+1} - p_t) I_t

Here Etpt+1\mathbb{E}_t \, p_{t+1} is the expectation of pt+1p_{t+1} taken at time tt.

28.5. Equilibrium#

In this section we define the equilibrium and discuss how to compute it.

28.5.1. Equilibrium conditions#

Speculators are assumed to be risk neutral, which means that they buy the commodity whenever expected profits are positive.

As a consequence, if expected profits are positive, then the market is not in equilibrium.

Hence, to be in equilibrium, prices must satisfy the “no-arbitrage” condition

(28.1)#αEtpt+1pt0 \alpha \mathbb{E}_t \, p_{t+1} - p_t \leq 0

This means that if the expected price is lower than the current price, there is no room for arbitrage.

Profit maximization gives the additional condition

(28.2)#αEtpt+1pt<0 implies It=0 \alpha \mathbb{E}_t \, p_{t+1} - p_t < 0 \text{ implies } I_t = 0

We also require that the market clears, with supply equaling demand in each period.

We assume that consumers generate demand quantity D(p)D(p) corresponding to price pp.

Let P:=D1P := D^{-1} be the inverse demand function.

Regarding quantities,

  • supply is the sum of carryover by speculators and the current harvest, and

  • demand is the sum of purchases by consumers and purchases by speculators.

Mathematically,

  • supply is given by Xt=αIt1+ZtX_t = \alpha I_{t-1} + Z_t, which takes values in S:=R+S := \mathbb R_+, while

  • demand =D(pt)+It = D(p_t) + I_t

Thus, the market equilibrium condition is

(28.3)#αIt1+Zt=D(pt)+It \alpha I_{t-1} + Z_t = D(p_t) + I_t

The initial condition X0SX_0 \in S is treated as given.

28.5.2. An equilibrium function#

How can we find an equilibrium?

Our path of attack will be to seek a system of prices that depend only on the current state.

(Our solution method involves using an ansatz, which is an educated guess — in this case for the price function.)

In other words, we take a function pp on SS and set pt=p(Xt)p_t = p(X_t) for every tt.

Prices and quantities then follow

(28.4)#pt=p(Xt),It=XtD(pt),Xt+1=αIt+Zt+1 p_t = p(X_t), \quad I_t = X_t - D(p_t), \quad X_{t+1} = \alpha I_t + Z_{t+1}

We choose pp so that these prices and quantities satisfy the equilibrium conditions above.

More precisely, we seek a pp such that (28.1) and (28.2) hold for the corresponding system (28.4).

(28.5)#p(x)=max{α0p(αI(x)+z)ϕ(z)dz,P(x)}(xS) p^*(x) = \max \left\{ \alpha \int_0^\infty p^*(\alpha I(x) + z) \phi(z)dz, P(x) \right\} \qquad (x \in S)

where

(28.6)#I(x):=xD(p(x))(xS) I(x) := x - D(p^*(x)) \qquad (x \in S)

It turns out that such a pp^* will suffice, in the sense that (28.1) and (28.2) hold for the corresponding system (28.4).

To see this, observe first that

Etpt+1=Etp(Xt+1)=Etp(αI(Xt)+Zt+1)=0p(αI(Xt)+z)ϕ(z)dz \mathbb{E}_t \, p_{t+1} = \mathbb{E}_t \, p^*(X_{t+1}) = \mathbb{E}_t \, p^*(\alpha I(X_t) + Z_{t+1}) = \int_0^\infty p^*(\alpha I(X_t) + z) \phi(z)dz

Thus (28.1) requires that

α0p(αI(Xt)+z)ϕ(z)dzp(Xt) \alpha \int_0^\infty p^*(\alpha I(X_t) + z) \phi(z)dz \leq p^*(X_t)

This inequality is immediate from (28.5).

Second, regarding (28.2), suppose that

α0p(αI(Xt)+z)ϕ(z)dz<p(Xt) \alpha \int_0^\infty p^*(\alpha I(X_t) + z) \phi(z)dz < p^*(X_t)

Then by (28.5) we have p(Xt)=P(Xt)p^*(X_t) = P(X_t)

But then D(p(Xt))=XtD(p^*(X_t)) = X_t and It=I(Xt)=0I_t = I(X_t) = 0.

As a consequence, both (28.1) and (28.2) hold.

We have found an equilibrium, which verifies the ansatz.

28.5.3. Computing the equilibrium#

We now know that an equilibrium can be obtained by finding a function pp^* that satisfies (28.5).

It can be shown that, under mild conditions there is exactly one function on SS satisfying (28.5).

Moreover, we can compute this function using successive approximation.

This means that we start with a guess of the function and then update it using (28.5).

This generates a sequence of functions p1,p2,p_1, p_2, \ldots

We continue until this process converges, in the sense that pkp_k and pk+1p_{k+1} are very close together.

Then we take the final pkp_k that we computed as our approximation of pp^*.

To implement our update step, it is helpful if we put (28.5) and (28.6) together.

This leads us to the update rule

(28.7)#pk+1(x)=max{α0pk(α(xD(pk+1(x)))+z)ϕ(z)dz,P(x)} p_{k+1}(x) = \max \left\{ \alpha \int_0^\infty p_k(\alpha ( x - D(p_{k+1}(x))) + z) \phi(z)dz, P(x) \right\}

In other words, we take pkp_k as given and, at each xx, solve for qq in

(28.8)#q=max{α0pk(α(xD(q))+z)ϕ(z)dz,P(x)} q = \max \left\{ \alpha \int_0^\infty p_k(\alpha ( x - D(q)) + z) \phi(z)dz, P(x) \right\}

Actually we can’t do this at every xx, so instead we do it on a grid of points x1,,xnx_1, \ldots, x_n.

Then we get the corresponding values q1,,qnq_1, \ldots, q_n.

Then we compute pk+1p_{k+1} as the linear interpolation of the values q1,,qnq_1, \ldots, q_n over the grid x1,,xnx_1, \ldots, x_n.

Then we repeat, seeking convergence.

28.6. Code#

The code below implements this iterative process, starting from p0=Pp_0 = P.

The distribution ϕ\phi is set to a shifted Beta distribution (although many other choices are possible).

The integral in (28.8) is computed via Monte Carlo.

α, a, c = 0.8, 1.0, 2.0
beta_a, beta_b = 5, 5
mc_draw_size = 250
gridsize = 150
grid_max = 35
grid = np.linspace(a, grid_max, gridsize)

beta_dist = beta(5, 5)
Z = a + beta_dist.rvs(mc_draw_size) * c    # Shock observations
D = P = lambda x: 1.0 / x
tol = 1e-4


def T(p_array):

    new_p = np.empty_like(p_array)

    # Interpolate to obtain p as a function.
    p = interp1d(grid,
                 p_array,
                 fill_value=(p_array[0], p_array[-1]),
                 bounds_error=False)

    # Update
    for i, x in enumerate(grid):

        h = lambda q: q - max(α * np.mean(p(α * (x - D(q)) + Z)), P(x))
        new_p[i] = brentq(h, 1e-8, 100)

    return new_p


fig, ax = plt.subplots()

price = P(grid)
ax.plot(grid, price, alpha=0.5, lw=1, label="inverse demand curve")
error = tol + 1
while error > tol:
    new_price = T(price)
    error = max(np.abs(new_price - price))
    price = new_price

ax.plot(grid, price, 'k-', alpha=0.5, lw=2, label=r'$p^*$')
ax.legend()
ax.set_xlabel('$x$')
ax.set_ylabel("prices")

plt.show()
_images/9489a183b76539a3771ea7e4617ab3c50b45eb8cc94e5eaa7b76fea022979abf.png

The figure above shows the inverse demand curve PP, which is also p0p_0, as well as our approximation of pp^*.

Once we have an approximation of pp^*, we can simulate a time series of prices.

# Turn the price array into a price function
p_star = interp1d(grid,
                  price,
                  fill_value=(price[0], price[-1]),
                  bounds_error=False)

def carry_over(x):
    return α * (x - D(p_star(x)))

def generate_cp_ts(init=1, n=50):
    X = np.empty(n)
    X[0] = init
    for t in range(n-1):
            Z = a + c * beta_dist.rvs()
            X[t+1] = carry_over(X[t]) + Z
    return p_star(X)

fig, ax = plt.subplots()
ax.plot(generate_cp_ts(), label="price")
ax.set_xlabel("time")
ax.legend()
plt.show()
_images/dc7a73212ad4f1295087540ce2189aeb6ec962f10936bcf10a0d43192fe2691d.png