Solving the competitive storage model using collocation

Dec. 26, 2015

The original Python notebook is on Github. Download it if you would like to follow along.

This IPython notebook solves the competitive storage model for the price of storable commodities using collocation as the numerical method. We reproduce an example from section 9.9.2 of Miranda and Fackler’s Applied Computational Economics and Finance using Python and independently of the authors’ CompEcon MATLAB toolbox.

The problem

The competitive storage model is a simple model for the price of a storable commodity. It has an uneven empirical record but it provides a good illustration of the computational problems which arise in modern macroecnomics.

We consider a commodity, say wheat, which can be stored subject to the condition that a fraction \(\lambda\) of the total stock at the end of a period disappears before the beginning of the next period and a cost per unit cost \(c\) has to be paid for storage. In each period there is an addition to stocks due to a stochastic harvest \(z\) and a reduction in stock from an outside demand given by a demand function \(P(q)\).

Storage decisions are taken by competitive speculators who must form expectations about future prices. These future prices depend on both the stochastic harvest in the next period as well as the total stock carried forward by all speculators taken together. We assume that speculators have rational expectations.The speculators are assumed to be risk neutral and discount their future profits by a discount factor \(\delta\).

Let \(s_t\) be the stock at the beginning of period \(t\) and \(x_t\) be the stock at the end of period \(t\) and \(p_t\) be the price that prevails in period \(t\). Then the equilibrium conditions are, first that outside consumers should demand the quantity \((s_t-x_t)\) that is not stored,

\[p_t = P(s_t-x_t)\]

that the stocks in consecutive periods be linked according to the storage technology, \[s_{t+1} = (1-\lambda)x_t+z_{t+1}\]

and most importantly the speculators maximize profits. \[p_t \ge \delta(1-\lambda)[E_t p_{t+1}]-c,\quad \text{with equality if $x_t>0$}\]

The last equation has the complementary slackness form it has becuase speculators may be at a corner solution where they do not carry forward any stocks at all. If this is the case then current prices can rise above the expected future return from carrying stocks. On the other hand current prices can never be lower than expected future revenues in equilibrium since such a condition would produce an infinite demand for wheat by speculators in the current period.

The solution algorithm

It is known that the above system of equations has a solution of the form \(p_t = \phi(s_t)\). However there is generally no analytical solution and the function \(\phi(\cdot)\) must be computed numerically.

Suppose we expect that prices next period are governed by the function \(f(s_{t+1})\). Then the price function this period, say \(g(s_{t})\), will be a function of this expectation through our equilibrium conditions above. We express this relationship in the form

\[g = Tf\]

where \(T\) is an abstract operator representing the working of the equilibrium conditions. It is precisely the property of a stationary rational equations solution that it is a fixed-point of this operator

\[\phi = T\phi\]

i.e., the expected law of prices is such as to produce behaviour which confirms the expect law of prices.

In order to find an approximate solution to this fixed-point equation we restrict our search to the finite dimensional space of Chebyshev polynomials of a given degree. The true solution may not lie in this space and thus no function in this space may satisfy the fixed-point equation at all points. In order to nonetheless get a good approximation to the true solution the collocation method we use in this notebook requires that the function come close to satising the fixed-point equation at a given finite set of points. In this notebook we take these points to be Chebyshev nodes. The task of finding the Chebyshev polynomial which satisfies this condition is carried out through successive approximation: starting from an initial guess \(f_0\) we calculate \(y_n = (Tf_0)(x_n)\) at our chosen Chebyshev nodes \(x_n\). Then \(f_1\) is then taken to be a Chebyshev polynomial that passes through the points \((x_n,y_n)\). This is continued until we reach a fixed point.

There are two more technical difficulties. The operator \(T\) involves computing a conditional expectation in order to find expected profits. For arbitrary functions this cannot be done exactly and we have to make a numerical approximation. In this notebook we use Gauss-Hermite quadrature to replace the continuous distribution of the shocks \(z_t\) by a discrete approximation.

