inventory_dynamics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "24d8829f",
   "metadata": {},
   "source": [
    "# Inventory Dynamics\n",
    "\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a630451c",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Inventory Dynamics](#Inventory-Dynamics)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Sample Paths](#Sample-Paths)  \n",
    "  - [Marginal Distributions](#Marginal-Distributions)  \n",
    "  - [Exercises](#Exercises)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4e36f05",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "In this lecture we will study the time path of inventories for firms that\n",
    "follow so-called s-S inventory dynamics.\n",
    "\n",
    "Such firms\n",
    "\n",
    "1. wait until inventory falls below some level $ s $ and then  \n",
    "1. order sufficient quantities to bring their inventory back up to capacity $ S $.  \n",
    "\n",
    "\n",
    "These kinds of policies are common in practice and also optimal in certain circumstances.\n",
    "\n",
    "A review of early literature and some macroeconomic implications can be found in [[Caplin, 1985](https://python.quantecon.org/zreferences.html#id50)].\n",
    "\n",
    "Here our main aim is to learn more about simulation, time series and Markov dynamics.\n",
    "\n",
    "While our Markov environment and many of the concepts we consider are related to those found in our [lecture on finite Markov chains](https://python.quantecon.org/finite_markov.html), the state space is a continuum in the current application.\n",
    "\n",
    "Let’s start with some imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dfd39ee",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"figure.figsize\"] = (11, 5)  #set default figure size\n",
    "import numpy as np\n",
    "from numba import njit, float64, prange\n",
    "from numba.experimental import jitclass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c38b410",
   "metadata": {},
   "source": [
    "## Sample Paths\n",
    "\n",
    "Consider a firm with inventory $ X_t $.\n",
    "\n",
    "The firm waits until $ X_t \\leq s $ and then restocks up to $ S $ units.\n",
    "\n",
    "It faces stochastic demand $ \\{ D_t \\} $, which we assume is IID.\n",
    "\n",
    "With notation $ a^+ := \\max\\{a, 0\\} $, inventory dynamics can be written\n",
    "as\n",
    "\n",
    "$$\n",
    "X_{t+1} =\n",
    "    \\begin{cases}\n",
    "      ( S - D_{t+1})^+ & \\quad \\text{if } X_t \\leq s \\\\\n",
    "      ( X_t - D_{t+1} )^+ &  \\quad \\text{if } X_t > s\n",
    "    \\end{cases}\n",
    "$$\n",
    "\n",
    "In what follows, we will assume that each $ D_t $ is lognormal, so that\n",
    "\n",
    "$$\n",
    "D_t = \\exp(\\mu + \\sigma Z_t)\n",
    "$$\n",
    "\n",
    "where $ \\mu $ and $ \\sigma $ are parameters and $ \\{Z_t\\} $ is IID\n",
    "and standard normal.\n",
    "\n",
    "Here’s a class that stores parameters and generates time paths for inventory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcaa5822",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "firm_data = [\n",
    "   ('s', float64),          # restock trigger level\n",
    "   ('S', float64),          # capacity\n",
    "   ('mu', float64),         # shock location parameter\n",
    "   ('sigma', float64)       # shock scale parameter\n",
    "]\n",
    "\n",
    "\n",
    "@jitclass(firm_data)\n",
    "class Firm:\n",
    "\n",
    "    def __init__(self, s=10, S=100, mu=1.0, sigma=0.5):\n",
    "\n",
    "        self.s, self.S, self.mu, self.sigma = s, S, mu, sigma\n",
    "\n",
    "    def update(self, x):\n",
    "        \"Update the state from t to t+1 given current state x.\"\n",
    "\n",
    "        Z = np.random.randn()\n",
    "        D = np.exp(self.mu + self.sigma * Z)\n",
    "        if x <= self.s:\n",
    "            return max(self.S - D, 0)\n",
    "        else:\n",
    "            return max(x - D, 0)\n",
    "\n",
    "    def sim_inventory_path(self, x_init, sim_length):\n",
    "\n",
    "        X = np.empty(sim_length)\n",
    "        X[0] = x_init\n",
    "\n",
    "        for t in range(sim_length-1):\n",
    "            X[t+1] = self.update(X[t])\n",
    "        return X"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35c62894",
   "metadata": {},
   "source": [
    "Let’s run a first simulation, of a single path:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d6cab1f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "firm = Firm()\n",
    "\n",
    "s, S = firm.s, firm.S\n",
    "sim_length = 100\n",
    "x_init = 50\n",
    "\n",
    "X = firm.sim_inventory_path(x_init, sim_length)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "bbox = (0., 1.02, 1., .102)\n",
    "legend_args = {'ncol': 3,\n",
    "               'bbox_to_anchor': bbox,\n",
    "               'loc': 3,\n",
    "               'mode': 'expand'}\n",
    "\n",
    "ax.plot(X, label=\"inventory\")\n",
    "ax.plot(np.full(sim_length, s), 'k--', label=\"$s$\")\n",
    "ax.plot(np.full(sim_length, S), 'k-', label=\"$S$\")\n",
    "ax.set_ylim(0, S+10)\n",
    "ax.set_xlabel(\"time\")\n",
    "ax.legend(**legend_args)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a83626a",
   "metadata": {},
   "source": [
    "Now let’s simulate multiple paths in order to build a more complete picture of\n",
    "the probabilities of different outcomes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cec28a6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "sim_length=200\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(np.full(sim_length, s), 'k--', label=\"$s$\")\n",
    "ax.plot(np.full(sim_length, S), 'k-', label=\"$S$\")\n",
    "ax.set_ylim(0, S+10)\n",
    "ax.legend(**legend_args)\n",
    "\n",
    "for i in range(400):\n",
    "    X = firm.sim_inventory_path(x_init, sim_length)\n",
    "    ax.plot(X, 'b', alpha=0.2, lw=0.5)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "110c97da",
   "metadata": {},
   "source": [
    "## Marginal Distributions\n",
    "\n",
    "Now let’s look at the marginal distribution $ \\psi_T $ of $ X_T $ for some\n",
    "fixed $ T $.\n",
    "\n",
    "We will do this by generating many draws of $ X_T $ given initial\n",
    "condition $ X_0 $.\n",
    "\n",
    "With these draws of $ X_T $ we can build up a picture of its distribution $ \\psi_T $.\n",
    "\n",
    "Here’s one visualization, with $ T=50 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4b79b2f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "T = 50\n",
    "M = 200  # Number of draws\n",
    "\n",
    "ymin, ymax = 0, S + 10\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(11, 6))\n",
    "\n",
    "for ax in axes:\n",
    "    ax.grid(alpha=0.4)\n",
    "\n",
    "ax = axes[0]\n",
    "\n",
    "ax.set_ylim(ymin, ymax)\n",
    "ax.set_ylabel('$X_t$', fontsize=16)\n",
    "ax.vlines((T,), -1.5, 1.5)\n",
    "\n",
    "ax.set_xticks((T,))\n",
    "ax.set_xticklabels((r'$T$',))\n",
    "\n",
    "sample = np.empty(M)\n",
    "for m in range(M):\n",
    "    X = firm.sim_inventory_path(x_init, 2 * T)\n",
    "    ax.plot(X, 'b-', lw=1, alpha=0.5)\n",
    "    ax.plot((T,), (X[T+1],), 'ko', alpha=0.5)\n",
    "    sample[m] = X[T+1]\n",
    "\n",
    "axes[1].set_ylim(ymin, ymax)\n",
    "\n",
    "axes[1].hist(sample,\n",
    "             bins=16,\n",
    "             density=True,\n",
    "             orientation='horizontal',\n",
    "             histtype='bar',\n",
    "             alpha=0.5)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9856b6d",
   "metadata": {},
   "source": [
    "We can build up a clearer picture by drawing more samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ce17ef2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "T = 50\n",
    "M = 50_000\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "sample = np.empty(M)\n",
    "for m in range(M):\n",
    "    X = firm.sim_inventory_path(x_init, T+1)\n",
    "    sample[m] = X[T]\n",
    "\n",
    "ax.hist(sample,\n",
    "         bins=36,\n",
    "         density=True,\n",
    "         histtype='bar',\n",
    "         alpha=0.75)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7eb7e302",
   "metadata": {},
   "source": [
    "Note that the distribution is bimodal\n",
    "\n",
    "- Most firms have restocked twice but a few have restocked only once (see figure with paths above).  \n",
    "- Firms in the second category have lower inventory.  \n",
    "\n",
    "\n",
    "We can also approximate the distribution using a [kernel density estimator](https://en.wikipedia.org/wiki/Kernel_density_estimation).\n",
    "\n",
    "Kernel density estimators can be thought of as smoothed histograms.\n",
    "\n",
    "They are preferable to histograms when the distribution being estimated is likely to be smooth.\n",
    "\n",
    "We will use a kernel density estimator from [scikit-learn](https://scikit-learn.org/stable/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb2a4acf",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KernelDensity\n",
    "\n",
    "def plot_kde(sample, ax, label=''):\n",
    "\n",
    "    xmin, xmax = 0.9 * min(sample), 1.1 * max(sample)\n",
    "    xgrid = np.linspace(xmin, xmax, 200)\n",
    "    kde = KernelDensity(kernel='gaussian').fit(sample[:, None])\n",
    "    log_dens = kde.score_samples(xgrid[:, None])\n",
    "\n",
    "    ax.plot(xgrid, np.exp(log_dens), label=label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19f080fe",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plot_kde(sample, ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69e7fbf5",
   "metadata": {},
   "source": [
    "The allocation of probability mass is similar to what was shown by the\n",
    "histogram just above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be4c24f3",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d34608c8",
   "metadata": {},
   "source": [
    "## Exercise 20.1\n",
    "\n",
    "This model is asymptotically stationary, with a unique stationary\n",
    "distribution.\n",
    "\n",
    "(See the discussion of stationarity in [our lecture on AR(1) processes](https://intro.quantecon.org/ar1_processes.html) for background — the fundamental concepts are the same.)\n",
    "\n",
    "In particular, the sequence of marginal distributions $ \\{\\psi_t\\} $\n",
    "is converging to a unique limiting distribution that does not depend on\n",
    "initial conditions.\n",
    "\n",
    "Although we will not prove this here, we can investigate it using simulation.\n",
    "\n",
    "Your task is to generate and plot the sequence $ \\{\\psi_t\\} $ at times\n",
    "$ t = 10, 50, 250, 500, 750 $ based on the discussion above.\n",
    "\n",
    "(The kernel density estimator is probably the best way to present each\n",
    "distribution.)\n",
    "\n",
    "You should see convergence, in the sense that differences between successive distributions are getting smaller.\n",
    "\n",
    "Try different initial conditions to verify that, in the long run, the distribution is invariant across initial conditions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1167ac0",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 20.1](https://python.quantecon.org/#id_ex1)\n",
    "\n",
    "Below is one possible solution:\n",
    "\n",
    "The computations involve a lot of CPU cycles so we have tried to write the\n",
    "code efficiently.\n",
    "\n",
    "This meant writing a specialized function rather than using the class above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f35c104",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "s, S, mu, sigma = firm.s, firm.S, firm.mu, firm.sigma\n",
    "\n",
    "@njit(parallel=True)\n",
    "def shift_firms_forward(current_inventory_levels, num_periods):\n",
    "\n",
    "    num_firms = len(current_inventory_levels)\n",
    "    new_inventory_levels = np.empty(num_firms)\n",
    "\n",
    "    for f in prange(num_firms):\n",
    "        x = current_inventory_levels[f]\n",
    "        for t in range(num_periods):\n",
    "            Z = np.random.randn()\n",
    "            D = np.exp(mu + sigma * Z)\n",
    "            if x <= s:\n",
    "                x = max(S - D, 0)\n",
    "            else:\n",
    "                x = max(x - D, 0)\n",
    "        new_inventory_levels[f] = x\n",
    "\n",
    "    return new_inventory_levels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62c671d5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "x_init = 50\n",
    "num_firms = 50_000\n",
    "\n",
    "sample_dates = 0, 10, 50, 250, 500, 750\n",
    "\n",
    "first_diffs = np.diff(sample_dates)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "X = np.full(num_firms, x_init)\n",
    "\n",
    "current_date = 0\n",
    "for d in first_diffs:\n",
    "    X = shift_firms_forward(X, d)\n",
    "    current_date += d\n",
    "    plot_kde(X, ax, label=f't = {current_date}')\n",
    "\n",
    "ax.set_xlabel('inventory')\n",
    "ax.set_ylabel('probability')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "806e4ec6",
   "metadata": {},
   "source": [
    "Notice that by $ t=500 $ or $ t=750 $ the densities are barely\n",
    "changing.\n",
    "\n",
    "We have reached a reasonable approximation of the stationary density.\n",
    "\n",
    "You can convince yourself that initial conditions don’t matter by\n",
    "testing a few of them.\n",
    "\n",
    "For example, try rerunning the code above with all firms starting at\n",
    "$ X_0 = 20 $ or $ X_0 = 80 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c77c79d6",
   "metadata": {},
   "source": [
    "## Exercise 20.2\n",
    "\n",
    "Using simulation, calculate the probability that firms that start with\n",
    "$ X_0 = 70 $ need to order twice or more in the first 50 periods.\n",
    "\n",
    "You will need a large sample size to get an accurate reading."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e361f53",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 20.2](https://python.quantecon.org/#id_ex2)\n",
    "\n",
    "Here is one solution.\n",
    "\n",
    "Again, the computations are relatively intensive so we have written a a\n",
    "specialized function rather than using the class above.\n",
    "\n",
    "We will also use parallelization across firms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea5b9c8a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@njit(parallel=True)\n",
    "def compute_freq(sim_length=50, x_init=70, num_firms=1_000_000):\n",
    "\n",
    "    firm_counter = 0  # Records number of firms that restock 2x or more\n",
    "    for m in prange(num_firms):\n",
    "        x = x_init\n",
    "        restock_counter = 0  # Will record number of restocks for firm m\n",
    "\n",
    "        for t in range(sim_length):\n",
    "            Z = np.random.randn()\n",
    "            D = np.exp(mu + sigma * Z)\n",
    "            if x <= s:\n",
    "                x = max(S - D, 0)\n",
    "                restock_counter += 1\n",
    "            else:\n",
    "                x = max(x - D, 0)\n",
    "\n",
    "        if restock_counter > 1:\n",
    "            firm_counter += 1\n",
    "\n",
    "    return firm_counter / num_firms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b8cd49e",
   "metadata": {},
   "source": [
    "Note the time the routine takes to run, as well as the output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d2f69b4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "freq = compute_freq()\n",
    "print(f\"Frequency of at least two stock outs = {freq}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d348a23",
   "metadata": {},
   "source": [
    "Try switching the `parallel` flag to `False` in the jitted function\n",
    "above.\n",
    "\n",
    "Depending on your system, the difference can be substantial.\n",
    "\n",
    "(On our desktop machine, the speed up is by a factor of 5.)"
   ]
  }
 ],
 "metadata": {
  "date": 1714442505.877428,
  "filename": "inventory_dynamics.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Inventory Dynamics"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}