This article originally appeared in Spanish as [SECMCA Regional Economic Notes No. 93](http://www.secmca.org/NOTAS_ECONOMICAS/articulo93NOV2017.pdf) in December 2017. This Jupyter notebook is a slightly modified version of the original, and was prepared on July 28, 2018.
The use of programming languages has become increasingly necessary in many of the tasks that we economist do, either because the theoretic models we develop lack analytic solutions, or because the new econometric estimator we wish to use is not yet available in a program with a graphical user interface (GUI, where we could just click with a mouse), or just because it is not efficient to analyze huge datasets with spreadsheets.
Firms profits are \begin{align*} \pi_1\left(q_1, q_2\right) &=P\left(q_1+q_2\right)q_1 - C_1\left(q_1\right) \\ \pi_2\left(q_1, q_2\right) &=P\left(q_1+q_2\right)q_2 - C_2\left(q_2\right) \end{align*}
In a Cournot equilibrium, each firm maximizes its profits taking as given the other firm's output. Thus, it must follow that \begin{align*} \frac{\partial \pi_1\left(q_1, q_2\right)}{\partial q_1} &= P\left(q_1+q_2\right) + P'\left(q_1+q_2\right)q_1 - C'_1\left(q_1\right) = 0\\ \frac{\partial \pi_2\left(q_1, q_2\right)}{\partial q_2} &= P\left(q_1+q_2\right) + P'\left(q_1+q_2\right)q_2 - C'_2\left(q_2\right) = 0 \end{align*}
Therefore, equilibrium output levels for this market are given by the solution to this nonlinear equation system \begin{equation} \label{eq:fzero} f\left(q_1, q_2\right) = \begin{bmatrix} \left(q_1+q_2\right)^{-\alpha} - \alpha q_1 \left(q_1+q_2\right)^{-\alpha - 1} - \beta_1q_1 \\ \left(q_1+q_2\right)^{-\alpha} - \alpha q_2 \left(q_1+q_2\right)^{-\alpha - 1} - \beta_2q_2\end{bmatrix} = \begin{bmatrix}0 \\ 0\end{bmatrix}\tag{1} \end{equation}
To find the root of the function defined in (1) we will use Newton's method. In general, this method is applied to the function $f: \Re^n \to \Re^n$ to find some (Notice that, depending on the function, there could be more than one solution, or no solution at all.) value $x^*$ such that $f(x^*)=0$. To that end, we start with a value $x_0 \in \Re^n$ and make the recursion \begin{equation}\label{eq:newton} x_{i+1} = x_i - J^{-1}(x_i) f(x_i) \tag{2} \end{equation} where $J(x_i)$ corresponds to the Jacobian of $f$ evaluated at $x_0$. In theory, following this recursion $x_i$ converges to $x^*$ as long as the $f$ function is continuously differentiable and the initial value $x_0$ is “sufficiently close” to the $x^*$ root.
First, we begin a Python session and import compecon
. This package is available in pip, and can be installed by typing
pip install compecon
in a command window.
import numpy as np
import matplotlib.pyplot as plt
from compecon import NLP, gridmake
from compecon.demos import demo
To solve this model computationally, we need to assign values to the parameters, so we set $\alpha=0.625$, $\beta_1=0.6$ and $\beta_2=0.8$.
alpha = 0.625
beta = np.array([0.6, 0.8])
The unknowns in our problem are the firms' output levels, $q_1$ and $q_2$. We define the market
function to tell us total output and resulting price, given the levels of $q_1$ y $q_2$. Notice that both quantities are passed to this function in the q
vector
def market(q):
quantity = q.sum()
price = quantity ** (-alpha)
return price, quantity
Then, we define the cournot
function, returning a two-element tuple: the objective function and its Jacobian matrix, both evaluated in a pair of quantities contained by the q
vector. To make the code easier, notice that the (1) function can be written more succinctly as
\begin{equation*}
f\left(q_1, q_2\right) = \begin{bmatrix}
P + \left(P' - c_1\right)q_1 \\
P + \left(P' - c_2\right)q_2\end{bmatrix}
= \begin{bmatrix}0 \\ 0\end{bmatrix}
\end{equation*}
and its Jacobian matrix is
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
Next, we compute the equilibrium using Newton's method (equation (2)) to find the root of the cournot
function. We set $q_0 = \begin{bmatrix}0.2 & 0.2\end{bmatrix}'$ as our initial value and iterate until the norm of the change between two successive values of the recursion is less than $10^{-10}$.
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'\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}')
Company 1 produces 0.8396 units, while company 2 produces 0.6888 units. Total production is 1.5284 and price is 0.7671
After just five iterations, Newton's method converges to the solution, which Python prints to screen:
Company 1 produces 0.8396 units, while company 2 produces 0.6888 units.
Total production is 1.5284 and price is 0.7671
q0 = np.array([0.2, 0.2])
cournot_problem = NLP(cournot)
q = cournot_problem.newton(q0)
price, quantity = market(q)
print(f'\nCompany 1 produces {q[0]:.4f} units, while' +
f' 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
After completing this code block, Python prints the following to screen:
Company 1 produces 0.8396 units, while company 2 produces 0.6888 units.
Total production is 1.5284 and price is 0.7671
Figure 1 illustrates the problem we just solved, where the axes represent the output levels of each firm. The quasi-vertical white line represents the profit-maximizing output level for firm 1, taking the output of firm 2 as given. Similarly, the quasi-horizontal line represents the profit maximizing output level for firm 2, given firm 1 output. The solution to the problem corresponds to the intersection of these two lines. See also the path to convergence (blue line) from the initial $q_0 = \begin{bmatrix}0.2 & 0.2\end{bmatrix}'$ point to the solution.
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
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
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')
where $D(p)$ is the market demand curve.
for every price $p>0$. In simple cases, this function can be found explicitly. However, in more complicated cases, there is no explicit solution. Suppose for example that demand and marginal cost are given by \begin{equation*} D(p) = p^{-\eta} \qquad\qquad CM(q) = \alpha\sqrt{q} + q^2 \end{equation*}
so that the functional equation to be solved for $S(p)$ is \begin{equation} \label{eq:funcional} \left[p - \frac{S(p)p^{\eta+1}}{\eta}\right] - \left[\alpha\sqrt{S(p)} + S(p)^2\right] = 0 \tag{5} \end{equation}
In equation (5), the unknown is the supply function $S(p)$, which makes (5) and infinite-dimension equation. Instead of solving the equation directly, we will approximate its solution using $n$ Chebyshev polynomials $\phi_i(x)$, which are defined recursively for $x \in [0,1]$ as: \begin{align*} \phi_0(x) & = 1 \\ \phi_1(x) & = x \\ \phi_{k + 1}(p_i) & = 2x \phi_k(x) - \phi_{k-1}(x), \qquad \text{for} k = 1,2, \dots \end{align*}
Substituting this last expression in (5) for each of the placement nodes (Chebyshev in this case) results in a non-linear system of $ n $ equations (one for each node) in $ n $ unknowns $ c_k $ (one for each polynomial of Cheybshev), which in principle can be solved by Newton's method, as in the last example. Thus, in practice, the system to be solved is
\begin{equation} \label{eq:collocation} \left[p_i - \frac{\left(\sum_{k=0}^{n-1}c_{k}\phi_k(p_i)\right)p_i^{\eta+1}}{\eta}\right] - \left[\alpha\sqrt{\sum_{k=0}^{n-1}c_{k}\phi_k(p_i)} + \left(\sum_{k=0}^{n-1}c_{k}\phi_k(p_i)\right)^2\right] = 0 \tag{7} \end{equation}for $i=1,2,\dots,n$ and for $k=1,2,\dots,n$.
To solve this model we start a new Python session:
import numpy as np
import matplotlib.pyplot as plt
from compecon import BasisChebyshev, NLP, nodeunif
from compecon.demos import demo
and set the $\alpha$ and $\eta$ parameters
alpha, eta = 1.0, 3.5
For convenience, we define a lambda
function to represent the demand
D = lambda p: p**(-eta)
We will approximate the solution for prices in the $p\in [\frac{1}{2}, 2]$ interval, using 25 collocation nodes. The compecon
library provides the BasisChebyshev
class to make computations with Chebyshev bases:
n, a, b = 25, 0.5, 2.0
S = BasisChebyshev(n, a, b, labels=['price'], l=['supply'])
Let's assume that our first guess is $S(p)=1$. To that end, we set the value of S
to one in each of the nodes
p = S.nodes
S.y = np.ones_like(p)
It is important to highlight that in this problem the unknowns are the $c_k$ coefficients from the Chebyshev basis; however, an object of BasisChebyshev
class automatically adjusts those coefficients so they are consistent with the values we set for the function at the nodes (here indicated by the .y
property).
def resid(c):
S.c = c # update interpolation coefficients
q = S(p) # compute quantity supplied at price nodes
return p - q * (p ** (eta+1) / eta) - alpha * np.sqrt(q) - q ** 2
Note that the resid
function takes a single argument (the coefficients for the Chebyshev basis). All other parameters (Q, p, eta, alpha
must be declared in the main script, where Python will find their values.
cournot = NLP(resid)
S.c = cournot.broyden(S.c, tol=1e-12, print=True)
Solving nonlinear equations by Broyden's method it bstep change -------------------- 0 0 4.08e-01 1 0 8.95e-02 2 0 1.37e-02 3 0 2.01e-03 4 0 3.36e-04 5 0 8.11e-05 6 0 1.28e-05 7 0 2.80e-06 8 0 4.88e-07 9 0 7.55e-08 10 0 1.84e-08 11 0 1.88e-09 12 0 4.63e-10 13 0 5.88e-11 14 0 1.17e-11 15 0 2.48e-12 16 0 2.91e-13
After 17 iterations, Broyden's method converges to the desired solution. We can visualize this in Figure 3, which shows the value of the function on 501 different points within the approximation interval. Notice that the residual plot crosses the horizontal axis 25 times; this occurs precisely at the collocation nodes (represented by red dots). This figure also shows the precision of the approximation: outside nodes, the function is within $5\times10^{-11}$ units from zero.
pplot = nodeunif(501, a, b)
demo.figure('Cournot Effective Firm Supply Function',
'Quantity', 'Price', [0, 4], [a, b])
plt.plot(5 * S(pplot), pplot, D(pplot), pplot)
plt.legend(('Supply','Demand'))
<matplotlib.legend.Legend at 0x1d4d5ffad30>
This block generates Figure 3.
p = pplot
demo.figure('Residual Function for Cournot Problem',
'Quantity', 'Residual')
plt.hlines(0, a, b, 'k', '--', lw=2)
plt.plot(pplot, resid(S.c))
plt.plot(S.nodes,np.zeros_like(S.nodes),'r*');
We now plot the effective supply for a varying number of firms; the result is shown in Figure 4.
m = np.array([1, 3, 5, 10, 15, 20])
demo.figure('Supply and Demand Functions', 'Quantity', 'Price', [0, 13])
plt.plot(np.outer(S(pplot), m), pplot)
plt.plot(D(pplot), pplot, linewidth=4, color='black')
plt.legend(['m=1', 'm=3', 'm=5', 'm=10', 'm=15', 'm=20', 'demand']);
In Figure 4 notice how the equilibrium price and quantity change as the number of firms increases.
The last figure in this example (Figure 5), shows the equilibrium price as a function of the number of firms.
pp = (b + a) / 2
dp = (b - a) / 2
m = np.arange(1, 26)
for i in range(50):
dp /= 2
pp = pp - np.sign(S(pp) * m - D(pp)) * dp
demo.figure('Cournot Equilibrium Price as Function of Industry Size',
'Number of Firms', 'Price')
plt.bar(m, pp);
from pandas_datareader import wb
To be able to download data from the World Bank, we first need to know the exact code of the indicator we want to read. The first time we do this task we will not know this code, but we can look for it in the World Bank website or more easily from Python itself. For example, to find data on real GDP per capita, we run the following using the .search
function:
wb.search('gdp.*capita.*const').iloc[:,:2]
id | name | |
---|---|---|
646 | 6.0.GDPpc_constant | GDP per capita, PPP (constant 2011 internation... |
7858 | NY.GDP.PCAP.KD | GDP per capita (constant 2010 US$) |
7860 | NY.GDP.PCAP.KN | GDP per capita (constant LCU) |
7862 | NY.GDP.PCAP.PP.KD | GDP per capita, PPP (constant 2011 internation... |
7863 | NY.GDP.PCAP.PP.KD.87 | GDP per capita, PPP (constant 1987 internation... |
where the dot followed by an asterisk (.*) indicates that any text in that position is a match. This function returns a data table with information about indicators that match the search criteria. In the preceding line, we use the code .iloc[:,:2]
so that Python only prints the first two columns from that table.
paises = ['CR', 'DO', 'GT', 'HN', 'NI', 'SV']
and we proceed to reed data from 1991:
datos = wb.download(indicator='NY.GDP.PCAP.KD',
country=paises,start=1991, end=2016)
It is also possible to read data for more than one indicator in a single call to the wb.download
function, writing their codes in a list (just like we did to read data on all six countries at once). In any case, we get a data table in panel format, where each columns corresponds to one of the indicators. For our example in particular, where we only read one indicator, it would be useful if the table was arranged so that each row correspond to a year and each column to a country. We can achieve it with this instruction:
GDP = datos.reset_index().pivot('year','country')
Once data is arrange this way, it is very easy to compute growth for all countries in a single step:
GROWTH = 100 * GDP.pct_change()
or to generate a formatted data table to be included in a \LaTeX document
#GROWTH.tail(6).round(2).to_latex('micuadro.tex')
print(GROWTH.tail(6).round(2).to_latex())
\begin{tabular}{lrrrrrr} \toprule {} & \multicolumn{6}{l}{NY.GDP.PCAP.KD} \\ country & Costa Rica & Dominican Republic & El Salvador & Guatemala & Honduras & Nicaragua \\ year & & & & & & \\ \midrule 2011 & 3.06 & 1.81 & 3.34 & 1.94 & 1.89 & 5.03 \\ 2012 & 3.59 & 1.42 & 2.34 & 0.80 & 2.24 & 5.24 \\ 2013 & 1.13 & 3.59 & 1.89 & 1.54 & 0.99 & 3.72 \\ 2014 & 2.40 & 6.35 & 1.49 & 2.03 & 1.29 & 3.60 \\ 2015 & 2.55 & 5.79 & 1.87 & 2.03 & 2.08 & 3.60 \\ 2016 & 3.10 & 5.41 & 2.06 & 1.04 & 2.02 & 3.50 \\ \bottomrule \end{tabular}
In last instruction, the .tail(6)
part indicates that we only want the last six observations, while the .to_latex('micuadro.tex')
part exports that table to a file named 'micuadro.tex', which can later te included in a document. The result of this code will look similar to this:
GROWTH.tail(6).round(2)
NY.GDP.PCAP.KD | ||||||
---|---|---|---|---|---|---|
country | Costa Rica | Dominican Republic | El Salvador | Guatemala | Honduras | Nicaragua |
year | ||||||
2011 | 3.06 | 1.81 | 3.34 | 1.94 | 1.89 | 5.03 |
2012 | 3.59 | 1.42 | 2.34 | 0.80 | 2.24 | 5.24 |
2013 | 1.13 | 3.59 | 1.89 | 1.54 | 0.99 | 3.72 |
2014 | 2.40 | 6.35 | 1.49 | 2.03 | 1.29 | 3.60 |
2015 | 2.55 | 5.79 | 1.87 | 2.03 | 2.08 | 3.60 |
2016 | 3.10 | 5.41 | 2.06 | 1.04 | 2.02 | 3.50 |
Finally, we plot the results in Figure 6. It is possible to enhance the aesthetic aspect of this figure by, for example, changing the legend position. Such enhancements are not presented here for space considerations.
GROWTH.columns = paises
GROWTH.plot();
It is also possible to plot each of the time series in a separate subplot, with the instruction
GROWTH.plot(subplots=True, layout=[2,3], sharey=True);
where we have specified that each time series should be plotted separately (subplots=True
), be arranged in two rows and three columns (layout=[2,3]
), and all subplots must have the same “y” axis (sharey=True
, to facilitate country comparisons).
where $c_t$ stands for consumption, $y_t$ income, $\epsilon$ a stochastic shock. In this case $\beta_1$ corresponds to the income elasticity of consumption.
import pandas_datareader.data as web
from numpy import log
import statsmodels.formula.api as smf
Once this is done, we are ready to import data. In this example, we use quarterly data on consumption and production in the United States, available in FRED, a database from the Federal Reserve Bank of Saint Louis. For “consumption” we use the “PCEC” (Personal Consumption Expenditures) series, and for “income” we use “GDP” (Gross Domestic Product).
usdata = web.DataReader(['PCEC','GDP'],'fred', 1947, 2017)
After executing this instuction, the usdata
variable points to a pandas
data table, in which each column corresponds to a variable and each row to a quarter. We now estimate the model by ordinary least squares (.ols
) and print a summary of the results
mod = smf.ols('PCEC ~ GDP', log(usdata)).fit()
print(mod.summary())
OLS Regression Results ============================================================================== Dep. Variable: PCEC R-squared: 1.000 Model: OLS Adj. R-squared: 1.000 Method: Least Squares F-statistic: 6.154e+05 Date: Sat, 28 Jul 2018 Prob (F-statistic): 0.00 Time: 18:27:34 Log-Likelihood: 584.98 No. Observations: 281 AIC: -1166. Df Residuals: 279 BIC: -1159. Df Model: 1 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ Intercept -0.6712 0.010 -64.053 0.000 -0.692 -0.651 GDP 1.0268 0.001 784.505 0.000 1.024 1.029 ============================================================================== Omnibus: 51.339 Durbin-Watson: 0.075 Prob(Omnibus): 0.000 Jarque-Bera (JB): 85.582 Skew: 1.027 Prob(JB): 2.61e-19 Kurtosis: 4.758 Cond. No. 47.1 ============================================================================== Warnings: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Notice that the .ols
function takes two arguments, the formula specifying the model, and the name of the data table containing the variables. In this code block we specify the data as log(usdata)
, which tells Python that we want the logarithm of the data, saving us the task of generating another data table with the transformed data beforehand (as would be necessary in, for example, Stata).
with we estimate in Python with
smf.ols('PCEC ~ GDP', log(usdata).diff()).fit().summary()
Dep. Variable: | PCEC | R-squared: | 0.490 |
---|---|---|---|
Model: | OLS | Adj. R-squared: | 0.488 |
Method: | Least Squares | F-statistic: | 267.2 |
Date: | Sat, 28 Jul 2018 | Prob (F-statistic): | 1.50e-42 |
Time: | 18:27:34 | Log-Likelihood: | 996.74 |
No. Observations: | 280 | AIC: | -1989. |
Df Residuals: | 278 | BIC: | -1982. |
Df Model: | 1 | ||
Covariance Type: | nonrobust |
coef | std err | t | P>|t| | [0.025 | 0.975] | |
---|---|---|---|---|---|---|
Intercept | 0.0063 | 0.001 | 8.747 | 0.000 | 0.005 | 0.008 |
GDP | 0.6130 | 0.038 | 16.346 | 0.000 | 0.539 | 0.687 |
Omnibus: | 102.218 | Durbin-Watson: | 2.550 |
---|---|---|---|
Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 1087.804 |
Skew: | -1.140 | Prob(JB): | 6.11e-237 |
Kurtosis: | 12.383 | Cond. No. | 90.9 |
We notice that now the $R^2$ is no longer close to one, and that the Durbin-Watson statistic is closer to 2.0, indicating lack of autocorrelation.
usdata['lPCEC'] = log(usdata['PCEC'])
usdata['lGDP'] = log(usdata['GDP'])
usdata['dlPCEC'] = usdata['lPCEC'].diff()
usdata['dlGDP'] = usdata['lGDP'].diff()
model = smf.ols('dlPCEC ~ dlGDP', usdata)
results = model.fit()
table = results.summary()
print(table)
With results from last Table at hand, we could predict that a one percentage point (p.p.) increase in GDP growth would lead to a 0.618 p.p. increase in consumption growth. However, given that the data sample covers such a long period (nearly 70 years of quarterly observation), it is reasonable to wonder whether the parameters in this model are constant, given that several structural changes could have occurred along these years. One way to evaluate such posibility is to estimate the model with a rolling sample. In particular, we are going to estimate this model with 24 quarterly observations rolling window, changing the sample by one quarter in every step.
In this case, since we are going to need growth data many times, it is more efficient to compute growth data only once and store it in a growth
variable. With the [1:]
code we are dropping the first observation, which we lose when we compute the first-order difference (.diff
). Furthermore, we use the .shape
property from the table to find out how many observations T
we have, and then we set the window range to h=24
observations:
growth = (100*log(usdata).diff())[1:]
T, nvar = growth.shape
h = 24
To faciliate next step, we define función window_beta1
, which takes as its only argument the number of the last observation to be included in the estimation, and returns the value of the estimated GDP coefficient
def window_beta1(k):
return smf.ols('PCEC~GDP',growth[k-h:k]).fit().params['GDP']
With this, we are ready to estimate the model many times, adding the results to the growth
table as the beta1
“indicator”. Plotting the results we get Figure 8, where we clearly see that the effect of GDP growth on consumption growth is quite unstable, and thus the predictions made with the simple model could be very poor.
growth.loc[h-1:,'beta1'] = [window_beta1(k) for k in range(h,T+1)]
growth[['beta1']].plot();