Finally, we need to work with a bounded state space for our calculations. If \(\lambda>0\) and \(\bar z\) the maximum possible realization of the shock then the beginning-of-period stock can never exceed \(\bar z/\lambda\). This however does not work when \(\lambda=0\). Also using this bound is harmful even when \(\lambda>0\) since this forces us to consider very large values of stock which are not reached in equilibrium, thereby wasting computational resources. We therefore follow the literature in arbitrarily assuming that there is an upper abound \(x_{max}\) on end-of-period stocks. The value of this parameter is adjusted to be large enought that it is never reached in equilibrium.

See Judd, Numerical Methods in Economics, Chapter 11 and 17 for more details as well as alternative algorithms. Section 17.4 discusses the competitive storage model.



Collocation in general

Collocation using Chebyshev polynomials and nodes

The competitive storage model

The problem data structure


# Solving the problem
class StorageSolution(object):
    A solution to the competitve storage model.
        hgdeg: degree of Gauss-Hermite 
               discretization of shocks
        cdeg: dgree of Chebychev polynomial to use
        tol: tolerance for stopping criteria
        maxiter: maximum no of iterations 
        log: if True print progress report
        The model is solved in the process of constructing 
        the object, so constructor calls are expensive.
        If the model cannot be solved in the give number of 
        iterations throws a ValueError. Computational 
        problems can cause an AssertionError to be thrown.
    Instance variables:
    Saved from constructor: prob, hgdeg, cdeg, tol
        smax       Maximum beginning-of-period stock
        smin       Minimum beginning-of-period stock
        zs         An array of discretized harvest shock values
        ws         An array of probability weights 
                   corresponding to `ws`
        phi        The computed price as a function of 
                   begining-of-period shocks
                    Discounted expected revenue net of storage 
                    costs per unit stored when the pricing 
                    function is `phi` and end-of-period stocks
                    in the current period is `x`.
                    The optimal amount to store when the pricing 
                    function is `phi` and beginning-of-period 
                    stocks is `s`. Rational expectations is
                    built in: the expected stocks tomorrow is 
                    consistent with storage decision taken.

        T(phi)      We want T(phi,s) = phi(s) for all s.
    def __init__(self,p,ghdeg=5,cdeg=50,
        # Store input parameters
        self.prob = p
        self.ghdeg = ghdeg
        self.cdeg = cdeg
        self.tol = tol
        # Compute Gauss-Hermition integration nodes
        ys, = hermgauss(ghdeg) =
        self.zs = np.exp(ys*p.yvol)
        zmin,zmax = np.min(self.zs), np.max(self.zs)
        # Bounds for beginning-of-period stocks
        # Note xmax is maximum end-of-period shocks
        self.smin = zmin
        self.smax = zmax + (1 - p.shrink)*p.xmax
        # Storage revenue computation
        def storage_revenue(phi,x):
            Etom = sum(
                        phi(z + (1 - p.shrink)*x)*w 
                        for (z,w) in zip(self.zs,
            return p.disc*(1 - p.shrink)*Etom - p.cost
        self.storage_revenue = storage_revenue
        # Optimal storage decision
        def opti_storage(phi,s):
            P = p.P
            xmax = p.xmax
            if s>xmax and storage_revenue(phi,xmax) > P(s-xmax):
                return xmax
            elif storage_revenue(phi,0) < P(s-0):
                return 0
                max_feas = min(xmax,s)
                x_star,r = \
                    spopt.brentq(lambda x:\
                assert r.converged
                return x_star
        self.opti_storage = opti_storage
        # Initial guess for pricing function
        def init(s):
            return p.P(s)
        # We want T(phi) = phi
        def T(phi,s):
            x_star = opti_storage(phi,s)
            ans = p.P(s-x_star)
            assert ans>=0
            return ans
        self.T = T
        self.phi = chebyshev_solve(T,init,cdeg,
        if self.phi is None:
            raise ValueError("Failed to converge")

Visualization and simulation

Example from Miranda & Fackler

Iteration 1: Error = 0.14
Iteration 11: Error = 1.6e-08
Proportion of stockouts = 0.71

<matplotlib.text.Text at 0x7f475ed685c0>