DEMSLV05

Cournot Equilibrium Model

In [1]:
import numpy as np
import matplotlib.pyplot as plt

from compecon import NLP, gridmake
from compecon.demos import demo
%matplotlib inline
C:\ProgramData\Anaconda3\lib\site-packages\ggplot\utils.py:81: FutureWarning: pandas.tslib is deprecated and will be removed in a future version.
You can access Timestamp as pandas.Timestamp
  pd.tslib.Timestamp,
C:\ProgramData\Anaconda3\lib\site-packages\ggplot\stats\smoothers.py:4: FutureWarning: The pandas.lib module is deprecated and will be removed in a future version. These are private functions and can be accessed from pandas._libs.lib instead
  from pandas.lib import Timestamp
C:\ProgramData\Anaconda3\lib\site-packages\statsmodels\compat\pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.
  from pandas.core import datetools

Parameters and initial value

In [2]:
alpha = 0.625
beta = np.array([0.6, 0.8])

Set up the Cournot function

In [3]:
def market(q):
    quantity = q.sum()
    price = quantity ** (-alpha)
    return price, quantity
In [4]:
def cournot(q):
    P, Q = market(q)
    P1 = -alpha * P/Q
    P2 = (-alpha - 1) * P1 / Q
    fval = P + (P1 - beta) * q
    fjac = np.diag(2 * P1 + P2 * q - beta) + np.fliplr(np.diag(P1 + P2 * q))
    return fval, fjac

Compute equilibrium using Newton method (explicitly)

In [5]:
q = np.array([0.2, 0.2])

for it in range(40):
    f, J = cournot(q)
    step = -np.linalg.solve(J, f)
    q += step
    if np.linalg.norm(step) < 1.e-10: break

price, quantity = market(q)
print(f'Company 1 produces {q[0]:.4f} units, while company 2 produces {q[1]:.4f} units.')
print(f'Total production is {quantity:.4f} and price is {price:.4f}')
Company 1 produces 0.8396 units, while company 2 produces 0.6888 units.
Total production is 1.5284 and price is 0.7671

Using a NLP object

In [11]:
q0 = np.array([0.2, 0.2])
cournot_problem = NLP(cournot)
q = cournot_problem.newton(q0, print=True)

price, quantity = market(q)
print(f'\nCompany 1 produces {q[0]:.4f} units, while company 2 produces {q[1]:.4f} units.')
print(f'Total production is {quantity:.4f} and price is {price:.4f}')
Solving nonlinear equations by Newton's method
it    bstep  change
--------------------
   0     0  4.64e-01
   1     0  9.53e-02
   2     0  3.47e-03
   3     0  4.20e-06
   4     0  5.77e-12

Company 1 produces 0.8396 units, while company 2 produces 0.6888 units.
Total production is 1.5284 and price is 0.7671

Generate data for contour plot

In [7]:
n = 100
q1 = np.linspace(0.1, 1.5, n)
q2 = np.linspace(0.1, 1.5, n)
z = np.array([cournot(q)[0] for q in gridmake(q1, q2).T]).T

Plot figures

In [8]:
steps_options = {'marker': 'o',
                 'color': (0.2, 0.2, .81),
                 'linewidth': 1.0,
                 'markersize': 9,
                 'markerfacecolor': 'white',
                 'markeredgecolor': 'red'}

contour_options = {'levels': [0.0],
                   'colors': '0.25',
                   'linewidths': 0.5}


Q1, Q2 = np.meshgrid(q1, q2)
Z0 = np.reshape(z[0], (n,n), order='F')
Z1 = np.reshape(z[1], (n,n), order='F')

methods = ['newton', 'broyden']
cournot_problem.opts['maxit', 'maxsteps', 'all_x'] = 10, 0, True

qmin, qmax = 0.1, 1.3

plt.figure(figsize=[12,6])
for it, method in enumerate(methods):
    x = cournot_problem.zero(method=method)
    demo.subplot(1, 2, it + 1, methods[it].capitalize() + "'s method",
                 '$q_1$', '$q_2$', [qmin, qmax], [qmin, qmax])
    plt.contour(Q1, Q2, Z0, **contour_options)
    plt.contour(Q1, Q2, Z1, **contour_options)
    plt.plot(*cournot_problem.x_sequence, **steps_options)

    demo.text(0.85, qmax, '$\pi_1 = 0$', 'left', 'top')
    demo.text(qmax, 0.55, '$\pi_2 = 0$', 'right', 'center')
plt.show()
Solving nonlinear equations by Newton's method
it    bstep  change
--------------------
   0     0  1.10e+00
   1     0  4.64e-01
   2     0  9.53e-02
   3     0  3.47e-03
   4     0  4.20e-06
Solving nonlinear equations by Broyden's method
it    bstep  change
--------------------
   0     0  4.64e-01
   1     0  2.17e-01
   2     0  5.45e-02
   3     0  1.53e-02
   4     0  1.04e-02
   5     0  4.55e-03
   6     0  1.33e-04
   7     0  2.94e-07
   8     0  1.34e-09
In [9]:
steps_options = {'marker': 'o',
                 'color': (0.2, 0.2, .81),
                 'linewidth': 2.5,
                 'markersize': 9,
                 'markerfacecolor': 'white',
                 'markeredgecolor': 'red'}

contour_options = {'levels': [0.0],
                   'colors': 'white',
                   'linewidths': 2.0}


Q1, Q2 = np.meshgrid(q1, q2)
Z0 = np.reshape(z[0], (n,n), order='F')
Z1 = np.reshape(z[1], (n,n), order='F')

methods = ['newton', 'broyden']
cournot_problem.opts['maxit', 'maxsteps', 'all_x'] = 10, 0, True

qmin, qmax = 0.1, 1.3

plt.figure(figsize=[12,12])

x = cournot_problem.zero(method='newton')
demo.figure("Convergence of Newton's method", '$q_1$', '$q_2$', [qmin, qmax], [qmin, qmax])
plt.contour(Q1, Q2, Z0, **contour_options)
plt.contour(Q1, Q2, Z1, **contour_options)

plt.plot(*cournot_problem.x_sequence, **steps_options)

demo.text(0.85, qmax, '$\pi_1 = 0$', 'left', 'top')
demo.text(qmax, 0.55, '$\pi_2 = 0$', 'right', 'center')

plt.savefig('c:/users/rromero/latex-aux/cournot.pdf',bbox_inches='tight')
plt.show()
Solving nonlinear equations by Newton's method
it    bstep  change
--------------------
   0     0  1.10e+00
   1     0  4.64e-01
   2     0  9.53e-02
   3     0  3.47e-03
   4     0  4.20e-06
<matplotlib.figure.Figure at 0x19d54143160>