Previous topic

Continuous State DP Continued

Name Resolution

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

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