{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2cf3a74c",
   "metadata": {},
   "source": [
    "# Commodity Prices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30d0740e",
   "metadata": {},
   "source": [
    "## Outline\n",
    "\n",
    "For more than half of all countries around the globe, [commodities](https://en.wikipedia.org/wiki/Commodity) account for [the majority of total exports](https://unctad.org/publication/commodities-and-development-report-2019).\n",
    "\n",
    "Examples of commodities include copper, diamonds, iron ore, lithium, cotton\n",
    "and coffee beans.\n",
    "\n",
    "In this lecture we give an introduction to the theory of commodity prices.\n",
    "\n",
    "The lecture is quite advanced relative to other lectures in this series.\n",
    "\n",
    "We need to compute an equilibrium, and that equilibrium is described by a\n",
    "price function.\n",
    "\n",
    "We will solve an equation where the price function is the unknown.\n",
    "\n",
    "This is harder than solving an equation for an unknown number, or vector.\n",
    "\n",
    "The lecture will discuss one way to solve a [functional equation](https://en.wikipedia.org/wiki/Functional_equation) (an equation where the unknown object is a function).\n",
    "\n",
    "For this lecture we need the `yfinance` library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66749a36",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "!pip install yfinance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1440818",
   "metadata": {},
   "source": [
    "We will use the following imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02a9dcd5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import yfinance as yf\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import interp1d\n",
    "from scipy.optimize import brentq\n",
    "from scipy.stats import beta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7db3f0c7",
   "metadata": {},
   "source": [
    "## Data\n",
    "\n",
    "The figure below shows the price of cotton in USD since the start of 2016."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd4113c1",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "s = yf.download('CT=F', '2016-1-1', '2023-4-1')['Close']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d87ed9e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(s, marker='o', alpha=0.5, ms=1)\n",
    "ax.set_ylabel('cotton price in USD', fontsize=12)\n",
    "ax.set_xlabel('date', fontsize=12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "926a9588",
   "metadata": {},
   "source": [
    "The figure shows surprisingly large movements in the price of cotton.\n",
    "\n",
    "What causes these movements?\n",
    "\n",
    "In general, prices depend on the choices and actions of\n",
    "\n",
    "1. suppliers,  \n",
    "1. consumers, and  \n",
    "1. speculators.  \n",
    "\n",
    "\n",
    "Our focus will be on the interaction between these parties.\n",
    "\n",
    "We will connect them together in a dynamic model of supply and demand, called\n",
    "the *competitive storage model*.\n",
    "\n",
    "This model was developed by\n",
    "[[Samuelson, 1971](https://intro.quantecon.org/zreferences.html#id22)],\n",
    "[[Wright and Williams, 1982](https://intro.quantecon.org/zreferences.html#id21)], [[Scheinkman and Schechtman, 1983](https://intro.quantecon.org/zreferences.html#id20)],\n",
    "[[Deaton and Laroque, 1992](https://intro.quantecon.org/zreferences.html#id19)], [[Deaton and Laroque, 1996](https://intro.quantecon.org/zreferences.html#id18)], and\n",
    "[[Chambers and Bailey, 1996](https://intro.quantecon.org/zreferences.html#id17)]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d85e8a19",
   "metadata": {},
   "source": [
    "## The competitive storage model\n",
    "\n",
    "In the competitive storage model, commodities are assets that\n",
    "\n",
    "1. can be traded by speculators and  \n",
    "1. have intrinsic value to consumers.  \n",
    "\n",
    "\n",
    "Total demand is the sum of consumer demand and demand by speculators.\n",
    "\n",
    "Supply is exogenous, depending on “harvests”.\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">These days, goods such as basic computer chips and integrated circuits are\n",
    "often treated as commodities in financial markets, being highly standardized,\n",
    "and, for these kinds of commodities, the word “harvest” is not\n",
    "appropriate.\n",
    "\n",
    "Nonetheless, we maintain it for simplicity.\n",
    "\n",
    "The equilibrium price is determined competitively.\n",
    "\n",
    "It is a function of the current state (which determines\n",
    "current harvests and predicts future harvests)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db271a49",
   "metadata": {},
   "source": [
    "## The model\n",
    "\n",
    "Consider a market for a single commodity, whose price is given at $ t $ by\n",
    "$ p_t $.\n",
    "\n",
    "The harvest of the commodity at time $ t $ is $ Z_t $.\n",
    "\n",
    "We assume that the sequence $ \\{ Z_t \\}_{t \\geq 1} $ is IID with common density function $ \\phi $, where $ \\phi $ is nonnegative.\n",
    "\n",
    "Speculators can store the commodity between periods, with $ I_t $ units\n",
    "purchased in the current period yielding $ \\alpha I_t $ units in the next.\n",
    "\n",
    "Here the parameter $ \\alpha \\in (0,1) $ is a depreciation rate for the commodity.\n",
    "\n",
    "For simplicity, the risk free interest rate is taken to be\n",
    "zero, so expected profit on purchasing $ I_t $ units is\n",
    "\n",
    "$$\n",
    "\\mathbb{E}_t \\, p_{t+1} \\cdot \\alpha I_t - p_t I_t\n",
    "   = (\\alpha \\mathbb{E}_t \\, p_{t+1} - p_t) I_t\n",
    "$$\n",
    "\n",
    "Here $ \\mathbb{E}_t \\, p_{t+1} $ is the expectation of $ p_{t+1} $ taken at time\n",
    "$ t $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bc11bd9",
   "metadata": {},
   "source": [
    "## Equilibrium\n",
    "\n",
    "In this section we define the equilibrium and discuss how to compute it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c862fc07",
   "metadata": {},
   "source": [
    "### Equilibrium conditions\n",
    "\n",
    "Speculators are assumed to be risk neutral, which means that they buy the\n",
    "commodity whenever expected profits are positive.\n",
    "\n",
    "As a consequence, if expected profits are positive, then the market is not in\n",
    "equilibrium.\n",
    "\n",
    "Hence, to be in equilibrium, prices must satisfy the “no-arbitrage”\n",
    "condition\n",
    "\n",
    "\n",
    "<a id='equation-eq-arbi'></a>\n",
    "$$\n",
    "\\alpha \\mathbb{E}_t \\, p_{t+1}  - p_t \\leq 0 \\tag{28.1}\n",
    "$$\n",
    "\n",
    "This means that if the expected price is lower than the current price, there is no room for arbitrage.\n",
    "\n",
    "Profit maximization gives the additional condition\n",
    "\n",
    "\n",
    "<a id='equation-eq-pmco'></a>\n",
    "$$\n",
    "\\alpha \\mathbb{E}_t \\, p_{t+1}  - p_t  < 0 \\text{ implies } I_t = 0 \\tag{28.2}\n",
    "$$\n",
    "\n",
    "We also require that the market clears, with supply equaling demand in each period.\n",
    "\n",
    "We assume that consumers generate demand quantity $ D(p) $ corresponding to\n",
    "price $ p $.\n",
    "\n",
    "Let $ P := D^{-1} $ be the inverse demand function.\n",
    "\n",
    "Regarding quantities,\n",
    "\n",
    "- supply is the sum of carryover by speculators and the current harvest, and  \n",
    "- demand is the sum of purchases by consumers and purchases by speculators.  \n",
    "\n",
    "\n",
    "Mathematically,\n",
    "\n",
    "- supply is given by $ X_t = \\alpha I_{t-1} + Z_t $, which takes values in $ S := \\mathbb R_+ $, while  \n",
    "- demand $ = D(p_t) + I_t $  \n",
    "\n",
    "\n",
    "Thus, the market equilibrium condition is\n",
    "\n",
    "\n",
    "<a id='equation-eq-mkeq'></a>\n",
    "$$\n",
    "\\alpha I_{t-1} + Z_t =  D(p_t) + I_t \\tag{28.3}\n",
    "$$\n",
    "\n",
    "The initial condition $ X_0 \\in S $ is treated as given."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c88b012",
   "metadata": {},
   "source": [
    "### An equilibrium function\n",
    "\n",
    "How can we find an equilibrium?\n",
    "\n",
    "Our path of attack will be to seek a system of prices that depend only on the\n",
    "current state.\n",
    "\n",
    "(Our solution method involves using an [ansatz](https://en.wikipedia.org/wiki/Ansatz), which is an educated guess — in this case for the price function.)\n",
    "\n",
    "In other words, we take a function $ p $ on $ S $ and set $ p_t = p(X_t) $ for every $ t $.\n",
    "\n",
    "Prices and quantities then follow\n",
    "\n",
    "\n",
    "<a id='equation-eq-eosy'></a>\n",
    "$$\n",
    "p_t = p(X_t), \\quad I_t = X_t - D(p_t), \\quad X_{t+1} = \\alpha I_t + Z_{t+1} \\tag{28.4}\n",
    "$$\n",
    "\n",
    "We choose $ p $ so that these prices and quantities satisfy the equilibrium\n",
    "conditions above.\n",
    "\n",
    "More precisely, we seek a $ p $ such that [(28.1)](#equation-eq-arbi) and [(28.2)](#equation-eq-pmco) hold for\n",
    "the corresponding system [(28.4)](#equation-eq-eosy).\n",
    "\n",
    "\n",
    "<a id='equation-eq-dopf'></a>\n",
    "$$\n",
    "p^*(x) = \\max\n",
    "    \\left\\{\n",
    "    \\alpha \\int_0^\\infty p^*(\\alpha I(x) + z) \\phi(z)dz, P(x)\n",
    "    \\right\\}\n",
    "    \\qquad (x \\in S) \\tag{28.5}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "\n",
    "<a id='equation-eq-einvf'></a>\n",
    "$$\n",
    "I(x) := x - D(p^*(x))\n",
    "    \\qquad (x \\in S) \\tag{28.6}\n",
    "$$\n",
    "\n",
    "It turns out that such a $ p^* $ will suffice, in the sense that [(28.1)](#equation-eq-arbi)\n",
    "and [(28.2)](#equation-eq-pmco) hold for the corresponding system [(28.4)](#equation-eq-eosy).\n",
    "\n",
    "To see this, observe first that\n",
    "\n",
    "$$\n",
    "\\mathbb{E}_t \\, p_{t+1}\n",
    "   = \\mathbb{E}_t \\, p^*(X_{t+1})\n",
    "   = \\mathbb{E}_t \\, p^*(\\alpha I(X_t) + Z_{t+1})\n",
    "   = \\int_0^\\infty p^*(\\alpha I(X_t) + z) \\phi(z)dz\n",
    "$$\n",
    "\n",
    "Thus [(28.1)](#equation-eq-arbi) requires that\n",
    "\n",
    "$$\n",
    "\\alpha \\int_0^\\infty p^*(\\alpha I(X_t) + z) \\phi(z)dz \\leq p^*(X_t)\n",
    "$$\n",
    "\n",
    "This inequality is immediate from [(28.5)](#equation-eq-dopf).\n",
    "\n",
    "Second, regarding [(28.2)](#equation-eq-pmco), suppose that\n",
    "\n",
    "$$\n",
    "\\alpha \\int_0^\\infty p^*(\\alpha I(X_t) + z) \\phi(z)dz < p^*(X_t)\n",
    "$$\n",
    "\n",
    "Then by [(28.5)](#equation-eq-dopf) we have $ p^*(X_t) = P(X_t) $\n",
    "\n",
    "But then $ D(p^*(X_t)) = X_t $ and $ I_t = I(X_t) = 0 $.\n",
    "\n",
    "As a consequence, both [(28.1)](#equation-eq-arbi) and [(28.2)](#equation-eq-pmco) hold.\n",
    "\n",
    "We have found an equilibrium, which verifies the ansatz."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25658f16",
   "metadata": {},
   "source": [
    "### Computing the equilibrium\n",
    "\n",
    "We now know that an equilibrium can be obtained by finding a function $ p^* $\n",
    "that satisfies [(28.5)](#equation-eq-dopf).\n",
    "\n",
    "It can be shown that, under mild conditions there is exactly one function on\n",
    "$ S $ satisfying [(28.5)](#equation-eq-dopf).\n",
    "\n",
    "Moreover, we can compute this function using successive approximation.\n",
    "\n",
    "This means that we start with a guess of the function and then update it using\n",
    "[(28.5)](#equation-eq-dopf).\n",
    "\n",
    "This generates a sequence of functions $ p_1, p_2, \\ldots $\n",
    "\n",
    "We continue until this process converges, in the sense that $ p_k $ and\n",
    "$ p_{k+1} $ are very close together.\n",
    "\n",
    "Then we take the final $ p_k $ that we computed as our approximation of $ p^* $.\n",
    "\n",
    "To implement our update step, it is helpful if we put [(28.5)](#equation-eq-dopf) and\n",
    "[(28.6)](#equation-eq-einvf) together.\n",
    "\n",
    "This leads us to the update rule\n",
    "\n",
    "\n",
    "<a id='equation-eq-dopf2'></a>\n",
    "$$\n",
    "p_{k+1}(x) = \\max\n",
    "    \\left\\{\n",
    "    \\alpha \\int_0^\\infty p_k(\\alpha ( x - D(p_{k+1}(x))) + z) \\phi(z)dz, P(x)\n",
    "    \\right\\} \\tag{28.7}\n",
    "$$\n",
    "\n",
    "In other words, we take $ p_k $ as given and, at each $ x $, solve for $ q $ in\n",
    "\n",
    "\n",
    "<a id='equation-eq-dopf3'></a>\n",
    "$$\n",
    "q = \\max\n",
    "    \\left\\{\n",
    "    \\alpha \\int_0^\\infty p_k(\\alpha ( x - D(q)) + z) \\phi(z)dz, P(x)\n",
    "    \\right\\} \\tag{28.8}\n",
    "$$\n",
    "\n",
    "Actually we can’t do this at every $ x $, so instead we do it on a grid of\n",
    "points $ x_1, \\ldots, x_n $.\n",
    "\n",
    "Then we get the corresponding values $ q_1, \\ldots, q_n $.\n",
    "\n",
    "Then we compute $ p_{k+1} $ as the linear interpolation of\n",
    "the values $ q_1, \\ldots, q_n $ over the grid $ x_1, \\ldots, x_n $.\n",
    "\n",
    "Then we repeat, seeking convergence."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "405226eb",
   "metadata": {},
   "source": [
    "## Code\n",
    "\n",
    "The code below implements this iterative process, starting from $ p_0 = P $.\n",
    "\n",
    "The distribution $ \\phi $ is set to a shifted Beta distribution (although many\n",
    "other choices are possible).\n",
    "\n",
    "The integral in [(28.8)](#equation-eq-dopf3) is computed via [Monte Carlo](https://intro.quantecon.org/monte_carlo.html#monte-carlo)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8610c2d2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "α, a, c = 0.8, 1.0, 2.0\n",
    "beta_a, beta_b = 5, 5\n",
    "mc_draw_size = 250\n",
    "gridsize = 150\n",
    "grid_max = 35\n",
    "grid = np.linspace(a, grid_max, gridsize)\n",
    "\n",
    "beta_dist = beta(5, 5)\n",
    "Z = a + beta_dist.rvs(mc_draw_size) * c    # Shock observations\n",
    "D = P = lambda x: 1.0 / x\n",
    "tol = 1e-4\n",
    "\n",
    "\n",
    "def T(p_array):\n",
    "\n",
    "    new_p = np.empty_like(p_array)\n",
    "\n",
    "    # Interpolate to obtain p as a function.\n",
    "    p = interp1d(grid,\n",
    "                 p_array,\n",
    "                 fill_value=(p_array[0], p_array[-1]),\n",
    "                 bounds_error=False)\n",
    "\n",
    "    # Update\n",
    "    for i, x in enumerate(grid):\n",
    "\n",
    "        h = lambda q: q - max(α * np.mean(p(α * (x - D(q)) + Z)), P(x))\n",
    "        new_p[i] = brentq(h, 1e-8, 100)\n",
    "\n",
    "    return new_p\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "price = P(grid)\n",
    "ax.plot(grid, price, alpha=0.5, lw=1, label=\"inverse demand curve\")\n",
    "error = tol + 1\n",
    "while error > tol:\n",
    "    new_price = T(price)\n",
    "    error = max(np.abs(new_price - price))\n",
    "    price = new_price\n",
    "\n",
    "ax.plot(grid, price, 'k-', alpha=0.5, lw=2, label=r'$p^*$')\n",
    "ax.legend()\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel(\"prices\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fc505ca",
   "metadata": {},
   "source": [
    "The figure above shows the inverse demand curve $ P $, which is also $ p_0 $, as\n",
    "well as our approximation of $ p^* $.\n",
    "\n",
    "Once we have an approximation of $ p^* $, we can simulate a time series of\n",
    "prices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a5c4b31",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Turn the price array into a price function\n",
    "p_star = interp1d(grid,\n",
    "                  price,\n",
    "                  fill_value=(price[0], price[-1]),\n",
    "                  bounds_error=False)\n",
    "\n",
    "def carry_over(x):\n",
    "    return α * (x - D(p_star(x)))\n",
    "\n",
    "def generate_cp_ts(init=1, n=50):\n",
    "    X = np.empty(n)\n",
    "    X[0] = init\n",
    "    for t in range(n-1):\n",
    "            Z = a + c * beta_dist.rvs()\n",
    "            X[t+1] = carry_over(X[t]) + Z\n",
    "    return p_star(X)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(generate_cp_ts(), label=\"price\")\n",
    "ax.set_xlabel(\"time\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "date": 1750037599.0890553,
  "filename": "commod_price.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Commodity Prices"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}