Quantitative Economics With Python
Quantitative Economics With Python
Python
I Introduction to Python 1
1 About Python 3
3 An Introductory Example 35
4 Python Essentials 55
6 NumPy 81
7 Matplotlib 99
8 SciPy 111
9 Numba 123
15 Debugging 237
16 Pandas 245
3
4 CONTENTS
48 Robustness 813
Introduction to Python
1
Chapter 1
About Python
1.1 Contents
β’ Overview 1.2
1.2 Overview
At this stage, itβs not our intention that you try to replicate all you see.
We will work through what follows at a slow pace later in the lecture series.
Our only objective for this lecture is to give you some feel of what Python is, and what it can
do.
3
4 CHAPTER 1. ABOUT PYTHON
β’ communications
β’ web development
β’ CGI and graphical user interfaces
β’ games
β’ multimedia, data processing, security, etc., etc., etc.
β’ Google
β’ Dropbox
β’ Reddit
β’ YouTube
β’ Walt Disney Animation, etc., etc.
The following chart, produced using Stack Overflow Trends, shows one measure of the relative
popularity of Python
The figure indicates not only that Python is widely used but also that adoption of Python
has accelerated significantly since 2012.
We suspect this is driven at least in part by uptake in the scientific domain, particularly in
rapidly growing fields like data science.
1.3. WHATβS PYTHON? 5
For example, the popularity of pandas, a library for data analysis with Python has exploded,
as seen here.
(The corresponding time path for MATLAB is shown for comparison)
Note that pandas takes off in 2012, which is the same year that we seek Pythonβs popularity
begin to spike in the first figure.
Overall, itβs clear that
1.3.3 Features
One nice feature of Python is its elegant syntax β weβll see many examples later on.
Elegant code might sound superfluous but in fact itβs highly beneficial because it makes the
syntax easy to read and easy to remember.
Remembering how to read from files, sort dictionaries and other such routine tasks means
that you donβt need to break your flow in order to hunt down correct syntax.
6 CHAPTER 1. ABOUT PYTHON
Fundamental matrix and array processing capabilities are provided by the excellent NumPy
library.
NumPy provides the basic array data type plus some simple processing operations.
For example, letβs build some arrays
import numpy as np # Load the library
[1]:
a = np.linspace(-np.pi, np.pi, 100) # Create even grid from -Ο to Ο
b = np.cos(a) # Apply cosine to each element of a
c = np.sin(a) # Apply sin to each element of a
The number you see here might vary slightly but itβs essentially zero.
(For older versions of Python and NumPy you need to use the np.dot function)
The SciPy library is built on top of NumPy and provides additional functionality.
2
For example, letβs calculate β«β2 π(π§)ππ§ where π is the standard normal density.
1.4. SCIENTIFIC PROGRAMMING 7
οΏ½ = norm()
value, error = quad(οΏ½.pdf, -2, 2) # Integrate using Gaussian quadrature
value
0.9544997361036417
[3]:
β’ linear algebra
β’ integration
β’ interpolation
β’ optimization
β’ distributions and random number generation
β’ signal processing
β’ etc., etc.
1.4.2 Graphics
The most popular and comprehensive Python library for creating figures and graphs is Mat-
plotlib.
Example 3D plot
β’ Plotly
β’ Bokeh
1.4. SCIENTIFIC PROGRAMMING 9
[4]:
3π₯ + π¦
We can manipulate expressions
expression = (x + y)**2
[5]: expression.expand()
[5]:
π₯2 + 2π₯π¦ + π¦2
solve polynomials
from sympy import solve
[6]:
solve(x**2 + x + 2)
[7]:
β
limit(sin(x) / x, x, 0)
[8]:
[8]:
1
diff(sin(x), x)
[9]:
[9]:
cos (π₯)
The beauty of importing this functionality into Python is that we are working within a fully
fledged programming language.
Can easily create tables of derivatives, generate LaTeX output, add it to figures, etc., etc.
1.4.4 Statistics
Pythonβs data manipulation and statistics libraries have improved rapidly over the last few
years.
Pandas
One of the most popular libraries for working with data is pandas.
Pandas is fast, efficient, flexible and well designed.
10 CHAPTER 1. ABOUT PYTHON
price weight
2010-12-28 0.471435 -1.190976
2010-12-29 1.432707 -0.312652
2010-12-30 -0.720589 0.887163
2010-12-31 0.859588 -0.636524
2011-01-01 0.015696 -2.242685
df.mean()
[11]:
price 0.411768
[11]: weight -0.699135
dtype: float64
Hereβs some example code that generates and plots a random graph, with node color deter-
mined by shortest path length from a central node.
import networkx as nx
[12]: import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(1234)
p = nx.single_source_shortest_path_length(g, ncenter)
plt.figure()
nx.draw_networkx_edges(g, pos, alpha=0.4)
nx.draw_networkx_nodes(g,
pos,
nodelist=list(p.keys()),
node_size=120, alpha=0.5,
node_color=list(p.values()),
cmap=plt.cm.jet_r)
plt.show()
/home/ubuntu/anaconda3/lib/python3.7/site-
packages/networkx/drawing/nx_pylab.py:579: MatplotlibDeprecationWarning:
The iterable function was deprecated in Matplotlib 3.1 and will be removed in
3.3. Use np.iterable instead.
if not cb.iterable(width):
Running your Python code on massive servers in the cloud is becoming easier and easier.
A nice example is Anaconda Enterprise.
See also
- Amazon Elastic Compute Cloud
- The Google App Engine (Python, Java, PHP or Go)
- Pythonanywhere
- Sagemath Cloud
12 CHAPTER 1. ABOUT PYTHON
Apart from the cloud computing options listed above, you might like to consider
- Parallel computing through IPython clusters.
- The Starcluster interface to Amazonβs EC2.
- GPU programming through PyCuda, PyOpenCL, Theano or similar.
There are many other interesting developments with scientific programming in Python.
Some representative examples include
- Jupyter β Python in your browser with code cells, embedded images, etc.
- Numba β Make Python run at the same speed as native machine code!
- Blaze β a generalization of NumPy.
- PyTables β manage large data sets.
- CVXPY β convex optimization in Python.
2.1 Contents
β’ Overview 2.2
β’ Anaconda 2.3
β’ Exercises 2.8
2.2 Overview
1. get a Python environment up and running with all the necessary tools
2. execute simple Python commands
3. run a sample program
4. install the code libraries that underpin these lectures
2.3 Anaconda
The core Python package is easy to install but not what you should choose for these lectures.
These lectures require the entire scientific programming ecosystem, which
13
14 CHAPTER 2. SETTING UP YOUR PYTHON ENVIRONMENT
Hence the best approach for our purposes is to install a free Python distribution that contains
β’ very popular
β’ cross platform
β’ comprehensive
β’ completely unrelated to the Nicki Minaj song of the same name
Anaconda also comes with a great package management system to organize your code li-
braries.
All of what follows assumes that you adopt this recommendation!.
Installing Anaconda is straightforward: download the binary and follow the instructions.
Important points:
Anaconda supplies a tool called conda to manage and upgrade your Anaconda packages.
One conda command you should execute regularly is the one that updates the whole Ana-
conda distribution.
As a practice run, please execute the following
1. Open up a terminal
2. Type conda update anaconda
Jupyter notebooks are one of the many possible ways to interact with Python and the scien-
tific libraries.
They use a browser-based interface to Python with
2.4. JUPYTER NOTEBOOKS 15
Because of these possibilities, Jupyter is fast turning into a major player in the scientific com-
puting ecosystem.
Hereβs an image showing execution of some code (borrowed from here) in a Jupyter notebook
You can find a nice example of the kinds of things you can do in a Jupyter notebook (such as
include maths and text) here.
While Jupyter isnβt the only way to code in Python, itβs great for when you wish to
Once you have installed Anaconda, you can start the Jupyter notebook.
Either
If you use the second option, you will see something like this (click to enlarge)
Thus, the Jupyter kernel is listening for Python commands on port 8888 of our local machine.
Hopefully, your default browser has also opened up with a web page that looks something like
this (click to enlarge)
2.4. JUPYTER NOTEBOOKS 17
The notebook displays an active cell, into which you can type Python commands.
Letβs start with how to edit code and run simple programs.
Running Cells
Notice that in the previous figure the cell is surrounded by a green border.
This means that the cell is in edit mode.
As a result, you can type in Python code and it will appear in the cell.
When youβre ready to execute the code in a cell, hit Shift-Enter instead of the usual En-
ter.
2.4. JUPYTER NOTEBOOKS 19
(Note: There are also menu and button options for running code in a cell that you can find
by exploring)
Modal Editing
The next thing to understand about the Jupyter notebook is that it uses a modal editing sys-
tem.
This means that the effect of typing at the keyboard depends on which mode you are in.
The two modes are
1. Edit mode
1. Command mode
To switch to
β’ command mode from edit mode, hit the Esc key or Ctrl-M
20 CHAPTER 2. SETTING UP YOUR PYTHON ENVIRONMENT
The modal behavior of the Jupyter notebook is a little tricky at first but very efficient when
you get used to it.
User Interface Tour
At this stage, we recommend you take your time to
β’ look at the various options in the menus and see what they do
β’ take the βuser interface tourβ, which can be accessed through the help menu
N = 20
ΞΈ = np.linspace(0.0, 2 * np.pi, N, endpoint=False)
radii = 10 * np.random.rand(N)
width = np.pi / 4 * np.random.rand(N)
ax = plt.subplot(111, polar=True)
bars = ax.bar(ΞΈ, radii, width=width, bottom=0.0)
plt.show()
2.4. JUPYTER NOTEBOOKS 21
Donβt worry about the details for now β letβs just run it and see what happens.
The easiest way to run this code is to copy and paste into a cell in the notebook.
(In older versions of Jupyter you might need to add the command %matplotlib inline
before you generate the figure)
Clicking on the top right of the lower split closes the on-line help.
Other Content
In addition to executing code, the Jupyter notebook allows you to embed text, equations, fig-
ures and even videos in the page.
For example, here we enter a mixture of plain text and LaTeX instead of code
24 CHAPTER 2. SETTING UP YOUR PYTHON ENVIRONMENT
Next we Esc to enter command mode and then type m to indicate that we are writing Mark-
down, a mark-up language similar to (but simpler than) LaTeX.
(You can also use your mouse to select Markdown from the Code drop-down box just below
the list of menu items)
Now we Shift+Enter to produce this
2.4. JUPYTER NOTEBOOKS 25
Notebook files are just text files structured in JSON and typically ending with .ipynb.
You can share them in the usual way that you share files β or by using web services such as
nbviewer.
The notebooks you see on that site are static html representations.
To run one, download it as an ipynb file by clicking on the download icon at the top right.
Save it somewhere, navigate to it from the Jupyter dashboard and then run as discussed
above.
QuantEcon has its own site for sharing Jupyter notebooks related to economics β QuantEcon
Notes.
Notebooks submitted to QuantEcon Notes can be shared with a link, and are open to com-
ments and votes by the community.
26 CHAPTER 2. SETTING UP YOUR PYTHON ENVIRONMENT
into a cell.
Alternatively, you can type the following into a terminal
Using the run command is often easier than copy and paste.
(You might find that the % is unnecessary β use %automagic to toggle the need for %)
Note that Jupyter only looks for test.py in the present working directory (PWD).
If test.py isnβt in that directory, you will get an error.
Letβs look at a successful example, where we run a file test.py with contents:
for i in range(5):
[2]: print('foobar')
foobar
foobar
foobar
foobar
foobar
Here
β’ pwd asks Jupyter to show the PWD (or %pwd β see the comment about automagic
above)
β Note that test.py is there (on our computer, because we saved it there earlier).
β’ cat test.py asks Jupyter to print the contents of test.py (or !type test.py on
Windows)
If youβre trying to run a file not in the present working directory, youβll get an error.
To fix this error you need to either
One way to achieve the first option is to use the Upload button
β’ The button is on the top level dashboard, where Jupyter first opened to
β’ Look where the pointer is in this picture
Note: You can type the first letter or two of each directory name and then use the tab key to
expand.
2.7. EDITORS AND IDES 29
Itβs often convenient to be able to see your code before you run it.
In the following example, we execute load white_noise_plot.py where
white_noise_plot.py is in the PWD.
(Use %load if automagic is off)
Now the code from the file appears in a cell ready to execute.
The preceding discussion covers most of what you need to know to interact with this website.
30 CHAPTER 2. SETTING UP YOUR PYTHON ENVIRONMENT
However, as you start to write longer programs, you might want to experiment with your
workflow.
There are many different options and we mention them only in passing.
2.7.1 JupyterLab
A text editor is an application that is specifically designed to work with text files β such as
Python programs.
Nothing beats the power and efficiency of a good text editor for working with program text.
A good text editor will provide
β’ efficient text editing commands (e.g., copy, paste, search and replace)
β’ syntax highlighting, etc.
The IPython shell has many of the features of the notebook: tab completion, color syntax,
etc.
It also has command history through the arrow key.
The up arrow key to brings previously typed commands to the prompt.
This saves a lot of typingβ¦
Hereβs one set up, on a Linux box, with
2.7.4 IDEs
IDEs are Integrated Development Environments, which allow you to edit, execute and inter-
act with code from an integrated environment.
One of the most popular in recent times is VS Code, which is now available via Anaconda.
We hear good things about VS Code β please tell us about your experiences on the forum.
2.8 Exercises
2.8.1 Exercise 1
If Jupyter is still running, quit by using Ctrl-C at the terminal where you started it.
Now launch again, but this time using jupyter notebook --no-browser.
This should start the kernel without launching the browser.
Note also the startup message: It should give you a URL such as
http://localhost:8888 where the notebook is running.
Now
2.8.2 Exercise 2
As an exercise, try
1. Installing Git.
2. Getting a copy of QuantEcon.py using Git.
For example, if youβve installed the command line version, open up a terminal and enter.
(This is just git clone in front of the URL for the repository)
Even better,
1. Sign up to GitHub.
2. Look into βforkingβ GitHub repositories (forking means making your own copy of a
GitHub repository, stored on GitHub).
3. Fork QuantEcon.py.
4. Clone your fork to some local directory, make edits, commit them, and push them back
up to your forked GitHub repo.
5. If you made a valuable improvement, send us a pull request!
An Introductory Example
3.1 Contents
β’ Overview 3.2
β’ Version 1 3.4
β’ Exercises 3.6
β’ Solutions 3.7
Note: These references offer help on installing Python but you should probably stick with the
method on our set up page.
Youβll then have an outstanding scientific computing environment (Anaconda) and be ready
to move on to the rest of our course.
3.2 Overview
In this lecture, we will write and then pick apart small Python programs.
35
36 CHAPTER 3. AN INTRODUCTORY EXAMPLE
The objective is to introduce you to basic Python syntax and data structures.
Deeper concepts will be covered in later lectures.
3.2.1 Prerequisites
Suppose we want to simulate and plot the white noise process π0 , π1 , β¦ , ππ , where each draw
ππ‘ is independent standard normal.
In other words, we want to generate figures that look something like this:
3.4 Version 1
Here are a few lines of code that perform the task we set
import numpy as np
[1]: import matplotlib.pyplot as plt
%matplotlib inline
x = np.random.randn(100)
plt.plot(x)
plt.show()
3.4. VERSION 1 37
After import numpy as np we have access to these attributes via the syntax np..
Hereβs another example
import numpy as np
[2]:
np.sqrt(4)
2.0
[2]:
2.0
[3]:
In fact, you can find and explore the directory for NumPy on your computer easily enough if
you look around.
On this machine, itβs located in
anaconda3/lib/python3.6/site-packages/numpy
Subpackages
Consider the line x = np.random.randn(100).
Here np refers to the package NumPy, while random is a subpackage of NumPy.
You can see the contents here.
Subpackages are just packages that are subdirectories of another package.
2.0
[4]:
2.0
[5]:
3.5. ALTERNATIVE VERSIONS 39
for i in range(ts_length):
e = np.random.randn()
οΏ½_values.append(e)
plt.plot(οΏ½_values)
plt.show()
40 CHAPTER 3. AN INTRODUCTORY EXAMPLE
In brief,
3.5.2 Lists
list
[7]:
The first element of x is an integer, the next is a string and the third is a Boolean value.
When adding a value to a list, we can use the syntax list_name.append(some_value)
x
[8]:
[10, 'foo', False]
[8]:
x.append(2.5)
[9]: x
Here append() is whatβs called a method, which is a function βattached toβ an objectβin
this case, the list x.
Weβll learn all about methods later on, but just to give you some idea,
β’ Python objects such as lists, strings, etc. all have methods that are used to manipulate
the data contained in the object.
β’ String objects have string methods, list objects have list methods, etc.
Now letβs consider the for loop from the program above, which was
for i in range(ts_length):
[16]: e = np.random.randn()
οΏ½_values.append(e)
Python executes the two indented lines ts_length times before moving on.
These two lines are called a code block, since they comprise the βblockβ of code that we
are looping over.
Unlike most other languages, Python knows the extent of the code block only from indenta-
tion.
In our program, indentation decreases after line οΏ½_values.append(e), telling Python that
this line marks the lower limit of the code block.
More on indentation belowβfor now, letβs look at another example of a for loop
animals = ['dog', 'cat', 'bird']
[17]: for animal in animals:
print("The plural of " + animal + " is " + animal + "s")
This example helps to clarify how the for loop works: When we execute a loop of the form
β’ For each element of the sequence, it βbindsβ the name variable_name to that ele-
ment and then executes the code block.
The sequence object can in fact be a very general object, as weβll see soon enough.
42 CHAPTER 3. AN INTRODUCTORY EXAMPLE
In discussing the for loop, we explained that the code blocks being looped over are delimited
by indentation.
In fact, in Python, all code blocks (i.e., those occurring inside loops, if clauses, function defi-
nitions, etc.) are delimited by indentation.
Thus, unlike most other languages, whitespace in Python code affects the output of the pro-
gram.
Once you get used to it, this is a good thing: It
On the other hand, it takes a bit of care to get right, so please remember:
β’ The line before the start of a code block always ends in a colon
β for i in range(10):
β if x > y:
β while x < 100:
β etc., etc.
β’ All lines in a code block must have the same amount of indentation.
β’ The Python standard is 4 spaces, and thatβs what you should use.
Tabs vs Spaces
One small βgotchaβ here is the mixing of tabs and spaces, which often leads to errors.
(Important: Within text files, the internal representation of tabs and spaces is not the same)
You can use your Tab key to insert 4 spaces, but you need to make sure itβs configured to do
so.
If you are using a Jupyter notebook you will have no problems here.
Also, good text editors will allow you to configure the Tab key to insert spaces instead of tabs
β trying searching online.
The for loop is the most common technique for iteration in Python.
But, for the purpose of illustration, letβs modify the program above to use a while loop in-
stead.
ts_length = 100
[18]: οΏ½_values = []
i = 0
while i < ts_length:
e = np.random.randn()
οΏ½_values.append(e)
i = i + 1
plt.plot(οΏ½_values)
3.5. ALTERNATIVE VERSIONS 43
plt.show()
Note that
β’ the code block for the while loop is again delimited only by indentation
β’ the statement i = i + 1 can be replaced by i += 1
Now letβs go back to the for loop, but restructure our program to make the logic clearer.
To this end, we will break our program into two parts:
data = generate_data(100)
plt.plot(data)
plt.show()
44 CHAPTER 3. AN INTRODUCTORY EXAMPLE
Letβs go over this carefully, in case youβre not familiar with functions and how they work.
We have defined a function called generate_data() as follows
This whole function definition is read by the Python interpreter and stored in memory.
When the interpreter gets to the expression generate_data(100), it executes the function
body with n set equal to 100.
The net result is that the name data is bound to the list οΏ½_values returned by the func-
tion.
3.5.7 Conditions
else:
e = np.random.randn()
οΏ½_values.append(e)
return οΏ½_values
Hopefully, the syntax of the if/else clause is self-explanatory, with indentation again delimit-
ing the extent of the code blocks.
Notes
β For example, the statement a = 10 assigns the name a to the value 10.
β The expression a == 10 evaluates to either True or False, depending on the
value of a.
Now, there are several ways that we can simplify the code above.
For example, we can get rid of the conditionals all together by just passing the desired gener-
ator type as a function.
To understand this, consider the following version.
def generate_data(n, generator_type):
[21]: οΏ½_values = []
for i in range(n):
e = generator_type()
οΏ½_values.append(e)
return οΏ½_values
46 CHAPTER 3. AN INTRODUCTORY EXAMPLE
This principle works more generallyβfor example, consider the following piece of code
max(7, 2, 4) # max() is a built-in Python function
[22]:
7
[22]:
m = max
[23]: m(7, 2, 4)
7
[23]:
Here we created another name for the built-in function max(), which could then be used in
identical ways.
In the context of our program, the ability to bind new names to functions means that there is
no problem passing a function as an argument to another functionβas we did above.
3.6. EXERCISES 47
We can also simplify the code for generating the list of random draws considerably by using
something called a list comprehension.
List comprehensions are an elegant Python tool for creating lists.
Consider the following example, where the list comprehension is on the right-hand side of the
second line
animals = ['dog', 'cat', 'bird']
[24]: plurals = [animal + 's' for animal in animals]
plurals
οΏ½_values = []
for i in range(n):
e = generator_type()
οΏ½_values.append(e)
into
3.6 Exercises
3.6.1 Exercise 1
3.6.2 Exercise 2
The binomial random variable π βΌ π΅ππ(π, π) represents the number of successes in π binary
trials, where each trial succeeds with probability π.
48 CHAPTER 3. AN INTRODUCTORY EXAMPLE
Without any import besides from numpy.random import uniform, write a function
binomial_rv such that binomial_rv(n, p) generates one draw of π .
Hint: If π is uniform on (0, 1) and π β (0, 1), then the expression U < p evaluates to True
with probability π.
3.6.3 Exercise 3
β’ If π is a bivariate uniform random variable on the unit square (0, 1)2 , then the proba-
bility that π lies in a subset π΅ of (0, 1)2 is equal to the area of π΅.
β’ If π1 , β¦ , ππ are IID copies of π , then, as π gets large, the fraction that falls in π΅, con-
verges to the probability of landing in π΅.
β’ For a circle, area = pi * radius^2.
3.6.4 Exercise 4
Write a program that prints one realization of the following random device:
3.6.5 Exercise 5
Your next task is to simulate and plot the correlated time series
3.6.6 Exercise 6
To do the next exercise, you will need to know how to produce a plot legend.
The following example should be sufficient to convey the idea
3.6. EXERCISES 49
import numpy as np
[29]: import matplotlib.pyplot as plt
Now, starting with your solution to exercise 5, plot three simulated time series, one for each
of the cases πΌ = 0, πΌ = 0.8 and πΌ = 0.98.
In particular, you should produce (modulo randomness) a figure that looks as follows
50 CHAPTER 3. AN INTRODUCTORY EXAMPLE
(The figure nicely illustrates how time series with the same one-step-ahead conditional volatil-
ities, as these three processes have, can have very different unconditional volatilities.)
Use a for loop to step through the πΌ values.
Important hints:
β’ If you call the plot() function multiple times before calling show(), all of the lines
you produce will end up on the same figure.
β And if you omit the argument 'b-' to the plot function, Matplotlib will automati-
cally select different colors for each line.
3.7 Solutions
3.7.1 Exercise 1
def factorial(n):
[30]: k = 1
for i in range(n):
k = k * (i + 1)
return k
factorial(4)
24
[30]:
3.7. SOLUTIONS 51
3.7.2 Exercise 2
from numpy.random import uniform
[31]:
def binomial_rv(n, p):
count = 0
for i in range(n):
U = uniform()
if U < p:
count = count + 1 # Or count += 1
return count
binomial_rv(10, 0.5)
5
[31]:
3.7.3 Exercise 3
area_estimate = count / n
3.1448
3.7.4 Exercise 4
from numpy.random import uniform
[33]:
payoff = 0
count = 0
for i in range(10):
U = uniform()
count = count + 1 if U < 0.5 else 0
if count == 3:
payoff = 1
print(payoff)
1
52 CHAPTER 3. AN INTRODUCTORY EXAMPLE
3.7.5 Exercise 5
The next line embeds all subsequent figures in the browser itself
Ξ± = 0.9
[34]: ts_length = 200
current_x = 0
x_values = []
for i in range(ts_length + 1):
x_values.append(current_x)
current_x = Ξ± * current_x + np.random.randn()
plt.plot(x_values)
plt.show()
3.7.6 Exercise 6
for Ξ± in Ξ±s:
x_values = []
current_x = 0
for i in range(ts_length):
x_values.append(current_x)
current_x = Ξ± * current_x + np.random.randn()
plt.plot(x_values, label=f'Ξ± = {Ξ±}')
plt.legend()
plt.show()
3.7. SOLUTIONS 53
54 CHAPTER 3. AN INTRODUCTORY EXAMPLE
Chapter 4
Python Essentials
4.1 Contents
In this lecture, weβll cover features of the language that are essential to reading and writing
Python code.
Weβve already met several built-in Python data types, such as strings, integers, floats and
lists.
Letβs learn a bit more about them.
One simple data type is Boolean values, which can be either True or False
x = True
[1]: x
True
[1]:
In the next line of code, the interpreter evaluates the expression on the right of = and binds y
to this value
55
56 CHAPTER 4. PYTHON ESSENTIALS
y = 100 < 10
[2]: y
False
[2]:
type(y)
[3]:
bool
[3]:
3
[7]:
The two most common data types used to represent numbers are integers and floats
a, b = 1, 2
[8]: c, d = 2.5, 10.0
type(a)
int
[8]:
type(c)
[9]:
float
[9]:
Computers distinguish between the two because, while floats are more informative, arithmetic
operations on integers are faster and more accurate.
As long as youβre using Python 3.x, division of integers yields floats
1 / 2
[10]:
0.5
[10]:
But be careful! If youβre still using Python 2.x, division of two integers returns only the inte-
ger part.
For integer division in Python 3.x use this syntax:
1 // 2
[11]:
0
[11]:
x = complex(1, 2)
[12]: y = complex(2, 1)
x * y
5j
[12]:
4.2.2 Containers
Python has several basic types for storing collections of (possibly heterogeneous) data.
Weβve already discussed lists.
A related data type is tuples, which are βimmutableβ lists
x = ('a', 'b') # Parentheses instead of the square brackets
[13]: x = 'a', 'b' # Or no brackets --- the meaning is identical
x
('a', 'b')
[13]:
type(x)
[14]:
tuple
[14]:
In Python, an object is called immutable if, once created, the object cannot be changed.
Conversely, an object is mutable if it can still be altered after creation.
Python lists are mutable
x = [1, 2]
[15]: x[0] = 10
x
[10, 2]
[15]:
---------------------------------------------------------------------------
<ipython-input-16-d1b2647f6c81> in <module>
1 x = (1, 2)
----> 2 x[0] = 10
Weβll say more about the role of mutable and immutable data a bit later.
Tuples (and lists) can be βunpackedβ as follows
integers = (10, 20, 30)
[17]: x, y, z = integers
x
10
[17]:
58 CHAPTER 4. PYTHON ESSENTIALS
y
[18]:
20
[18]:
[4, 6, 8]
[19]:
a[1:3]
[20]:
[4, 6]
[20]:
'bar'
[22]:
dict
[23]:
d['age']
[24]:
33
[24]:
set
[25]:
4.3. INPUT AND OUTPUT 59
s2 = {'b', 'c'}
[26]: s1.issubset(s2)
False
[26]:
s1.intersection(s2)
[27]:
{'b'}
[27]:
{'bar', 'foo'}
[28]:
Letβs briefly review reading and writing to text files, starting with writing
f = open('newfile.txt', 'w') # Open 'newfile.txt' for writing
[29]: f.write('Testing\n') # Here '\n' means new line
f.write('Testing again')
f.close()
Here
β’ The built-in function open() creates a file object for writing to.
β’ Both write() and close() are methods of file objects.
'Testing\nTesting again'
[31]:
print(out)
[32]:
Testing
Testing again
4.3.1 Paths
Note that if newfile.txt is not in the present working directory then this call to open()
fails.
60 CHAPTER 4. PYTHON ESSENTIALS
In this case, you can shift the file to the pwd or specify the full path to the file
f = open('insert_full_path_to_file/newfile.txt', 'r')
4.4 Iterating
One of the most important tasks in computing is stepping through a sequence of data and
performing a given action.
One of Pythonβs strengths is its simple, flexible interface to this kind of iteration via the for
loop.
Many Python objects are βiterableβ, in the sense that they can be looped over.
To give an example, letβs write the file us_cities.txt, which lists US cities and their popula-
tion, to the present working directory.
%%file us_cities.txt
[33]: new york: 8244910
los angeles: 3819702
chicago: 2707120
houston: 2145146
philadelphia: 1536471
phoenix: 1469471
san antonio: 1359758
san diego: 1326179
dallas: 1223229
Overwriting us_cities.txt
Suppose that we want to make the information more readable, by capitalizing names and
adding commas to mark thousands.
The program below reads the data in and makes the conversion:
data_file = open('us_cities.txt', 'r')
[34]: for line in data_file:
city, population = line.split(':') # Tuple unpacking
city = city.title() # Capitalize city names
population = f'{int(population):,}' # Add commas to numbers
print(city.ljust(15) + population)
data_file.close()
Here format() is a string method used for inserting variables into strings.
The reformatting of each line is the result of three different string methods, the details of
which can be left till later.
4.4. ITERATING 61
The interesting part of this program for us is line 2, which shows that
1. The file object f is iterable, in the sense that it can be placed to the right of in within
a for loop.
2. Iteration steps through each line in the file.
One thing you might have noticed is that Python tends to favor looping without explicit in-
dexing.
For example,
x_values = [1, 2, 3] # Some iterable x
[35]: for x in x_values:
print(x * x)
1
4
9
is preferred to
for i in range(len(x_values)):
[36]: print(x_values[i] * x_values[i])
1
4
9
When you compare these two alternatives, you can see why the first one is preferred.
Python provides some facilities to simplify looping without indices.
One is zip(), which is used for stepping through pairs from two sequences.
For example, try running the following code
countries = ('Japan', 'Korea', 'China')
[37]: cities = ('Tokyo', 'Seoul', 'Beijing')
for country, city in zip(countries, cities):
print(f'The capital of {country} is {city}')
The zip() function is also useful for creating dictionaries β for example
names = ['Tom', 'John']
[38]: marks = ['E', 'F']
dict(zip(names, marks))
62 CHAPTER 4. PYTHON ESSENTIALS
If we actually need the index from a list, one option is to use enumerate().
To understand what enumerate() does, consider the following example
letter_list = ['a', 'b', 'c']
[39]: for index, letter in enumerate(letter_list):
print(f"letter_list[{index}] = '{letter}'")
letter_list[0] = 'a'
letter_list[1] = 'b'
letter_list[2] = 'c'
4.5.1 Comparisons
Many different kinds of expressions evaluate to one of the Boolean values (i.e., True or
False).
A common type is comparisons, such as
x, y = 1, 2
[41]: x < y
True
[41]:
x > y
[42]:
False
[42]:
False
[45]:
Note that when testing conditions, we can use any valid Python expression
x = 'yes' if 42 else 'no'
[47]: x
'yes'
[47]:
x = 'yes' if [] else 'no'
[48]: x
'no'
[48]:
β’ Expressions that evaluate to zero, empty sequences or containers (strings, lists, etc.)
and None are all equivalent to False.
Remember
Letβs talk a bit more about functions, which are all important for good programming style.
Python has a number of built-in functions that are available without import.
We have already met some
max(19, 20)
[54]:
20
[54]:
range(4) # in python3 this returns a range iterator object
[55]:
range(0, 4)
[55]:
list(range(4)) # will evaluate the range iterator and create a list
[56]:
[0, 1, 2, 3]
[56]:
str(22)
[57]:
'22'
[57]:
type(22)
[58]:
int
[58]:
False
[59]:
any(bools) # False if all are False and True otherwise
[]:
True
[]:
User-defined functions are important for improving the clarity of your code by
Functions without a return statement automatically return the special Python object None.
4.6.3 Docstrings
Python has a system for adding comments to functions, modules, etc. called docstrings.
The nice thing about docstrings is that they are available at run-time.
Try running this
def f(x):
[62]: """
This function squares its argument
"""
return x**2
Type: function
String Form:<function f at 0x2223320>
File: /home/john/temp/temp.py
Definition: f(x)
Docstring: This function squares its argument
f??
[64]:
Type: function
String Form:<function f at 0x2223320>
File: /home/john/temp/temp.py
Definition: f(x)
Source:
def f(x):
"""
This function squares its argument
"""
return x**2
66 CHAPTER 4. PYTHON ESSENTIALS
With one question mark we bring up the docstring, and with two we get the source code as
well.
and
f = lambda x: x**3
[66]:
(4.0, 4.440892098500626e-14)
[67]:
Here the function created by lambda is said to be anonymous because it was never given a
name.
If you did the exercises in the previous lecture, you would have come across the statement
In this call to Matplotlibβs plot function, notice that the last argument is passed in
name=argument syntax.
This is called a keyword argument, with label being the keyword.
Non-keyword arguments are called positional arguments, since their meaning is determined by
order
Keyword arguments are particularly useful when a function has a lot of arguments, in which
case itβs hard to remember the right order.
4.7. CODING STYLE AND PEP8 67
The keyword argument values we supplied in the definition of f become the default values
f(2)
[69]:
3
[69]:
To learn more about the Python programming philosophy type import this at the
prompt.
Among other things, Python strongly favors consistency in programming style.
Weβve all heard the saying about consistency and little minds.
In programming, as in mathematics, the opposite is true
β’ A mathematical paper where the symbols βͺ and β© were reversed would be very hard to
read, even if the author told you so on the first page.
4.8 Exercises
4.8.1 Exercise 1
Part 1: Given two numeric lists or tuples x_vals and y_vals of equal length, compute their
inner product using zip().
Part 2: In one line, count the number of even numbers in 0,β¦,99.
Part 3: Given pairs = ((2, 5), (4, 2), (9, 8), (12, 10)), count the number of
pairs (a, b) such that both a and b are even.
68 CHAPTER 4. PYTHON ESSENTIALS
4.8.2 Exercise 2
π
π(π₯) = π0 + π1 π₯ + π2 π₯2 + β― ππ π₯π = β ππ π₯π (1)
π=0
Write a function p such that p(x, coeff) that computes the value in Eq. (1) given a point
x and a list of coefficients coeff.
Try to use enumerate() in your loop.
4.8.3 Exercise 3
Write a function that takes a string as an argument and returns the number of capital letters
in the string.
Hint: 'foo'.upper() returns 'FOO'.
4.8.4 Exercise 4
Write a function that takes two sequences seq_a and seq_b as arguments and returns True
if every element in seq_a is also an element of seq_b, else False.
4.8.5 Exercise 5
When we cover the numerical libraries, we will see they include many alternatives for interpo-
lation and function approximation.
Nevertheless, letβs write our own function approximation routine as an exercise.
In particular, without using any imports, write a function linapprox that takes as argu-
ments
and returns the piecewise linear interpolation of f at x, based on n evenly spaced grid points
a = point[0] < point[1] < ... < point[n-1] = b.
Aim for clarity, not efficiency.
4.9. SOLUTIONS 69
4.9 Solutions
4.9.1 Exercise 1
Part 1 Solution:
Hereβs one possible solution
x_vals = [1, 2, 3]
[71]: y_vals = [1, 1, 1]
sum([x * y for x, y in zip(x_vals, y_vals)])
6
[71]:
Part 2 Solution:
One solution is
sum([x % 2 == 0 for x in range(100)])
[73]:
50
[73]:
Some less natural alternatives that nonetheless help to illustrate the flexibility of list compre-
hensions are
len([x for x in range(100) if x % 2 == 0])
[75]:
50
[75]:
and
sum([1 for x in range(100) if x % 2 == 0])
[76]:
50
[76]:
Part 3 Solution
Hereβs one possibility
pairs = ((2, 5), (4, 2), (9, 8), (12, 10))
[77]: sum([x % 2 == 0 and y % 2 == 0 for x, y in pairs])
2
[77]:
4.9.2 Exercise 2
def p(x, coeff):
[78]: return sum(a * x**i for i, a in enumerate(coeff))
6
[79]:
4.9.3 Exercise 3
3
[80]:
3
[81]:
4.9.4 Exercise 4
Hereβs a solution:
def f(seq_a, seq_b):
[82]: is_subset = True
for a in seq_a:
if a not in seq_b:
is_subset = False
return is_subset
# == test == #
True
False
Of course, if we use the sets data type then the solution is easier
def f(seq_a, seq_b):
[83]: return set(seq_a).issubset(set(seq_b))
4.9.5 Exercise 5
def linapprox(f, a, b, n, x):
[84]: """
Evaluates the piecewise linear interpolant of f at x on the interval
[a, b], with n evenly spaced grid points.
Parameters
==========
f : function
The function to approximate
4.9. SOLUTIONS 71
n : integer
Number of grid points
Returns
=======
A float. The interpolant evaluated at x
"""
length_of_interval = b - a
num_subintervals = n - 1
step = length_of_interval / num_subintervals
# === x must lie between the gridpoints (point - step) and point === #
u, v = point - step, point
5.1 Contents
β’ Overview 5.2
β’ Objects 5.3
β’ Summary 5.4
5.2 Overview
Python is a pragmatic language that blends object-oriented and procedural styles, rather than
taking a purist approach.
However, at a foundational level, Python is object-oriented.
73
74 CHAPTER 5. OOP I: INTRODUCTION TO OBJECT ORIENTED PROGRAMMING
5.3 Objects
In Python, an object is a collection of data and instructions held in computer memory that
consists of
1. a type
2. a unique identity
3. data (i.e., content)
4. methods
5.3.1 Type
Python provides for different types of objects, to accommodate different categories of data.
For example
s = 'This is a string'
[1]: type(s)
str
[1]:
x = 42 # Now let's create an integer
[2]: type(x)
int
[2]:
---------------------------------------------------------------------------
<ipython-input-5-263a89d2d982> in <module>
----> 1 '300' + 400
5.3. OBJECTS 75
Here we are mixing types, and itβs unclear to Python whether the user wants to
To avoid the error, you need to clarify by changing the relevant type.
For example,
int('300') + 400 # To add as numbers, change the string to an integer
[6]:
700
[6]:
5.3.2 Identity
In Python, each object has a unique identifier, which helps Python (and us) keep track of the
object.
The identity of an object can be obtained via the id() function
y = 2.5
[7]: z = 2.5
id(y)
140473280865624
[7]:
id(z)
[8]:
140473280865648
[8]:
In this example, y and z happen to have the same value (i.e., 2.5), but they are not the
same object.
The identity of an object is in fact just the address of the object in memory.
If we set x = 42 then we create an object of type int that contains the data 42.
In fact, it contains more, as the following example shows
x = 42
[9]: x
42
[9]:
x.imag
[10]:
76 CHAPTER 5. OOP I: INTRODUCTION TO OBJECT ORIENTED PROGRAMMING
0
[10]:
x.__class__
[11]:
int
[11]:
When Python creates this integer object, it stores with it various auxiliary information, such
as the imaginary part, and the type.
Any name following a dot is called an attribute of the object to the left of the dot.
We see from this example that objects have attributes that contain auxiliary information.
They also have attributes that act like functions, called methods.
These attributes are important, so letβs discuss them in-depth.
5.3.4 Methods
True
[12]:
callable(x.__doc__)
[13]:
False
[13]:
Methods typically act on the data contained in the object they belong to, or combine that
data with other data
x = ['a', 'b']
[14]: x.append('c')
s = 'This is a string'
s.upper()
'THIS IS A STRING'
[14]:
s.lower()
[15]:
'this is a string'
[15]:
s.replace('This', 'That')
[16]:
'That is a string'
[16]:
['aa', 'b']
[17]:
5.4. SUMMARY 77
It doesnβt look like there are any methods used here, but in fact the square bracket assign-
ment notation is just a convenient interface to a method call.
What actually happens is that Python calls the __setitem__ method, as follows
x = ['a', 'b']
[18]: x.__setitem__(0, 'aa') # Equivalent to x[0] = 'aa'
x
['aa', 'b']
[18]:
(If you wanted to you could modify the __setitem__ method, so that square bracket as-
signment does something totally different)
5.4 Summary
<function __main__.f(x)>
[19]:
type(f)
[20]:
function
[20]:
id(f)
[21]:
140473061744848
[21]:
f.__name__
[22]:
'f'
[22]:
We can see that f has type, identity, attributes and so onβjust like any other object.
It also has methods.
One example is the __call__ method, which just evaluates the function
f.__call__(3)
[23]:
9
[23]:
140473369783576
[24]:
This uniform treatment of data in Python (everything is an object) helps keep the language
simple and consistent.
Part II
79
Chapter 6
NumPy
6.1 Contents
β’ Overview 6.2
β’ Exercises 6.7
β’ Solutions 6.8
βLetβs be clear: the work of science has nothing whatever to do with consensus.
Consensus is the business of politics. Science, on the contrary, requires only one
investigator who happens to be right, which means that he or she has results that
are verifiable by reference to the real world. In science consensus is irrelevant.
What is relevant is reproducible results.β β Michael Crichton
6.2 Overview
In this lecture, we introduce NumPy arrays and the fundamental array processing operations
provided by NumPy.
6.2.1 References
81
82 CHAPTER 6. NUMPY
β’ Loops in Python over Python data types like lists carry significant overhead.
β’ C and Fortran code contains a lot of type information that can be used for optimiza-
tion.
β’ Various optimizations can be carried out during compilation when the compiler sees the
instructions as a whole.
However, for a task like the one described above, thereβs no need to switch back to C or For-
tran.
Instead, we can use NumPy, where the instructions look like this:
import numpy as np
[1]:
x = np.random.uniform(0, 1, size=1000000)
x.mean()
0.49990572835210734
[1]:
The operations of creating the array and computing its mean are both passed out to carefully
optimized machine code compiled from C.
More generally, NumPy sends operations in batches to optimized C and Fortran code.
This is similar in spirit to Matlab, which provides an interface to fast Fortran routines.
In a later lecture, weβll discuss code that isnβt easy to vectorize and how such routines can
also be optimized.
The most important thing that NumPy defines is an array data type formally called a
numpy.ndarray.
NumPy arrays power a large proportion of the scientific Python ecosystem.
6.4. NUMPY ARRAYS 83
NumPy arrays are somewhat like native Python lists, except that
There are also dtypes to represent complex numbers, unsigned integers, etc.
On modern machines, the default dtype for arrays is float64
a = np.zeros(3)
[4]: type(a[0])
numpy.float64
[4]:
numpy.int64
[5]:
Here z is a flat array with no dimension β neither row nor column vector.
The dimension is recorded in the shape attribute, which is a tuple
z.shape
[7]:
(10,)
[7]:
Here the shape tuple has only one element, which is the length of the array (tuples with one
element end with a comma).
To give it dimension, we can change the shape attribute
z.shape = (10, 1)
[8]: z
84 CHAPTER 6. NUMPY
array([[0.],
[8]: [0.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.]])
z = np.zeros(4)
[9]: z.shape = (2, 2)
z
array([[0., 0.],
[9]: [0., 0.]])
In the last case, to make the 2 by 2 array, we could also pass a tuple to the zeros() func-
tion, as in z = np.zeros((2, 2)).
array([[1., 0.],
[12]: [0., 1.]])
In addition, NumPy arrays can be created from Python lists, tuples, etc. using np.array
z = np.array([10, 20]) # ndarray from Python list
[13]: z
array([10, 20])
[13]:
type(z)
[14]:
numpy.ndarray
[14]:
z = np.array((10, 20), dtype=float) # Here 'float' is equivalent to 'np.float64'
[15]: z
6.4. NUMPY ARRAYS 85
array([10., 20.])
[15]:
z = np.array([[1, 2], [3, 4]]) # 2D array from a list of lists
[16]: z
array([[1, 2],
[16]: [3, 4]])
See also np.asarray, which performs a similar function, but does not make a distinct copy
of data already in a NumPy array.
na = np.linspace(10, 20, 2)
[17]: na is np.asarray(na) # Does not copy NumPy arrays
True
[17]:
na is np.array(na) # Does make a new copy --- perhaps unnecessarily
[18]:
False
[18]:
To read in the array data from a text file containing numeric data use np.loadtxt or
np.genfromtxtβsee the documentation for details.
array([[1, 2],
[23]: [3, 4]])
z[0, 0]
[24]:
1
[24]:
z[0, 1]
[25]:
2
[25]:
And so on.
86 CHAPTER 6. NUMPY
Note that indices are still zero-based, to maintain compatibility with Python sequences.
Columns and rows can be extracted as follows
z[0, :]
[26]:
array([1, 2])
[26]:
z[:, 1]
[27]:
array([2, 4])
[27]:
array([2. , 3. , 3.5])
[29]:
array([2. , 3. , 3.5])
[33]:
z[:] = 42
[34]: z
array([4, 3, 2, 1])
[35]:
6.4. NUMPY ARRAYS 87
array([1, 2, 3, 4])
[36]:
a.sum() # Sum
[37]:
10
[37]:
a.mean() # Mean
[38]:
2.5
[38]:
a.max() # Max
[39]:
4
[39]:
a.argmax() # Returns the index of the maximal element
[40]:
3
[40]:
a.cumsum() # Cumulative sum of the elements of a
[41]:
array([ 1, 3, 6, 10])
[41]:
a.cumprod() # Cumulative product of the elements of a
[42]:
array([ 1, 2, 6, 24])
[42]:
a.var() # Variance
[43]:
1.25
[43]:
a.std() # Standard deviation
[44]:
1.118033988749895
[44]:
a.shape = (2, 2)
[45]: a.T # Equivalent to a.transpose()
array([[1, 3],
[45]: [2, 4]])
Many of the methods discussed above have equivalent functions in the NumPy namespace
a = np.array((4, 3, 2, 1))
[48]:
np.sum(a)
[49]:
10
[49]:
88 CHAPTER 6. NUMPY
np.mean(a)
[50]:
2.5
[50]:
array([[2., 2.],
[55]: [2., 2.]])
A + 10
[56]:
array([[11., 11.],
[56]: [11., 11.]])
A * B
[57]:
array([[1., 1.],
[57]: [1., 1.]])
With Anacondaβs scientific Python package based around Python 3.5 and above, one can use
the @ symbol for matrix multiplication, as follows:
6.5. OPERATIONS ON ARRAYS 89
A = np.ones((2, 2))
[58]: B = np.ones((2, 2))
A @ B
array([[2., 2.],
[58]: [2., 2.]])
(For older versions of Python and NumPy you need to use the np.dot function)
We can also use @ to take the inner product of two flat arrays
A = np.array((1, 2))
[59]: B = np.array((10, 20))
A @ B
50
[59]:
array([[1, 2],
[]: [3, 4]])
A @ (0, 1)
[61]:
array([2, 4])
[61]:
array([42, 44])
[62]:
a[-1] = 0 # Change last element to 0
[63]: a
array([42, 0])
[63]:
Mutability leads to the following behavior (which can be shocking to MATLAB program-
mersβ¦)
a = np.random.randn(3)
[64]: a
NumPy provides versions of the standard functions log, exp, sin, etc. that act element-
wise on arrays
z = np.array([1, 2, 3])
[70]: np.sin(z)
y[i] = np.sin(z[i])
Because they act element-wise on arrays, these functions are called vectorized functions.
In NumPy-speak, they are also called ufuncs, which stands for βuniversal functionsβ.
As we saw above, the usual arithmetic operations (+, *, etc.) also work element-wise, and
combining these with the ufuncs gives a very large set of fast element-wise functions.
z
[72]:
array([1, 2, 3])
[72]:
(1 / np.sqrt(2 * np.pi)) * np.exp(- 0.5 * z**2)
[73]:
array([0.24197072, 0.05399097, 0.00443185])
[73]:
array([1, 0, 0, 1])
[77]:
However, this approach doesnβt always obtain the same speed as a more carefully crafted vec-
torized function.
6.6.2 Comparisons
array([False, True])
[79]:
92 CHAPTER 6. NUMPY
z != y
[80]:
array([ True, False])
[80]:
6.6.3 Sub-packages
NumPy provides some additional functionality related to scientific programming through its
sub-packages.
Weβve already seen how we can generate random variables using np.random
z = np.random.randn(10000) # Generate standard normals
[86]: y = np.random.binomial(10, 0.5, size=1000) # 1,000 draws from Bin(10, 0.5)
y.mean()
5.013
[86]:
-2.0000000000000004
[87]:
np.linalg.inv(A) # Compute the inverse
[88]:
array([[-2. , 1. ],
[88]: [ 1.5, -0.5]])
Much of this functionality is also available in SciPy, a collection of modules that are built on
top of NumPy.
6.7. EXERCISES 93
6.7 Exercises
6.7.1 Exercise 1
π
π(π₯) = π0 + π1 π₯ + π2 π₯2 + β― ππ π₯π = β ππ π₯π (1)
π=0
Earlier, you wrote a simple function p(x, coeff) to evaluate Eq. (1) without considering
efficiency.
Now write a new function that does the same job, but uses NumPy arrays and array opera-
tions for its computations, rather than any form of Python loop.
(Such functionality is already implemented as np.poly1d, but for the sake of the exercise
donβt use this class)
6.7.2 Exercise 2
β’ Divide the unit interval [0, 1] into π subintervals πΌ0 , πΌ1 , β¦ , πΌπβ1 such that the length of
πΌπ is ππ .
β’ Draw a uniform random variable π on [0, 1] and return the π such that π β πΌπ .
If you canβt see how this works, try thinking through the flow for a simple example, such as q
= [0.25, 0.75] It helps to sketch the intervals on paper.
94 CHAPTER 6. NUMPY
If you can, write the method so that draw(k) returns k draws from q.
6.7.3 Exercise 3
6.8 Solutions
import matplotlib.pyplot as plt
[90]: %matplotlib inline
6.8.1 Exercise 1
Letβs test it
x = 2
[92]: coef = np.linspace(2, 4, 3)
print(coef)
print(p(x, coef))
# For comparison
q = np.poly1d(np.flip(coef))
print(q(x))
[2. 3. 4.]
24.0
24.0
6.8.2 Exercise 2
class DiscreteRV:
"""
Generates an array of draws from a discrete random variable with vector of
probabilities given by q.
"""
The logic is not obvious, but if you take your time and read it slowly, you will understand.
There is a problem here, however.
Suppose that q is altered after an instance of discreteRV is created, for example by
q = (0.1, 0.9)
[94]: d = DiscreteRV(q)
d.q = (0.5, 0.5)
The problem is that Q does not change accordingly, and Q is the data used in the draw
method.
To deal with this, one option is to compute Q every time the draw method is called.
But this is inefficient relative to computing Q once-off.
A better option is to use descriptors.
A solution from the quantecon library using descriptors that behaves as we desire can be
found here.
6.8.3 Exercise 3
"""
class ECDF:
"""
One-dimensional empirical distribution function given a vector of
observations.
Parameters
----------
observations : array_like
An array of observations
96 CHAPTER 6. NUMPY
Attributes
----------
observations : array_like
An array of observations
"""
Parameters
----------
x : scalar(float)
The x at which the ecdf is evaluated
Returns
-------
scalar(float)
Fraction of the sample less than x
"""
return np.mean(self.observations <= x)
Parameters
----------
a : scalar(float), optional(default=None)
Lower endpoint of the plot interval
b : scalar(float), optional(default=None)
Upper endpoint of the plot interval
"""
Matplotlib
7.1 Contents
β’ Overview 7.2
β’ Exercises 7.6
β’ Solutions 7.7
7.2 Overview
Weβve already generated quite a few figures in these lectures using Matplotlib.
Matplotlib is an outstanding graphics library, designed for scientific computing, with
99
100 CHAPTER 7. MATPLOTLIB
Hereβs the kind of easy example you might find in introductory treatments
import matplotlib.pyplot as plt
[1]: %matplotlib inline
import numpy as np
This is simple and convenient, but also somewhat limited and un-Pythonic.
For example, in the function calls, a lot of objects get created and passed around without
making themselves known to the programmer.
Python programmers tend to prefer a more explicit style of programming (run import this
in a code block and look at the second line).
This leads us to the alternative, object-oriented Matplotlib API.
Hereβs the code corresponding to the preceding figure using the object-oriented API
fig, ax = plt.subplots()
[2]: ax.plot(x, y, 'b-', linewidth=2)
plt.show()
7.3. THE APIS 101
7.3.3 Tweaks
Weβve also used alpha to make the line slightly transparentβwhich makes it look smoother.
The location of the legend can be changed by replacing ax.legend() with
ax.legend(loc='upper center').
fig, ax = plt.subplots()
[4]: ax.plot(x, y, 'r-', linewidth=2, label='sine function', alpha=0.6)
ax.legend(loc='upper center')
plt.show()
fig, ax = plt.subplots()
[5]: ax.plot(x, y, 'r-', linewidth=2, label='$y=\sin(x)$', alpha=0.6)
ax.legend(loc='upper center')
plt.show()
Matplotlib has a huge array of functions and features, which you can discover over time as
you have need for them.
We mention just a few.
fig, ax = plt.subplots()
x = np.linspace(-4, 4, 150)
for i in range(3):
m, s = uniform(-1, 1), uniform(1, 2)
y = norm.pdf(x, loc=m, scale=s)
current_label = f'$\mu = {m:.2}$'
ax.plot(x, y, linewidth=2, alpha=0.6, label=current_label)
ax.legend()
plt.show()
7.4.3 3D Plots
Perhaps you will find a set of customizations that you regularly use.
Suppose we usually prefer our axes to go through the origin, and to have a grid.
7.5. FURTHER READING 107
Hereβs a nice example from Matthew Doty of how the object-oriented API can be used to
build a custom subplots function that implements these changes.
Read carefully through the code and see if you can follow whatβs going on
def subplots():
[10]: "Custom subplots with axes through the origin"
fig, ax = plt.subplots()
ax.grid()
return fig, ax
1. calls the standard plt.subplots function internally to generate the fig, ax pair,
2. makes the desired customizations to ax, and
3. passes the fig, ax pair back to the calling code.
7.6 Exercises
7.6.1 Exercise 1
7.7 Solutions
7.7.1 Exercise 1
for ΞΈ in ΞΈ_vals:
ax.plot(x, np.cos(np.pi * ΞΈ * x) * np.exp(- x))
plt.show()
7.7. SOLUTIONS 109
110 CHAPTER 7. MATPLOTLIB
Chapter 8
SciPy
8.1 Contents
β’ Statistics 8.3
β’ Optimization 8.5
β’ Integration 8.6
β’ Exercises 8.8
β’ Solutions 8.9
SciPy builds on top of NumPy to provide common tools for scientific programming such as
β’ linear algebra
β’ numerical integration
β’ interpolation
β’ optimization
β’ distributions and random number generation
β’ signal processing
β’ etc., etc
111
112 CHAPTER 8. SCIPY
SciPy is a package that contains various tools that are built on top of NumPy, using its array
data type and related functionality.
In fact, when we import SciPy we also get NumPy, as can be seen from the SciPy initializa-
tion file
# Import numpy symbols to scipy namespace
[1]: import numpy as _num
linalg = None
from numpy import *
from numpy.random import rand, randn
from numpy.fft import fft, ifft
from numpy.lib.scimath import *
__all__ = []
__all__ += _num.__all__
__all__ += ['randn', 'rand', 'fft', 'ifft']
del _num
# Remove the linalg imported from numpy so that the scipy.linalg package
# can be imported.
del linalg
__all__.remove('linalg')
However, itβs more common and better practice to use NumPy functionality explicitly
import numpy as np
[2]:
a = np.identity(3)
8.3 Statistics
π₯(πβ1) (1 β π₯)(πβ1)
π(π₯; π, π) = 1
(0 β€ π₯ β€ 1) (1)
β«0 π’(πβ1) (1 β π’)(πβ1) ππ’
Sometimes we need access to the density itself, or the cdf, the quantiles, etc.
For this, we can use scipy.stats, which provides all of this functionality as well as random
number generation in a single consistent interface.
Hereβs an example of usage
from scipy.stats import beta
[5]: import matplotlib.pyplot as plt
%matplotlib inline
In this code, we created a so-called rv_frozen object, via the call q = beta(5, 5).
The βfrozenβ part of the notation implies that q represents a particular distribution with a
particular set of parameters.
Once weβve done so, we can then generate random numbers, evaluate the density, etc., all
from this fixed distribution
q.cdf(0.4) # Cumulative distribution function
[6]:
114 CHAPTER 8. SCIPY
0.26656768000000003
[6]:
q.pdf(0.4) # Density function
[7]:
2.0901888000000013
[7]:
q.ppf(0.8) # Quantile (inverse cdf) function
[8]:
0.6339134834642708
[8]:
q.mean()
[9]:
0.5
[9]:
identifier = scipy.stats.distribution_name(shape_parameters)
identifier = scipy.stats.distribution_name(shape_parameters,
loc=c, scale=d)
fig, ax = plt.subplots()
ax.hist(obs, bins=40, density=True)
ax.plot(grid, beta.pdf(grid, 5, 5), 'k-', linewidth=2)
plt.show()
8.4. ROOTS AND FIXED POINTS 115
(2.00290964225675, -0.003101209098539087)
[11]:
8.4.1 Bisection
And so on.
This is bisection.
Hereβs a fairly simplistic implementation of the algorithm in Python.
It works for all sufficiently well behaved increasing continuous functions with π(π) < 0 < π(π)
def bisect(f, a, b, tol=10e-5):
[13]: """
Implements the bisection root finding algorithm, assuming that f is a
real-valued function on [a, b] satisfying f(a) < 0 < f(b).
"""
8.4. ROOTS AND FIXED POINTS 117
lower, upper = a, b
In fact, SciPy provides its own bisection function, which we now test using the function π de-
fined in Eq. (2)
from scipy.optimize import bisect
[14]:
bisect(f, 0, 1)
0.4082935042806639
[14]:
β’ When the function is well-behaved, the Newton-Raphson method is faster than bisec-
tion.
β’ When the function is less well-behaved, the Newton-Raphson might fail.
Letβs investigate this using the same function π, first looking at potential instability
from scipy.optimize import newton
[15]:
newton(f, 0.2) # Start the search at initial condition x = 0.2
0.40829350427935673
[15]:
newton(f, 0.7) # Start the search at x = 0.7 instead
[16]:
0.7001700000000279
[16]:
98.3 Β΅s Β± 359 ns per loop (mean Β± std. dev. of 7 runs, 10000 loops each)
198 Β΅s Β± 20.4 Β΅s per loop (mean Β± std. dev. of 7 runs, 10000 loops each)
118 CHAPTER 8. SCIPY
So far we have seen that the Newton-Raphson method is fast but not robust.
This bisection algorithm is robust but relatively slow.
This illustrates a general principle
β’ If you have specific knowledge about your function, you might be able to exploit it to
generate efficiency.
β’ If not, then the algorithm choice involves a trade-off between the speed of convergence
and robustness.
In practice, most default algorithms for root-finding, optimization and fixed points use hybrid
methods.
These methods typically combine a fast method with a robust method in the following man-
ner:
In scipy.optimize, the function brentq is such a hybrid method and a good default
brentq(f, 0, 1)
[19]:
0.40829350427936706
[19]:
%timeit brentq(f, 0, 1)
[20]:
24.9 Β΅s Β± 71.8 ns per loop (mean Β± std. dev. of 7 runs, 10000 loops each)
Here the correct solution is found and the speed is almost the same as newton.
array(1.)
[21]:
If you donβt get good results, you can always switch back to the brentq root finder, since
the fixed point of a function π is the root of π(π₯) βΆ= π₯ β π(π₯).
8.5. OPTIMIZATION 119
8.5 Optimization
0.0
[22]:
8.6 Integration
0.33333333333333337
[23]:
In fact, quad is an interface to a very standard numerical integration routine in the Fortran
library QUADPACK.
It uses Clenshaw-Curtis quadrature, based on expansion in terms of Chebychev polynomials.
There are other options for univariate integrationβa useful one is fixed_quad, which is fast
and hence works well inside for loops.
120 CHAPTER 8. SCIPY
We saw that NumPy provides a module for linear algebra called linalg.
SciPy also provides a module for linear algebra with the same name.
The latter is not an exact superset of the former, but overall it has more functionality.
We leave you to investigate the set of available routines.
8.8 Exercises
8.8.1 Exercise 1
8.9 Solutions
8.9.1 Exercise 1
0.408294677734375
[26]:
122 CHAPTER 8. SCIPY
Chapter 9
Numba
9.1 Contents
β’ Overview 9.2
β’ Vectorization 9.4
β’ Numba 9.5
In addition to whatβs in Anaconda, this lecture will need the following libraries:
!pip install --upgrade quantecon
[1]:
9.2 Overview
In our lecture on NumPy, we learned one method to improve speed and efficiency in numeri-
cal work.
That method, called vectorization, involved sending array processing operations in batch to
efficient low-level code.
This clever idea dates back to Matlab, which uses it extensively.
Unfortunately, vectorization is limited and has several weaknesses.
One weakness is that it is highly memory-intensive.
Another problem is that only some algorithms can be vectorized.
In the last few years, a new Python library called Numba has appeared that solves many of
these problems.
It does so through something called just in time (JIT) compilation.
JIT compilation is effective in many numerical settings and can generate extremely fast, effi-
cient code.
It can also do other tricks such as facilitate multithreading (a form of parallelization well
suited to numerical work).
123
124 CHAPTER 9. NUMBA
To understand what Numba does and why, we need some background knowledge.
Letβs start by thinking about higher-level languages, such as Python.
These languages are optimized for humans.
This means that the programmer can leave many details to the runtime environment
The upside is that, compared to low-level languages, Python is typically faster to write, less
error-prone and easier to debug.
The downside is that Python is harder to optimize β that is, turn into fast machine code β
than languages like C or Fortran.
Indeed, the standard implementation of Python (called CPython) cannot match the speed of
compiled languages such as C or Fortran.
Does that mean that we should just switch to C or Fortran for everything?
The answer is no, no and one hundred times no.
High productivity languages should be chosen over high-speed languages for the great major-
ity of scientific computing tasks.
This is because
1. Of any given program, relatively few lines are ever going to be time-critical.
2. For those lines of code that are time-critical, we can achieve C-like speed using a combi-
nation of NumPy and Numba.
Letβs start by trying to understand why high-level languages like Python are slower than com-
piled code.
20
[2]:
Even for this simple operation, the Python interpreter has a fair bit of work to do.
For example, in the statement a + b, the interpreter has to know which operation to invoke.
If a and b are strings, then a + b requires string concatenation
9.3. WHERE ARE THE BOTTLENECKS? 125
a, b = 'foo', 'bar'
[3]: a + b
'foobar'
[3]:
['foo', 'bar']
[4]:
(We say that the operator + is overloaded β its action depends on the type of the objects on
which it acts)
As a result, Python must check the type of the objects and then call the correct operation.
This involves substantial overheads.
Static Types
Compiled languages avoid these overheads with explicit, static types.
For example, consider the following C code, which sums the integers from 1 to 10
#include <stdio.h>
int main(void) {
int i;
int sum = 0;
for (i = 1; i <= 10; i++) {
sum = sum + i;
}
printf("sum = %d\n", sum);
return 0;
}
β’ In modern computers, memory addresses are allocated to each byte (one byte = 8 bits).
β’ For example, a 64 bit integer is stored in 8 bytes of memory.
β’ An array of π such integers occupies 8π consecutive memory slots.
126 CHAPTER 9. NUMBA
Moreover, the compiler is made aware of the data type by the programmer.
Hence, each successive data point can be accessed by shifting forward in memory space by a
known and fixed amount.
9.4 Vectorization
β’ The machine code itself is typically compiled from carefully optimized C or Fortran.
This can greatly accelerate many (but not all) numerical computations.
0.646047830581665
[6]:
0.022841215133666992
[7]:
The second code block β which achieves the same thing as the first β runs much faster.
The reason is that in the second implementation we have broken the loop down into three
basic operations
1. draw n uniforms
2. square them
3. sum them
Many functions provided by NumPy are so-called universal functions β also called ufuncs.
This means that they
For example, consider the problem of maximizing a function π of two variables (π₯, π¦) over the
square [βπ, π] Γ [βπ, π].
For π and π letβs choose
cos(π₯2 + π¦2 )
π(π₯, π¦) = and π = 3
1 + π₯2 + π¦ 2
Hereβs a plot of π
import matplotlib.pyplot as plt
[10]: %matplotlib inline
from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib import cm
qe.tic()
for x in grid:
for y in grid:
z = f(x, y)
if z > m:
m = z
qe.toc()
4.219594240188599
[11]:
qe.tic()
np.max(f(x, y))
qe.toc()
0.02581310272216797
[12]:
In the vectorized version, all the looping takes place in compiled code.
As you can see, the second version is much faster.
(Weβll make it even faster again below when we discuss Numba)
9.5 Numba
9.5.1 Prerequisites
9.5.2 An Example
π₯π‘+1 = 4π₯π‘ (1 β π₯π‘ )
Hereβs the plot of a typical trajectory, starting from π₯0 = 0.1, with π‘ on the x-axis
def qm(x0, n):
[13]: x = np.empty(n+1)
x[0] = x0
for t in range(n):
x[t+1] = 4 * x[t] * (1 - x[t])
return x
x = qm(0.1, 250)
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(x, 'b-', lw=2, alpha=0.8)
ax.set_xlabel('time', fontsize=16)
plt.show()
9.5. NUMBA 131
Letβs time and compare identical function calls across these two versions:
qe.util.tic()
[15]: qm(0.1, int(10**5))
time1 = qe.util.toc()
qe.util.tic()
[16]: qm_numba(0.1, int(10**5))
time2 = qe.util.toc()
The first execution is relatively slow because of JIT compilation (see below).
Next time and all subsequent times it runs much faster:
qe.util.tic()
[17]: qm_numba(0.1, int(10**5))
time2 = qe.util.toc()
Numba attempts to generate fast machine code using the infrastructure provided by the
LLVM Project.
It does this by inferring type information on the fly.
As you can imagine, this is easier for simple Python objects (simple scalar data types, such as
floats, integers, etc.).
Numba also plays well with NumPy arrays, which it treats as typed memory regions.
In an ideal setting, Numba can infer all necessary type information.
This allows it to generate native machine code, without having to call the Python runtime
environment.
In such a setting, Numba will be on par with machine code from low-level languages.
When Numba cannot infer all type information, some Python objects are given generic
object status, and some code is generated using the Python runtime.
In this second setting, Numba typically provides only minor speed gains β or none at all.
Hence, itβs prudent when using Numba to focus on speeding up small, time-critical snippets of
code.
This will give you much better performance than blanketing your Python programs with
@jit statements.
A Gotcha: Global Variables
Consider the following example
a = 1
[20]:
@jit
def add_x(x):
return a + x
print(add_x(10))
9.5. NUMBA 133
11
a = 2
[21]:
print(add_x(10))
11
Notice that changing the global had no effect on the value returned by the function.
When Numba compiles machine code for functions, it treats global variables as constants to
ensure type stability.
Numba can also be used to create custom ufuncs with the @vectorize decorator.
To illustrate the advantage of using Numba to vectorize a function, we return to a maximiza-
tion problem discussed above
from numba import vectorize
[22]:
@vectorize
def f_vec(x, y):
return np.cos(x**2 + y**2) / (1 + x**2 + y**2)
qe.tic()
np.max(f_vec(x, y))
qe.toc()
0.025510549545288086
[22]:
qe.tic()
np.max(f_vec(x, y))
qe.toc()
0.01980733871459961
[23]:
10.1 Contents
β’ Overview 10.2
β’ Cython 10.3
β’ Joblib 10.4
β’ Exercises 10.6
β’ Solutions 10.7
In addition to whatβs in Anaconda, this lecture will need the following libraries:
!pip install --upgrade quantecon
[1]:
10.2 Overview
In this lecture, we review some other scientific libraries that are useful for economic research
and analysis.
We have, however, already picked most of the low hanging fruit in terms of economic re-
search.
Hence you should feel free to skip this lecture on first pass.
10.3 Cython
Like Numba, Cython provides an approach to generating fast compiled code that can be used
from Python.
As was the case with Numba, a key problem is the fact that Python is dynamically typed.
As youβll recall, Numba solves this problem (where possible) by inferring type.
135
136 CHAPTER 10. OTHER SCIENTIFIC LIBRARIES
Cythonβs approach is different β programmers add type definitions directly to their βPythonβ
code.
As such, the Cython language can be thought of as Python with type definitions.
In addition to a language specification, Cython is also a language translator, transforming
Cython code into optimized C and C++ code.
Cython also takes care of building language extensions β the wrapper code that interfaces
between the resulting compiled code and Python.
Important Note:
In what follows code is executed in a Jupyter notebook.
This is to take advantage of a Cython cell magic that makes Cython particularly easy to use.
Some modifications are required to run the code outside a notebook.
π
1 β πΌπ+1
β πΌπ =
π=0
1βπΌ
return sum;
}
If youβre not familiar with C, the main thing you should take notice of is the type definitions
Here cdef is a Cython keyword indicating a variable declaration and is followed by a type.
The %%cython line at the top is not actually Cython code β itβs a Jupyter cell magic indi-
cating the start of Cython code.
After executing the cell, you can now call the function geo_prog_cython from within
Python.
What you are in fact calling is compiled C code with a Python call interface
import quantecon as qe
[5]: qe.util.tic()
geo_prog(0.99, int(10**6))
qe.util.toc()
0.17586088180541992
[5]:
qe.util.tic()
[6]: geo_prog_cython(0.99, int(10**6))
qe.util.toc()
0.04195237159729004
[6]:
138 CHAPTER 10. OTHER SCIENTIFIC LIBRARIES
Letβs go back to the first problem that we worked with: generating the iterates of the
quadratic map
π₯π‘+1 = 4π₯π‘ (1 β π₯π‘ )
The problem of computing iterates and returning a time series requires us to work with ar-
rays.
The natural array type to work with is NumPy arrays.
Hereβs a Cython implementation that initializes, populates and returns a NumPy array
%%cython
[7]: import numpy as np
If you run this code and time it, you will see that its performance is disappointing β nothing
like the speed gain we got from Numba
qe.util.tic()
[8]: qm_cython_first_pass(0.1, int(10**5))
qe.util.toc()
0.06516456604003906
[8]:
This example was also computed in the Numba lecture, and you can see Numba is around 90
times faster.
The reason is that working with NumPy arrays incurs substantial Python overheads.
We can do better by using Cythonβs typed memoryviews, which provide more direct access to
arrays in memory.
When using them, the first step is to create a NumPy array.
Next, we declare a memoryview and bind it to the NumPy array.
Hereβs an example:
%%cython
[9]: import numpy as np
from numpy cimport float_t
Here
0.0008153915405273438
[10]:
10.3.3 Summary
Cython requires more expertise than Numba, and is a little more fiddly in terms of getting
good performance.
In fact, itβs surprising how difficult it is to beat the speed improvements provided by Numba.
Nonetheless,
10.4 Joblib
10.4.1 Caching
Perhaps, like us, you sometimes run a long computation that simulates a model at a given set
of parameters β to generate a figure, say, or a table.
140 CHAPTER 10. OTHER SCIENTIFIC LIBRARIES
20 minutes later you realize that you want to tweak the figure and now you have to do it all
again.
What caching will do is automatically store results at each parameterization.
With Joblib, results are compressed and stored on file, and automatically served back up to
you when you repeat the calculation.
10.4.2 An Example
Letβs look at a toy example, related to the quadratic map model discussed above.
Letβs say we want to generate a long trajectory from a certain initial condition π₯0 and see
what fraction of the sample is below 0.1.
(Weβll omit JIT compilation or other speedups for simplicity)
Hereβs our code
from joblib import Memory
[12]: location = './cachedir'
memory = Memory(location='./joblib_cache')
@memory.cache
def qm(x0, n):
x = np.empty(n+1)
x[0] = x0
for t in range(n):
x[t+1] = 4 * x[t] * (1 - x[t])
return np.mean(x < 0.1)
We are using joblib to cache the result of calling qm at a given set of parameters.
With the argument location=β./joblib_cacheβ, any call to this function results in both the in-
put values and output values being stored a subdirectory joblib_cache of the present working
directory.
(In UNIX shells, . refers to the present working directory)
The first time we call the function with a given set of parameters we see some extra output
that notes information being cached
qe.util.tic()
[13]: n = int(1e7)
qm(0.2, n)
qe.util.toc()
________________________________________________________________________________
[Memory] Calling __main__--home-ubuntu-repos-lecture-source-py-_build-pdf-
jupyter-executed-__ipython-input__.qmβ¦
qm(0.2, 10000000)
______________________________________________________________qm - 16.9s, 0.3min
TOC: Elapsed: 0:00:16.87
16.87077283859253
[13]:
The next time we call the function with the same set of parameters, the result is returned
almost instantaneously
qe.util.tic()
[14]: n = int(1e7)
qm(0.2, n)
qe.util.toc()
10.5. OTHER OPTIONS 141
0.001150369644165039
[14]:
There are in fact many other approaches to speeding up your Python code.
One is interfacing with Fortran.
If you are comfortable writing Fortran you will find it very easy to create extension modules
from Fortran code using F2Py.
F2Py is a Fortran-to-Python interface generator that is particularly simple to use.
Robert Johansson provides a very nice introduction to F2Py, among other things.
Recently, a Jupyter cell magic for Fortran has been developed β you might want to give it a
try.
10.6 Exercises
10.6.1 Exercise 1
For example, let the period length be one month, and suppose the current state is high.
We see from the graph that the state next month will be
Your task is to simulate a sequence of monthly volatility states according to this rule.
Set the length of the sequence to n = 100000 and start in the high state.
Implement a pure Python version, a Numba version and a Cython version, and compare
speeds.
142 CHAPTER 10. OTHER SCIENTIFIC LIBRARIES
To test your code, evaluate the fraction of time that the chain spends in the low state.
If your code is correct, it should be about 2/3.
10.7 Solutions
10.7.1 Exercise 1
We let
β’ 0 represent βlowβ
β’ 1 represent βhighβ
Letβs run this code and check that the fraction of time spent in the low state is about 0.666
n = 100000
[17]: x = compute_series(n)
print(np.mean(x == 0)) # Fraction of time x is in state 0
0.66673
0.14121103286743164
[18]:
0.66615
0.0012984275817871094
[21]:
%%cython
[23]: import numpy as np
from numpy cimport int_t, float_t
compute_series_cy(10)
[24]:
array([1, 1, 1, 1, 1, 0, 0, 0, 1, 0])
[24]:
x = compute_series_cy(n)
[25]: print(np.mean(x == 0))
0.65561
qe.util.tic()
[26]: compute_series_cy(n)
qe.util.toc()
0.0029497146606445312
[26]:
145
Chapter 11
11.1 Contents
β’ Overview 11.2
β’ Summary 11.6
11.2 Overview
When computer programs are small, poorly written code is not overly costly.
But more data, more sophisticated models, and more computer power are enabling us to take
on more challenging problems that involve writing longer programs.
For such programs, investment in good coding practices will pay high returns.
The main payoffs are higher productivity and faster code.
In this lecture, we review some elements of good coding practice.
We also touch on modern developments in scientific computing β such as just in time compi-
lation β and how they affect good program design.
Here
147
148 CHAPTER 11. WRITING GOOD CODE
1. sets π0 = 1
2. iterates using Eq. (1) to produce a sequence π0 , π1 , π2 β¦ , ππ
3. plots the sequence
for j in range(3):
k[0] = 1
for t in range(49):
k[t+1] = s * k[t]**Ξ±[j] + (1 - Ξ΄) * k[t]
axes[0].plot(k, 'o-', label=rf"$\alpha = {Ξ±[j]},\; s = {s},\; \delta={Ξ΄}$")
axes[0].grid(lw=0.2)
axes[0].set_ylim(0, 18)
axes[0].set_xlabel('time')
axes[0].set_ylabel('capital')
axes[0].legend(loc='upper left', frameon=True, fontsize=14)
for j in range(3):
k[0] = 1
for t in range(49):
k[t+1] = s[j] * k[t]**Ξ± + (1 - Ξ΄) * k[t]
axes[1].plot(k, 'o-', label=rf"$\alpha = {Ξ±},\; s = {s},\; \delta={Ξ΄}$")
axes[1].grid(lw=0.2)
axes[1].set_xlabel('time')
axes[1].set_ylabel('capital')
axes[1].set_ylim(0, 18)
axes[1].legend(loc='upper left', frameon=True, fontsize=14)
for j in range(3):
k[0] = 1
for t in range(49):
k[t+1] = s * k[t]**Ξ± + (1 - Ξ΄[j]) * k[t]
axes[2].plot(k, 'o-', label=rf"$\alpha = {Ξ±},\; s = {s},\; \delta={Ξ΄[j]}$")
axes[2].set_ylim(0, 18)
11.3. AN EXAMPLE OF BAD CODE 149
axes[2].set_xlabel('time')
axes[2].set_ylabel('capital')
axes[2].grid(lw=0.2)
axes[2].legend(loc='upper left', frameon=True, fontsize=14)
plt.show()
There are usually many different ways to write a program that accomplishes a given task.
For small programs, like the one above, the way you write code doesnβt matter too much.
But if you are ambitious and want to produce useful things, youβll write medium to large pro-
grams too.
In those settings, coding style matters a great deal.
Fortunately, lots of smart people have thought about the best way to write code.
Here are some basic precepts.
If you look at the code above, youβll see numbers like 50 and 49 and 3 scattered through the
code.
These kinds of numeric literals in the body of your code are sometimes called βmagic num-
bersβ.
This is not a complement.
While numeric literals are not all evil, the numbers shown in the program above should cer-
tainly be replaced by named constants.
For example, the code above could declare the variable time_series_length = 50.
Then in the loops, 49 should be replaced by time_series_length - 1.
The advantages are:
Yes, we realize that you can just cut and paste and change a few symbols.
But as a programmer, your aim should be to automate repetition, not do it yourself.
More importantly, repeating the same logic in different places means that eventually one of
them will likely be wrong.
If you want to know more, read the excellent summary found on this page.
Weβll talk about how to avoid repetition below.
11.4. GOOD CODING PRACTICE 151
Sure, global variables (i.e., names assigned to values outside of any function or class) are con-
venient.
Rookie programmers typically use global variables with abandon β as we once did ourselves.
But global variables are dangerous, especially in medium to large size programs, since
This makes it much harder to be certain about what some small part of a given piece of code
actually commands.
Hereβs a useful discussion on the topic.
While the odd global in small scripts is no big deal, we recommend that you teach yourself to
avoid them.
(Weβll discuss how just below).
JIT Compilation
In fact, thereβs now another good reason to avoid global variables.
In scientific computing, weβre witnessing the rapid growth of just in time (JIT) compilation.
JIT compilation can generate excellent performance for scripting languages like Python and
Julia.
But the task of the compiler used for JIT compilation becomes much harder when many
global variables are present.
(This is because data type instability hinders the generation of efficient machine code β weβll
learn more about such topics later on)
Fortunately, we can easily avoid the evils of global variables and WET code.
β’ WET stands for βwe love typingβ and is the opposite of DRY.
Hereβs some code that reproduces the plot above with better coding style.
It uses a function to avoid repetition.
Note also that
β’ global variables are quarantined by collecting together at the end, not the start of the
program
β’ magic numbers are avoided
β’ the loop at the end where the actual work is done is short and relatively simple
ax.grid(lw=0.2)
ax.set_xlabel('time')
ax.set_ylabel('capital')
ax.set_ylim(0, 18)
ax.legend(loc='upper left', frameon=True, fontsize=14)
plt.show()
11.6. SUMMARY 153
11.6 Summary
12.1 Contents
β’ Overview 12.2
β’ Exercises 12.6
β’ Solutions 12.7
12.2 Overview
So imagine now you want to write a program with consumers, who can
155
156 CHAPTER 12. OOP II: BUILDING CLASSES
As discussed an earlier lecture, in the OOP paradigm, data and functions are bundled to-
gether into βobjectsβ.
An example is a Python list, which not only stores data but also knows how to sort itself, etc.
x = [1, 5, 4]
[2]: x.sort()
x
[1, 4, 5]
[2]:
As we now know, sort is a function that is βpart ofβ the list object β and hence called a
method.
If we want to make our own types of objects we need to use class definitions.
A class definition is a blueprint for a particular class of objects (e.g., lists, strings or complex
numbers).
It describes
β’ Methods set out in the class definition act on this (and other) data.
In Python, the data and methods of an object are collectively referred to as attributes.
Attributes are accessed via βdotted attribute notationβ
β’ object_name.data
β’ object_name.method_name()
In the example
x = [1, 5, 4]
[3]: x.sort()
x.__class__
list
[3]:
β’ x is an object or instance, created from the definition for Python lists, but with its own
particular data.
β’ x.sort() and x.__class__ are two attributes of x.
β’ dir(x) can be used to view all the attributes of x.
OOP is useful for the same reason that abstraction is useful: for recognizing and exploiting
the common structure.
For example,
β’ a Markov chain consists of a set of states and a collection of transition probabilities for
moving across states
β’ a general equilibrium theory consists of a commodity space, preferences, technologies,
and an equilibrium definition
β’ a game consists of a list of players, lists of actions available to each player, player pay-
offs as functions of all playersβ actions, and a timing protocol
These are all abstractions that collect together βobjectsβ of the same βtypeβ.
Recognizing common structure allows us to employ common tools.
In economic theory, this might be a proposition that applies to all games of a certain type.
In Python, this might be a method thatβs useful for all Markov chains (e.g., simulate).
When we use OOP, the simulate method is conveniently bundled together with the Markov
chain object.
Admittedly a little contrived, this example of a class helps us internalize some new syntax.
Hereβs one implementation
class Consumer:
[4]:
def __init__(self, w):
"Initialize consumer with w dollars of wealth"
self.wealth = w
This class defines instance data wealth and three methods: __init__, earn and spend
β’ wealth is instance data because each consumer we create (each instance of the
Consumer class) will have its own separate wealth data.
The ideas behind the earn and spend methods were discussed above.
Both of these act on the instance data wealth.
The __init__ method is a constructor method.
Whenever we create an instance of the class, this method will be called automatically.
Calling __init__ sets up a βnamespaceβ to hold the instance data β more on this soon.
Weβll also discuss the role of self just below.
Usage
Hereβs an example of usage
c1 = Consumer(10) # Create instance with initial wealth 10
[5]: c1.spend(5)
c1.wealth
5
[5]:
12.4. DEFINING YOUR OWN CLASSES 159
c1.earn(15)
[6]: c1.spend(100)
Insufficent funds
We can of course create multiple instances each with its own data
c1 = Consumer(10)
[7]: c2 = Consumer(12)
c2.spend(4)
c2.wealth
8
[7]:
c1.wealth
[8]:
10
[8]:
When we access or set attributes weβre actually just modifying the dictionary maintained by
the instance.
Self
If you look at the Consumer class definition again youβll see the word self throughout the
code.
The rules with self are that
β e.g., the earn method references self.wealth rather than just wealth
β’ Any method defined within the class should have self as its first argument
There are no examples of the last rule in the preceding code but we will see some shortly.
Details
In this section, we look at some more formal details related to classes and self
β’ You might wish to skip to the next section on first pass of this lecture.
β’ You can return to these details after youβve familiarized yourself with more examples.
Methods actually live inside a class object formed when the interpreter reads the class defini-
tion
160 CHAPTER 12. OOP II: BUILDING CLASSES
Note how the three methods __init__, earn and spend are stored in the class object.
Consider the following code
c1 = Consumer(10)
[12]: c1.earn(10)
c1.wealth
20
[12]:
When you call earn via c1.earn(10) the interpreter passes the instance c1 and the argu-
ment 10 to Consumer.earn.
In fact, the following are equivalent
β’ c1.earn(10)
β’ Consumer.earn(c1, 10)
In the function call Consumer.earn(c1, 10) note that c1 is the first argument.
Recall that in the definition of the earn method, self is the first parameter
def earn(self, y):
[13]: "The consumer earns y dollars"
self.wealth += y
The end result is that self is bound to the instance c1 inside the function call.
Thatβs why the statement self.wealth += y inside earn ends up modifying c1.wealth.
For our next example, letβs write a simple class to implement the Solow growth model.
The Solow growth model is a neoclassical growth model where the amount of capital stock
per capita ππ‘ evolves according to the rule
π π§ππ‘πΌ + (1 β πΏ)ππ‘
ππ‘+1 = (1)
1+π
Here
The steady state of the model is the π that solves Eq. (1) when ππ‘+1 = ππ‘ = π.
Hereβs a class that implements this model.
Some points of interest in the code are
β’ An instance maintains a record of its current capital stock in the variable self.k.
β Notice how inside update the reference to the local method h is self.h.
"""
def __init__(self, n=0.05, # population growth rate
s=0.25, # savings rate
Ξ΄=0.1, # depreciation rate
Ξ±=0.3, # share of labor
z=2.0, # productivity
k=1.0): # current capital stock
def h(self):
"Evaluate the h function"
# Unpack parameters (get rid of self to simplify notation)
n, s, Ξ΄, Ξ±, z = self.n, self.s, self.Ξ΄, self.Ξ±, self.z
# Apply the update rule
return (s * z * self.k**Ξ± + (1 - Ξ΄) * self.k) / (1 + n)
def update(self):
"Update the current state (i.e., the capital stock)."
self.k = self.h()
def steady_state(self):
"Compute the steady state value of capital."
# Unpack parameters (get rid of self to simplify notation)
n, s, Ξ΄, Ξ±, z = self.n, self.s, self.Ξ΄, self.Ξ±, self.z
# Compute and return steady state
return ((s * z) / (n + Ξ΄))**(1 / (1 - Ξ±))
Hereβs a little program that uses the class to compute time series from two different initial
conditions.
The common steady state is also plotted for comparison
s1 = Solow()
[15]: s2 = Solow(k=8.0)
T =
162 CHAPTER 12. OOP II: BUILDING CLASSES
ax.legend()
plt.show()
Next, letβs write a class for a simple one good market where agents are price takers.
The market consists of the following objects:
Here
The class provides methods to compute various values of interest, including competitive equi-
librium price and quantity, tax revenue raised, consumer surplus and producer surplus.
Hereβs our implementation
12.4. DEFINING YOUR OWN CLASSES 163
"""
self.ad, self.bd, self.az, self.bz, self.tax = ad, bd, az, bz, tax
if ad < az:
raise ValueError('Insufficient demand.')
def price(self):
"Return equilibrium price"
return (self.ad - self.az + self.bz * self.tax) / (self.bd + self.bz)
def quantity(self):
"Compute equilibrium quantity"
return self.ad - self.bd * self.price()
def consumer_surp(self):
"Compute consumer surplus"
# == Compute area under inverse demand function == #
integrand = lambda x: (self.ad / self.bd) - (1 / self.bd) * x
area, error = quad(integrand, 0, self.quantity())
return area - self.price() * self.quantity()
def producer_surp(self):
"Compute producer surplus"
# == Compute area above inverse supply curve, excluding tax == #
integrand = lambda x: -(self.az / self.bz) + (1 / self.bz) * x
area, error = quad(integrand, 0, self.quantity())
return (self.price() - self.tax) * self.quantity() - area
def taxrev(self):
"Compute tax revenue"
return self.tax * self.quantity()
Hereβs a short program that uses this class to plot an inverse demand curve together with in-
verse supply curves with and without taxes
164 CHAPTER 12. OOP II: BUILDING CLASSES
q_max = m.quantity() * 2
q_grid = np.linspace(0.0, q_max, 100)
pd = m.inverse_demand(q_grid)
ps = m.inverse_supply(q_grid)
psno = m.inverse_supply_no_tax(q_grid)
fig, ax = plt.subplots()
ax.plot(q_grid, pd, lw=2, alpha=0.6, label='demand')
ax.plot(q_grid, ps, lw=2, alpha=0.6, label='supply')
ax.plot(q_grid, psno, '--k', lw=2, alpha=0.6, label='supply without tax')
ax.set_xlabel('quantity', fontsize=14)
ax.set_xlim(0, q_max)
ax.set_ylabel('price', fontsize=14)
ax.legend(loc='lower right', frameon=False, fontsize=14)
plt.show()
def deadw(m):
[20]: "Computes deadweight loss for market m."
# == Create analogous market with no tax == #
m_no_tax = Market(m.ad, m.bd, m.az, m.bz, 0)
# == Compare surplus, return difference == #
surp1 = m_no_tax.consumer_surp() + m_no_tax.producer_surp()
surp2 = m.consumer_surp() + m.producer_surp() + m.taxrev()
return surp1 - surp2
1.125
[21]:
Letβs look at one more example, related to chaotic dynamics in nonlinear systems.
One simple transition rule that can generate complex dynamics is the logistic map
Letβs write a class for generating time series from this model.
Hereβs one implementation
class Chaos:
[22]: """
Models the dynamical system with :math:`x_{t+1} = r x_t (1 - x_t)`
"""
def __init__(self, x0, r):
"""
Initialize with state x0 and parameter r
"""
self.x, self.r = x0, r
def update(self):
"Apply the map to update state."
self.x = self.r * self.x *(1 - self.x)
fig, ax = plt.subplots()
ax.set_xlabel('$t$', fontsize=14)
ax.set_ylabel('$x_t$', fontsize=14)
x = ch.generate_sequence(ts_length)
ax.plot(range(ts_length), x, 'bo-', alpha=0.5, lw=2, label='$x_t$')
plt.show()
166 CHAPTER 12. OOP II: BUILDING CLASSES
ax.set_xlabel('$r$', fontsize=16)
plt.show()
12.5. SPECIAL METHODS 167
Python provides special methods with which some neat tricks can be performed.
For example, recall that lists and tuples have a notion of length and that this length can be
queried via the len function
x = (10, 20)
[26]: len(x)
2
[26]:
168 CHAPTER 12. OOP II: BUILDING CLASSES
If you want to provide a return value for the len function when applied to your user-defined
object, use the __len__ special method
class Foo:
[27]:
def __len__(self):
return 42
Now we get
f = Foo()
[28]: len(f)
42
[28]:
50
[30]:
12.6 Exercises
12.6.1 Exercise 1
The empirical cumulative distribution function (ecdf) corresponding to a sample {ππ }ππ=1 is
defined as
1 π
πΉπ (π₯) βΆ= β 1{ππ β€ π₯} (π₯ β R) (3)
π π=1
Here 1{ππ β€ π₯} is an indicator function (one if ππ β€ π₯ and zero otherwise) and hence πΉπ (π₯)
is the fraction of the sample that falls below π₯.
The GlivenkoβCantelli Theorem states that, provided that the sample is IID, the ecdf πΉπ con-
verges to the true distribution function πΉ .
Implement πΉπ as a class called ECDF, where
β’ A given sample {ππ }ππ=1 are the instance data, stored as self.observations.
β’ The class implements a __call__ method that returns πΉπ (π₯) for any π₯.
12.6.2 Exercise 2
π
2
π(π₯) = π0 + π1 π₯ + π2 π₯ + β― ππ π₯ π
= β ππ π₯π (π₯ β R) (4)
π=0
The instance data for the class Polynomial will be the coefficients (in the case of Eq. (4),
the numbers π0 , β¦ , ππ ).
Provide methods that
12.7 Solutions
12.7.1 Exercise 1
class ECDF:
[31]:
def __init__(self, observations):
self.observations = observations
# == test == #
[32]:
from random import uniform
print(F(0.5))
0.4
0.48
12.7.2 Exercise 2
class Polynomial:
[33]:
def __init__(self, coefficients):
"""
Creates an instance of the Polynomial class representing
def differentiate(self):
"Reset self.coefficients to those of p' instead of p."
new_coefficients = []
for i, a in enumerate(self.coefficients):
new_coefficients.append(i * a)
# Remove the first element, which is zero
del new_coefficients[0]
# And reset coefficients data to new values
self.coefficients = new_coefficients
return new_coefficients
Chapter 13
13.1 Contents
β’ Overview 13.2
β’ Details 13.3
β’ Implementation 13.4
β’ Summary 13.10
13.2 Overview
This lecture creates non-stochastic and stochastic versions of Paul Samuelsonβs celebrated
multiplier accelerator model [118].
In doing so, we extend the example of the Solow model class in our second OOP lecture.
Our objectives are to
171
172 CHAPTER 13. OOP III: SAMUELSON MULTIPLIER ACCELERATOR
Weβll also use the following for various tasks described below:
from quantecon import LinearStateSpace
[3]: import cmath
import math
import sympy
from sympy import Symbol, init_printing
from cmath import sqrt
Samuelson used a second-order linear difference equation to represent a model of national out-
put based on three components:
β’ a national output identity asserting that national outcome is the sum of consumption
plus investment plus government purchases.
β’ a Keynesian consumption function asserting that consumption at time π‘ is equal to a
constant times national output at time π‘ β 1.
β’ an investment accelerator asserting that investment at time π‘ equals a constant called
the accelerator coefficient times the difference in output between period π‘ β 1 and π‘ β 2.
β’ the idea that consumption plus investment plus government purchases constitute aggre-
gate demand, which automatically calls forth an equal amount of aggregate supply.
(To read about linear difference equations see here or chapter IX of [121])
Samuelson used the model to analyze how particular values of the marginal propensity to
consume and the accelerator coefficient might give rise to transient business cycles in national
output.
Possible dynamic properties include
Later we present an extension that adds a random shock to the right side of the national in-
come identity representing random fluctuations in aggregate demand.
This modification makes national output become governed by a second-order stochastic linear
difference equation that, with appropriate parameter values, gives rise to recurrent irregular
business cycles.
(To read about stochastic linear difference equations see chapter XI of [121])
13.3 Details
πΆπ‘ = πππ‘β1 + πΎ (1)
ππ‘ = πΆπ‘ + πΌπ‘ + πΊπ‘ (3)
Equations Eq. (1), Eq. (2), and Eq. (3) imply the following second-order linear difference
equation for national income:
ππ‘ = (π + π)ππ‘β1 β πππ‘β2 + (πΎ + πΊπ‘ )
or
Μ ,
πβ1 = πβ1 Μ
πβ2 = πβ2
Weβll ordinarily set the parameters (π, π) so that starting from an arbitrary pair of initial con-
Μ , πβ2
ditions (πβ1 Μ ), national income π _π‘ converges to a constant value as π‘ becomes large.
The deterministic version of the model described so far β meaning that no random shocks
hit aggregate demand β has only transient fluctuations.
We can convert the model to one that has persistent irregular fluctuations by adding a ran-
dom shock to aggregate demand.
ππ‘ = π1 ππ‘β1 + π2 ππ‘β2
or
To discover the properties of the solution of Eq. (6), it is useful first to form the characteris-
tic polynomial for Eq. (6):
π§ 2 β π1 π§ β π 2 (7)
These are two special values of π§, say π§ = π1 and π§ = π2 , such that if we set π§ equal to one of
these values in expression Eq. (7), the characteristic polynomial Eq. (7) equals zero:
π§2 β π1 π§ β π2 = (π§ β π1 )(π§ β π2 ) = 0 (8)
π1 = ππππ , π2 = ππβππ
where π is the amplitude of the complex number and π is its angle or phase.
These can also be represented as
π1 = π(πππ (π) + π sin(π))
π2 = π(πππ (π) β π sin(π))
ππ‘ = ππ‘1 π1 + ππ‘2 π2
where π1 and π2 are constants that depend on the two initial conditions and on π1 , π2 .
When the roots are complex, it is useful to pursue the following calculations.
Notice that
ππ‘ = π1 (ππππ )π‘ + π2 (ππβππ )π‘
= π1 ππ‘ ππππ‘ + π2 ππ‘ πβπππ‘
= π1 ππ‘ [cos(ππ‘) + π sin(ππ‘)] + π2 ππ‘ [cos(ππ‘) β π sin(ππ‘)]
= (π1 + π2 )ππ‘ cos(ππ‘) + π(π1 β π2 )ππ‘ sin(ππ‘)
The only way that ππ‘ can be a real number for each π‘ is if π1 + π2 is a real number and π1 β π2
is an imaginary number.
This happens only when π1 and π2 are complex conjugates, in which case they can be written
in the polar forms
π1 = π£πππ , π2 = π£πβππ
So we can write
176 CHAPTER 13. OOP III: SAMUELSON MULTIPLIER ACCELERATOR
where π£ and π are constants that must be chosen to satisfy initial conditions for πβ1 , πβ2 .
This formula shows that when the roots are complex, ππ‘ displays oscillations with period
πΜ = 2π
π and damping factor π.
We say that πΜ is the period because in that amount of time the cosine wave cos(ππ‘ + π) goes
through exactly one complete cycles.
(Draw a cosine function to convince yourself of this please)
Remark: Following [118], we want to choose the parameters π, π of the model so that the ab-
solute values (of the possibly complex) roots π1 , π2 of the characteristic polynomial are both
strictly less than one:
Remark: When both roots π1 , π2 of the characteristic polynomial have absolute values
strictly less than one, the absolute value of the larger one governs the rate of convergence to
the steady state of the non stochastic version of the model.
We want to have a method in the class that automatically generates a simulation, either non-
stochastic (π = 0) or stochastic (π > 0).
We also show how to map the Samuelson model into a simple instance of the
LinearStateSpace class described here.
We can use a LinearStateSpace instance to do various things that we did above with our
homemade function and class.
Among other things, we show by example that the eigenvalues of the matrix π΄ that we use to
form the instance of the LinearStateSpace class for the Samuelson model equal the roots
of the characteristic polynomial Eq. (7) for the Samuelson multiplier accelerator model.
Here is the formula for the matrix π΄ in the linear state space system in the case that govern-
ment expenditures are a constant πΊ:
1 0 0
π΄=β‘
β’πΎ + πΊ π β€
1 π2 β₯
β£ 0 1 0β¦
13.4 Implementation
# Set axis
xmin, ymin = -3, -2
xmax, ymax = -xmin, -ymin
plt.axis([xmin, xmax, ymin, ymax])
arrowprops=plot_arrowsr, fontsize='12')
ax.annotate(r'$\rho_1 + \rho_2 = 1$', xy=(0.38, 0.6), xytext=(0.6, 0.8),
arrowprops=plot_arrowsr, fontsize='12')
ax.annotate(r'$\rho_2 < 1 + \rho_1$', xy=(-0.5, 0.3), xytext=(-1.3, 0.6),
arrowprops=plot_arrowsl, fontsize='12')
ax.annotate(r'$\rho_2 = 1 + \rho_1$', xy=(-0.38, 0.6), xytext=(-1, 0.8),
arrowprops=plot_arrowsl, fontsize='12')
ax.annotate(r'$\rho_2 = -1$', xy=(1.5, -1), xytext=(1.8, -1.3),
arrowprops=plot_arrowsl, fontsize='12')
ax.annotate(r'${\rho_1}^2 + 4\rho_2 = 0$', xy=(1.15, -0.35),
xytext=(1.5, -0.3), arrowprops=plot_arrowsr, fontsize='12')
ax.annotate(r'${\rho_1}^2 + 4\rho_2 < 0$', xy=(1.4, -0.7),
xytext=(1.8, -0.6), arrowprops=plot_arrowsr, fontsize='12')
return fig
param_plot()
plt.show()
The graph portrays regions in which the (π1 , π2 ) root pairs implied by the (π1 = (π + π), π2 =
βπ) difference equation parameter pairs in the Samuelson model are such that:
β’ (π1 , π2 ) are complex with modulus less than 1 - in this case, the {ππ‘ } sequence displays
damped oscillations.
β’ (π1 , π2 ) are both real, but one is strictly greater than 1 - this leads to explosive growth.
13.4. IMPLEMENTATION 179
β’ (π1 , π2 ) are both real, but one is strictly less than β1 - this leads to explosive oscilla-
tions.
β’ (π1 , π2 ) are both real and both are less than 1 in absolute value - in this case, there is
smooth convergence to the steady state without damped cycles.
Later weβll present the graph with a red mark showing the particular point implied by the
setting of (π, π).
discriminant = Ο1 ** 2 + 4 * Ο2
if Ο2 > 1 + Ο1 or Ο2 < -1:
print('Explosive oscillations')
elif Ο1 + Ο2 > 1:
print('Explosive growth')
elif discriminant < 0:
print('Roots are complex with modulus less than one; \
therefore damped oscillations')
else:
print('Roots are real and absolute values are less than one; \
therefore get smooth convergence to a steady state')
Roots are real and absolute values are less than one; therefore get smooth
convergence to a steady state
plt.subplots(figsize=(10, 6))
plt.plot(function)
plt.xlabel('Time $t$')
plt.ylabel('$Y_t$', rotation=0)
plt.grid()
plt.show()
The following function calculates roots of the characteristic polynomial using high school al-
gebra.
(Weβll calculate the roots in other ways later)
The function also plots a ππ‘ starting from initial conditions that we set
180 CHAPTER 13. OOP III: SAMUELSON MULTIPLIER ACCELERATOR
roots = []
Ο1 = Ξ± + Ξ²
Ο2 = -Ξ²
print(f'Ο_1 is {Ο1}')
print(f'Ο_2 is {Ο2}')
discriminant = Ο1 ** 2 + 4 * Ο2
if discriminant == 0:
roots.append(-Ο1 / 2)
print('Single real root: ')
print(''.join(str(roots)))
elif discriminant > 0:
roots.append((-Ο1 + sqrt(discriminant).real) / 2)
roots.append((-Ο1 - sqrt(discriminant).real) / 2)
print('Two real roots: ')
print(''.join(str(roots)))
else:
roots.append((-Ο1 + sqrt(discriminant)) / 2)
roots.append((-Ο1 - sqrt(discriminant)) / 2)
print('Two complex roots: ')
print(''.join(str(roots)))
return y_t
plot_y(y_nonstochastic())
Ο_1 is 1.42
Ο_2 is -0.5
Two real roots:
[-0.6459687576256715, -0.7740312423743284]
Absolute values of roots are less than one
13.4. IMPLEMENTATION 181
The next cell writes code that takes as inputs the modulus π and phase π of a conjugate pair
of complex numbers in polar form
π1 = π exp(ππ), π2 = π exp(βππ)
β’ The code assumes that these two complex numbers are the roots of the characteristic
polynomial
β’ It then reverse-engineers (π, π) and (π1 , π2 ), pairs that would generate those roots
r = .95
period = 10 # Length of cycle in units of time
οΏ½ = 2 * math.pi/period
a, b = (0.6346322893124001+0j), (0.9024999999999999-0j)
Ο1, Ο2 = (1.5371322893124+0j), (-0.9024999999999999+0j)
Ο1, Ο2
(1.5371322893124, -0.9024999999999999)
[10]:
Here weβll use numpy to compute the roots of the characteristic polynomial
r1, r2 = np.roots([1, -Ο1, -Ο2])
[11]:
p1 = cmath.polar(r1)
p2 = cmath.polar(r2)
r, οΏ½ = 0.95, 0.6283185307179586
p1, p2 = (0.95, 0.6283185307179586), (0.95, -0.6283185307179586)
a, b = (0.6346322893124001+0j), (0.9024999999999999-0j)
Ο1, Ο2 = 1.5371322893124, -0.9024999999999999
# Useful constants
Ο1 = Ξ± + Ξ²
Ο2 = -Ξ²
categorize_solution(Ο1, Ο2)
return y_t
plot_y(y_nonstochastic())
Roots are complex with modulus less than one; therefore damped oscillations
Roots are [0.85+0.27838822j 0.85-0.27838822j]
Roots are complex
Roots are less than one
a, b = 0.6180339887498949, 1.0
Roots are complex with modulus less than one; therefore damped oscillations
Roots are [0.80901699+0.58778525j 0.80901699-0.58778525j]
Roots are complex
Roots are less than one
We can also use sympy to compute analytic formulas for the roots
init_printing()
[14]:
r1 = Symbol("Ο_1")
r2 = Symbol("Ο_2")
z = Symbol("z")
[14]:
13.5. STOCHASTIC SHOCKS 185
π1 1 π1 1
[ β βπ12 + 4π2 , + βπ12 + 4π2 ]
2 2 2 2
a = Symbol("Ξ±")
[15]: b = Symbol("Ξ²")
r1 = a + b
r2 = -b
[15]:
πΌ π½ βπΌ2 + 2πΌπ½ + π½ 2 β 4π½ πΌ π½ βπΌ2 + 2πΌπ½ + π½ 2 β 4π½
[ + β , + + ]
2 2 2 2 2 2
πΌ π½ 1 πΌ π½ 1
[ + β βπΌ2 + 2πΌπ½ + π½ 2 β 4π½, + + βπΌ2 + 2πΌπ½ + π½ 2 β 4π½]
2 2 2 2 2 2
Now weβll construct some code to simulate the stochastic version of the model that emerges
when we add a random shock process to aggregate demand
def y_stochastic(y_0=0, y_1=0, Ξ±=0.8, Ξ²=0.2, Ξ³=10, n=100, Ο=5):
[16]:
"""This function takes parameters of a stochastic version of
the model and proceeds to analyze the roots of the characteristic
polynomial and also generate a simulation.
"""
# Useful constants
Ο1 = Ξ± + Ξ²
Ο2 = -Ξ²
# Categorize solution
categorize_solution(Ο1, Ο2)
# Generate shocks
οΏ½ = np.random.normal(0, 1, n)
return y_t
plot_y(y_stochastic())
Roots are real and absolute values are less than one; therefore get smooth
convergence to a steady state
[0.7236068 0.2763932]
Roots are real
Roots are less than one
Letβs do a simulation in which there are shocks and the characteristic polynomial has complex
roots
r = .97
[17]:
period = 10 # Length of cycle in units of time
οΏ½ = 2 * math.pi/period
a, b = 0.6285929690873979, 0.9409000000000001
Roots are complex with modulus less than one; therefore damped oscillations
[0.78474648+0.57015169j 0.78474648-0.57015169j]
Roots are complex
Roots are less than one
13.6. GOVERNMENT SPENDING 187
# Useful constants
Ο1 = Ξ± + Ξ²
Ο2 = -Ξ²
# Categorize solution
categorize_solution(Ο1, Ο2)
# Generate shocks
οΏ½ = np.random.normal(0, 1, n)
# Stochastic
else:
οΏ½ = np.random.normal(0, 1, n)
return Ο1 * x[t - 1] + Ο2 * x[t - 2] + Ξ³ + g + Ο * οΏ½[t]
# No government spending
if g == 0:
y_t.append(transition(y_t, t))
Roots are real and absolute values are less than one; therefore get smooth
convergence to a steady state
[0.7236068 0.2763932]
Roots are real
Roots are less than one
13.6. GOVERNMENT SPENDING 189
We can also see the response to a one time jump in government expenditures
plot_y(y_stochastic_g(g=500, g_t=50, duration='one-off'))
[20]:
Roots are real and absolute values are less than one; therefore get smooth
convergence to a steady state
[0.7236068 0.2763932]
Roots are real
Roots are less than one
190 CHAPTER 13. OOP III: SAMUELSON MULTIPLIER ACCELERATOR
.. math::
Parameters
----------
y_0 : scalar
Initial condition for Y_0
y_1 : scalar
Initial condition for Y_1
Ξ± : scalar
Marginal propensity to consume
Ξ² : scalar
Accelerator coefficient
n : int
Number of iterations
Ο : scalar
Volatility parameter. It must be greater than or equal to 0. Set
equal to 0 for a non-stochastic model.
g : scalar
Government spending shock
g_t : int
Time at which government spending shock occurs. Must be specified
when duration != None.
duration : {None, 'permanent', 'one-off'}
Specifies type of government spending shock. If none, government
spending equal to g for all t.
"""
def __init__(self,
y_0=100,
y_1=50,
Ξ±=1.3,
Ξ²=0.2,
Ξ³=10,
n=100,
Ο=0,
g=0,
g_t=0,
duration=None):
def root_type(self):
if all(isinstance(root, complex) for root in self.roots):
return 'Complex conjugate'
elif len(self.roots) > 1:
return 'Double real'
else:
return 'Single real'
13.7. WRAPPING EVERYTHING INTO A CLASS 191
def root_less_than_one(self):
if all(abs(root) < 1 for root in self.roots):
return True
def solution_type(self):
Ο1, Ο2 = self.Ο1, self.Ο2
discriminant = Ο1 ** 2 + 4 * Ο2
if Ο2 >= 1 + Ο1 or Ο2 <= -1:
return 'Explosive oscillations'
elif Ο1 + Ο2 >= 1:
return 'Explosive growth'
elif discriminant < 0:
return 'Damped oscillations'
else:
return 'Steady state'
# Stochastic
else:
οΏ½ = np.random.normal(0, 1, self.n)
return self.Ο1 * x[t - 1] + self.Ο2 * x[t - 2] + self.Ξ³ + g \
+ self.Ο * οΏ½[t]
def generate_series(self):
# No government spending
if self.g == 0:
y_t.append(self._transition(y_t, t))
def summary(self):
print('Summary\n' + '-' * 50)
print(f'Root type: {self.root_type()}')
print(f'Solution type: {self.solution_type()}')
print(f'Roots: {str(self.roots)}')
if self.root_less_than_one() == True:
print('Absolute value of roots is less than one')
else:
print('Absolute value of roots is not less than one')
if self.Ο > 0:
192 CHAPTER 13. OOP III: SAMUELSON MULTIPLIER ACCELERATOR
if self.g != 0:
print('Government spending equal to ' + str(self.g))
if self.duration != None:
print(self.duration.capitalize() +
' government spending shock at t = ' + str(self.g_t))
def plot(self):
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(self.generate_series())
ax.set(xlabel='Iteration', xlim=(0, self.n))
ax.set_ylabel('$Y_t$', rotation=0)
ax.grid()
return fig
def param_plot(self):
fig = param_plot()
ax = fig.gca()
plt.legend(fontsize=12, loc=3)
return fig
Summary
--------------------------------------------------
Root type: Complex conjugate
Solution type: Damped oscillations
Roots: [0.65+0.27838822j 0.65-0.27838822j]
Absolute value of roots is less than one
Stochastic series with Ο = 2
13.7. WRAPPING EVERYTHING INTO A CLASS 193
sam.plot()
[23]: plt.show()
Weβll use our graph to show where the roots lie and how their location is consistent with the
behavior of the path just graphed.
The red + sign shows the location of the roots
sam.param_plot()
[24]: plt.show()
194 CHAPTER 13. OOP III: SAMUELSON MULTIPLIER ACCELERATOR
It turns out that we can use the QuantEcon.py LinearStateSpace class to do much of the
work that we have done from scratch above.
Here is how we map the Samuelson model into an instance of a LinearStateSpace class
"""This script maps the Samuelson model in the the
[25]: ``LinearStateSpace`` class
"""
Ξ± = 0.8
Ξ² = 0.9
Ο1 = Ξ± + Ξ²
Ο2 = -Ξ²
Ξ³ = 10
Ο = 1
g = 10
n = 100
A = [[1, 0, 0],
[Ξ³ + g, Ο1, Ο2],
[0, 1, 0]]
x, y = sam_t.simulate(ts_length=n)
13.8. USING THE LINEARSTATESPACE CLASS 195
axes[-1].set_xlabel('Iteration')
plt.show()
Letβs plot impulse response functions for the instance of the Samuelson model using a
method in the LinearStateSpace class
imres = sam_t.impulse_response()
[26]: imres = np.asarray(imres)
y1 = imres[:, :, 0]
y2 = imres[:, :, 1]
y1.shape
[26]:
(2, 6, 1)
(2, 6, 1)
Now letβs compute the zeros of the characteristic polynomial by simply calculating the eigen-
values of π΄
196 CHAPTER 13. OOP III: SAMUELSON MULTIPLIER ACCELERATOR
A = np.asarray(A)
[27]: w, v = np.linalg.eig(A)
print(w)
We could also create a subclass of LinearStateSpace (inheriting all its methods and at-
tributes) to add more functions to use
class SamuelsonLSS(LinearStateSpace):
[28]:
"""
This subclass creates a Samuelson multiplier-accelerator model
as a linear state space system.
"""
def __init__(self,
y_0=100,
y_1=100,
Ξ±=0.8,
Ξ²=0.9,
Ξ³=10,
Ο=1,
g=10):
self.Ξ±, self.Ξ² = Ξ±, Ξ²
self.y_0, self.y_1, self.g = y_0, y_1, g
self.Ξ³, self.Ο = Ξ³, Ο
self.Ο1 = Ξ± + Ξ²
self.Ο2 = -Ξ²
x, y = self.simulate(ts_length)
axes[-1].set_xlabel('Iteration')
return fig
x, y = self.impulse_response(j)
return fig
13.8.3 Illustrations
samlss.plot_simulation(100, stationary=True)
[31]: plt.show()
samlss.plot_irf(100)
[32]: plt.show()
13.9. PURE MULTIPLIER MODEL 199
samlss.multipliers()
[33]:
array([7.414389, 6.835896, 0.578493])
[33]:
Letβs shut down the accelerator by setting π = 0 to get a pure multiplier model
β’ the absence of cycles gives an idea about why Samuelson included the accelerator
[35]:
200 CHAPTER 13. OOP III: SAMUELSON MULTIPLIER ACCELERATOR
pure_multiplier.plot_irf(100)
[38]:
[38]:
202 CHAPTER 13. OOP III: SAMUELSON MULTIPLIER ACCELERATOR
13.10 Summary
In this lecture, we wrote functions and classes to represent non-stochastic and stochastic ver-
sions of the Samuelson (1939) multiplier-accelerator model, described in [118].
13.10. SUMMARY 203
We saw that different parameter values led to different output paths, which could either be
stationary, explosive, or oscillating.
We also were able to represent the model using the QuantEcon.py LinearStateSpace class.
204 CHAPTER 13. OOP III: SAMUELSON MULTIPLIER ACCELERATOR
Chapter 14
14.1 Contents
β’ Overview 14.2
β’ Iterables and Iterators 14.3
β’ Names and Name Resolution 14.4
β’ Handling Errors 14.5
β’ Decorators and Descriptors 14.6
β’ Generators 14.7
β’ Recursive Function Calls 14.8
β’ Exercises 14.9
β’ Solutions 14.10
14.2 Overview
With this last lecture, our advice is to skip it on first pass, unless you have a burning de-
sire to read it.
Itβs here
A variety of topics are treated in the lecture, including generators, exceptions and descriptors.
205
206 CHAPTER 14. MORE LANGUAGE FEATURES
14.3.1 Iterators
Writing us_cities.txt
f = open('us_cities.txt')
[2]: f.__next__()
We see that file objects do indeed have a __next__ method, and that calling this method
returns the next line in the file.
The next method can also be accessed via the builtin function next(), which directly calls
this method
next(f)
[4]:
'chicago: 2707120\n'
[4]:
(0, 'foo')
[5]:
next(e)
[6]:
(1, 'bar')
[6]:
2009-05-18,9167.05,9167.82,8997.74,9038.69,147800,9038.69
2009-05-15,9150.21,9272.08,9140.90,9265.02,172000,9265.02
2009-05-14,9212.30,9223.77,9052.41,9093.73,169400,9093.73
2009-05-13,9305.79,9379.47,9278.89,9340.49,176000,9340.49
2009-05-12,9358.25,9389.61,9298.61,9298.61,188400,9298.61
2009-05-11,9460.72,9503.91,9342.75,9451.98,230800,9451.98
2009-05-08,9351.40,9464.43,9349.57,9432.83,220200,9432.83
Writing test_table.csv
All iterators can be placed to the right of the in keyword in for loop statements.
In fact this is how the for loop works: If we write
for x in iterator:
<code block>
f = open('somefile.txt', 'r')
for line in f:
# do something
14.3.3 Iterables
You already know that we can put a Python list to the right of in in a for loop
208 CHAPTER 14. MORE LANGUAGE FEATURES
spam
eggs
list
[11]:
next(x)
[12]:
---------------------------------------------------------------------------
<ipython-input-12-92de4e9f6b1e> in <module>
----> 1 next(x)
list
[13]:
y = iter(x)
[14]: type(y)
list_iterator
[14]:
next(y)
[15]:
'foo'
[15]:
next(y)
[16]:
'bar'
[16]:
next(y)
[17]:
---------------------------------------------------------------------------
<ipython-input-17-81b9d2f0f16a> in <module>
----> 1 next(y)
StopIteration:
---------------------------------------------------------------------------
<ipython-input-18-ef50b48e4398> in <module>
----> 1 iter(42)
Some built-in functions that act on sequences also work with iterables
For example
x = [10, -10]
[19]: max(x)
10
[19]:
y = iter(x)
[20]: type(y)
list_iterator
[20]:
max(y)
[21]:
10
[21]:
One thing to remember about iterators is that they are depleted by use
x = [10, -10]
[22]: y = iter(x)
max(y)
10
[22]:
210 CHAPTER 14. MORE LANGUAGE FEATURES
max(y)
[23]:
---------------------------------------------------------------------------
<ipython-input-23-062424e6ec08> in <module>
----> 1 max(y)
We now know that when this statement is executed, Python creates an object of type int in
your computerβs memory, containing
β’ the value 42
β’ some associated attributes
g = f
id(g) == id(f)
True
[25]:
g('test')
[26]:
test
In the first step, a function object is created, and the name f is bound to it.
After binding the name g to the same object, we can use it anywhere we would use f.
What happens when the number of names bound to an object goes to zero?
Hereβs an example of this situation, where the name x is first bound to one object and then
rebound to another
14.4. NAMES AND NAME RESOLUTION 211
x = 'foo'
[27]: id(x)
140341457605440
[27]:
x = 'bar' # No names bound to the first object
[28]:
14.4.2 Namespaces
Writing math2.py
Next letβs import the math module from the standard library
import math
[32]:
These two different bindings of pi exist in different namespaces, each one implemented as a
dictionary.
We can look at the dictionary directly, using module_name.__dict__
212 CHAPTER 14. MORE LANGUAGE FEATURES
import math
[35]:
math.__dict__.items()
import math2
[36]:
math2.__dict__.items()
function max>, 'min': <built-in function min>, 'next': <built-in function next>,
'oct': <built-in function oct>, 'ord': <built-in function ord>, 'pow': <built-in
function pow>, 'print': <built-in function print>, 'repr': <built-in function
repr>, 'round': <built-in function round>, 'setattr': <built-in function
setattr>, 'sorted': <built-in function sorted>, 'sum': <built-in function sum>,
'vars': <built-in function vars>, 'None': None, 'Ellipsis': Ellipsis,
'NotImplemented': NotImplemented, 'False': False, 'True': True, 'bool': <class
'bool'>, 'memoryview': <class 'memoryview'>, 'bytearray': <class 'bytearray'>,
'bytes': <class 'bytes'>, 'classmethod': <class 'classmethod'>, 'complex':
<class 'complex'>, 'dict': <class 'dict'>, 'enumerate': <class 'enumerate'>,
'filter': <class 'filter'>, 'float': <class 'float'>, 'frozenset': <class
'frozenset'>, 'property': <class 'property'>, 'int': <class 'int'>, 'list':
<class 'list'>, 'map': <class 'map'>, 'object': <class 'object'>, 'range':
<class 'range'>, 'reversed': <class 'reversed'>, 'set': <class 'set'>, 'slice':
<class 'slice'>, 'staticmethod': <class 'staticmethod'>, 'str': <class 'str'>,
'super': <class 'super'>, 'tuple': <class 'tuple'>, 'type': <class 'type'>,
'zip': <class 'zip'>, '__debug__': True, 'BaseException': <class
'BaseException'>, 'Exception': <class 'Exception'>, 'TypeError': <class
'TypeError'>, 'StopAsyncIteration': <class 'StopAsyncIteration'>,
'StopIteration': <class 'StopIteration'>, 'GeneratorExit': <class
'GeneratorExit'>, 'SystemExit': <class 'SystemExit'>, 'KeyboardInterrupt':
<class 'KeyboardInterrupt'>, 'ImportError': <class 'ImportError'>,
'ModuleNotFoundError': <class 'ModuleNotFoundError'>, 'OSError': <class
'OSError'>, 'EnvironmentError': <class 'OSError'>, 'IOError': <class 'OSError'>,
'EOFError': <class 'EOFError'>, 'RuntimeError': <class 'RuntimeError'>,
'RecursionError': <class 'RecursionError'>, 'NotImplementedError': <class
'NotImplementedError'>, 'NameError': <class 'NameError'>, 'UnboundLocalError':
<class 'UnboundLocalError'>, 'AttributeError': <class 'AttributeError'>,
'SyntaxError': <class 'SyntaxError'>, 'IndentationError': <class
'IndentationError'>, 'TabError': <class 'TabError'>, 'LookupError': <class
'LookupError'>, 'IndexError': <class 'IndexError'>, 'KeyError': <class
'KeyError'>, 'ValueError': <class 'ValueError'>, 'UnicodeError': <class
'UnicodeError'>, 'UnicodeEncodeError': <class 'UnicodeEncodeError'>,
'UnicodeDecodeError': <class 'UnicodeDecodeError'>, 'UnicodeTranslateError':
<class 'UnicodeTranslateError'>, 'AssertionError': <class 'AssertionError'>,
'ArithmeticError': <class 'ArithmeticError'>, 'FloatingPointError': <class
'FloatingPointError'>, 'OverflowError': <class 'OverflowError'>,
'ZeroDivisionError': <class 'ZeroDivisionError'>, 'SystemError': <class
'SystemError'>, 'ReferenceError': <class 'ReferenceError'>, 'MemoryError':
<class 'MemoryError'>, 'BufferError': <class 'BufferError'>, 'Warning': <class
'Warning'>, 'UserWarning': <class 'UserWarning'>, 'DeprecationWarning': <class
'DeprecationWarning'>, 'PendingDeprecationWarning': <class
'PendingDeprecationWarning'>, 'SyntaxWarning': <class 'SyntaxWarning'>,
'RuntimeWarning': <class 'RuntimeWarning'>, 'FutureWarning': <class
'FutureWarning'>, 'ImportWarning': <class 'ImportWarning'>, 'UnicodeWarning':
<class 'UnicodeWarning'>, 'BytesWarning': <class 'BytesWarning'>,
'ResourceWarning': <class 'ResourceWarning'>, 'ConnectionError': <class
'ConnectionError'>, 'BlockingIOError': <class 'BlockingIOError'>,
'BrokenPipeError': <class 'BrokenPipeError'>, 'ChildProcessError': <class
'ChildProcessError'>, 'ConnectionAbortedError': <class
'ConnectionAbortedError'>, 'ConnectionRefusedError': <class
'ConnectionRefusedError'>, 'ConnectionResetError': <class
'ConnectionResetError'>, 'FileExistsError': <class 'FileExistsError'>,
'FileNotFoundError': <class 'FileNotFoundError'>, 'IsADirectoryError': <class
'IsADirectoryError'>, 'NotADirectoryError': <class 'NotADirectoryError'>,
'InterruptedError': <class 'InterruptedError'>, 'PermissionError': <class
'PermissionError'>, 'ProcessLookupError': <class 'ProcessLookupError'>,
'TimeoutError': <class 'TimeoutError'>, 'open': <built-in function open>,
'copyright': Copyright (c) 2001-2019 Python Software Foundation.
All Rights Reserved.
As you know, we access elements of the namespace using the dotted attribute notation
math.pi
[37]:
3.141592653589793
[37]:
['__doc__',
[40]: '__file__',
'__loader__',
'__name__',
'__package__',
'__spec__',
'acos',
'acosh',
'asin',
'asinh']
print(math.__doc__)
[41]:
math.__name__
[42]:
'math'
[42]:
__main__
When we run a script using IPythonβs run command, the contents of the file are executed as
part of __main__ too.
To see this, letβs create a file mod.py that prints its own __name__ attribute
%%file mod.py
[44]: print(__name__)
Writing mod.py
mod
216 CHAPTER 14. MORE LANGUAGE FEATURES
__main__
In the second case, the code is executed as part of __main__, so __name__ is equal to
__main__.
To see the contents of the namespace of __main__ we use vars() rather than
vars(__main__) .
If you do this in IPython, you will see a whole lot of variables that IPython needs, and has
initialized when you started up your session.
If you prefer to see only the variables you have initialized, use whos
x = 2
[47]: y = 3
import numpy as np
%whos
import amodule
At this point, the interpreter creates a namespace for the module amodule and starts exe-
cuting commands in the module.
While this occurs, the namespace amodule.__dict__ is the global namespace.
14.4. NAMES AND NAME RESOLUTION 217
Once execution of the module finishes, the interpreter returns to the module from where the
import statement was made.
In this case itβs __main__, so the namespace of __main__ again becomes the global names-
pace.
Important fact: When we call a function, the interpreter creates a local namespace for that
function, and registers the variables in that namespace.
The reason for this will be explained in just a moment.
Variables in the local namespace are called local variables.
After the function returns, the namespace is deallocated and lost.
While the function is executing, we can view the contents of the local namespace with
locals().
For example, consider
def f(x):
[48]: a = 2
print(locals())
return a * x
{'x': 1, 'a': 2}
2
[49]:
We have been using various built-in functions, such as max(), dir(), str(), list(),
len(), range(), type(), etc.
How does access to these names work?
dir()[0:10]
[50]:
['In', 'Out', '_', '_11', '_13', '_14', '_15', '_16', '_19', '_2']
[50]:
dir(__builtins__)[0:10]
[51]:
['ArithmeticError',
[51]: 'AssertionError',
'AttributeError',
'BaseException',
'BlockingIOError',
218 CHAPTER 14. MORE LANGUAGE FEATURES
'BrokenPipeError',
'BufferError',
'BytesWarning',
'ChildProcessError',
'ConnectionAbortedError']
But __builtins__ is special, because we can always access them directly as well
max
[53]:
<function max>
[53]:
__builtins__.max == max
[54]:
True
[54]:
If the interpreter is executing a function, then the directly accessible namespaces are
Here f is the enclosing function for g, and each function gets its own namespaces.
14.4. NAMES AND NAME RESOLUTION 219
Now we can give the rule for how namespace resolution works:
The order in which the interpreter searches for names is
If the name is not in any of these namespaces, the interpreter raises a NameError.
This is called the LEGB rule (local, enclosing, global, builtin).
Hereβs an example that helps to illustrate .
Consider a script test.py that looks as follows
%%file test.py
[56]: def g(x):
a = 1
x = x + a
return x
a = 0
y = g(10)
print("a = ", a, "y = ", y)
Writing test.py
a = 0 y = 11
x
[58]:
2
[58]:
First,
This is a good time to say a little more about mutable vs immutable objects.
Consider the code segment
def f(x):
[59]: x = x + 1
return x
x = 1
print(f(x), x)
2 1
We now understand what will happen here: The code prints 2 as the value of f(x) and 1 as
the value of x.
First f and x are registered in the global namespace.
The call f(x) creates a local namespace and adds x to it, bound to 1.
Next, this local x is rebound to the new integer object 2, and this value is returned.
None of this affects the global x.
However, itβs a different story when we use a mutable data type such as a list
def f(x):
[]: x[0] = x[0] + 1
return x
x = [1]
print(f(x), x)
[2] [2]
π
1
π 2 βΆ= β(π¦π β π¦)Μ 2 π¦ Μ = sample mean
π β 1 π=1
β’ Because the debugging information provided by the interpreter is often less useful than
the information on possible errors you have in your head when writing code.
β’ Because errors causing execution to stop are frustrating if youβre in the middle of a
large computation.
β’ Because itβs reduces confidence in your code on the part of your users (if you are writing
for others).
14.5.1 Assertions
If we run this with an array of length one, the program will terminate and print our error
message
var([1])
[62]:
---------------------------------------------------------------------------
<ipython-input-62-8419b6ab38ec> in <module>
----> 1 var([1])
<ipython-input-61-e6ffb16a7098> in var(y)
1 def var(y):
2 n = len(y)
----> 3 assert n > 1, 'Sample size must be greater than one.'
4 return np.sum((y - y.mean())**2) / float(n-1)
The approach used above is a bit limited, because it always leads to termination.
Sometimes we can handle errors more gracefully, by treating special cases.
Letβs look at how this is done.
Exceptions
Hereβs an example of a common error type
def f:
[63]:
Since illegal syntax cannot be executed, a syntax error terminates execution of the program.
Hereβs a different kind of error, unrelated to syntax
1 / 0
[64]:
---------------------------------------------------------------------------
<ipython-input-64-bc757c3fda29> in <module>
----> 1 1 / 0
Hereβs another
x1 = y1
[65]:
---------------------------------------------------------------------------
<ipython-input-65-a7b8d65e9e45> in <module>
----> 1 x1 = y1
And another
14.5. HANDLING ERRORS 223
'foo' + 6
[66]:
---------------------------------------------------------------------------
<ipython-input-66-216809d6e6fe> in <module>
----> 1 'foo' + 6
And another
X = []
[67]: x = X[0]
---------------------------------------------------------------------------
<ipython-input-67-082a18d7a0aa> in <module>
1 X = []
----> 2 x = X[0]
f(0.0)
[71]:
f('foo')
[75]:
f('foo')
[79]:
def f(x):
[80]: try:
return 1.0 / x
except:
print('Error. Returned None')
return None
Letβs look at some special syntax elements that are routinely used by Python developers.
You might not need the following concepts immediately, but you will see them in other peo-
pleβs code.
Hence you need to understand them at some stage of your Python education.
14.6.1 Decorators
Decorators are a bit of syntactic sugar that, while easily avoided, have turned out to be popu-
lar.
Itβs very easy to say what decorators do.
On the other hand it takes a bit of effort to explain why you might use them.
An Example
Suppose we are working on a program that looks something like this
import numpy as np
[81]:
def f(x):
return np.log(np.log(x))
def g(x):
return np.sqrt(42 * x)
Now suppose thereβs a problem: occasionally negative numbers get fed to f and g in the cal-
culations that follow.
If you try it, youβll see that when these functions are called with negative numbers they re-
turn a NumPy object called nan .
This stands for βnot a numberβ (and indicates that you are trying to evaluate a mathematical
function at a point where it is not defined).
Perhaps this isnβt what we want, because it causes other problems that are hard to pick up
later on.
Suppose that instead we want the program to terminate whenever this happens, with a sensi-
ble error message.
This change is easy enough to implement
import numpy as np
[82]:
def f(x):
226 CHAPTER 14. MORE LANGUAGE FEATURES
def g(x):
assert x >= 0, "Argument must be nonnegative"
return np.sqrt(42 * x)
Notice however that there is some repetition here, in the form of two identical lines of code.
Repetition makes our code longer and harder to maintain, and hence is something we try
hard to avoid.
Here itβs not a big deal, but imagine now that instead of just f and g, we have 20 such func-
tions that we need to modify in exactly the same way.
This means we need to repeat the test logic (i.e., the assert line testing nonnegativity) 20
times.
The situation is still worse if the test logic is longer and more complicated.
In this kind of scenario the following approach would be neater
import numpy as np
[83]:
def check_nonneg(func):
def safe_function(x):
assert x >= 0, "Argument must be nonnegative"
return func(x)
return safe_function
def f(x):
return np.log(np.log(x))
def g(x):
return np.sqrt(42 * x)
f = check_nonneg(f)
g = check_nonneg(g)
# Program continues with various calculations using f and g
def g(x):
return np.sqrt(42 * x)
f = check_nonneg(f)
g = check_nonneg(g)
with
@check_nonneg
[86]: def f(x):
return np.log(np.log(x))
@check_nonneg
def g(x):
return np.sqrt(42 * x)
14.6.2 Descriptors
One potential problem we might have here is that a user alters one of these variables but not
the other
car = Car()
[88]: car.miles
1000
[88]:
car.kms
[89]:
228 CHAPTER 14. MORE LANGUAGE FEATURES
1610.0
[89]:
car.miles = 6000
[90]: car.kms
1610.0
[90]:
In the last two lines we see that miles and kms are out of sync.
What we really want is some mechanism whereby each time a user sets one of these variables,
the other is automatically updated.
A Solution
In Python, this issue is solved using descriptors.
A descriptor is just a Python object that implements certain methods.
These methods are triggered when the object is accessed through dotted attribute notation.
The best way to understand this is to see it in action.
Consider this alternative version of the Car class
class Car:
[91]:
def __init__(self, miles=1000):
self._miles = miles
self._kms = miles * 1.61
def get_miles(self):
return self._miles
def get_kms(self):
return self._kms
1000
[92]:
car.miles = 6000
[93]: car.kms
9660.0
[93]:
The methods get_miles, set_miles, get_kms and set_kms define what happens when
you get (i.e. access) or set (bind) these variables
The builtin Python function property takes getter and setter methods and creates a prop-
erty.
For example, after car is created as an instance of Car, the object car.miles is a property.
Being a property, when we set its value via car.miles = 6000 its setter method is trig-
gered β in this case set_miles.
Decorators and Properties
These days its very common to see the property function used via a decorator.
Hereβs another version of our Car class that works as before but now uses decorators to set
up the properties
class Car:
[94]:
def __init__(self, miles=1000):
self._miles = miles
self._kms = miles * 1.61
@property
def miles(self):
return self._miles
@property
def kms(self):
return self._kms
@miles.setter
def miles(self, value):
self._miles = value
self._kms = value * 1.61
@kms.setter
def kms(self, value):
self._kms = value
self._miles = value / 1.61
14.7 Generators
tuple
[95]:
plural = [string + 's' for string in singular]
[96]: plural
generator
[98]:
next(plural)
[99]:
'dogs'
[99]:
next(plural)
[100]:
'cats'
[100]:
next(plural)
[101]:
'birds'
[101]:
The function sum() calls next() to get the items, adds successive terms.
In fact, we can omit the outer brackets in this case
sum(x * x for x in range(10))
[103]:
285
[103]:
The most flexible way to create generator objects is to use generator functions.
Letβs look at some examples.
Example 1
Hereβs a very simple example of a generator function
def f():
[104]: yield 'start'
yield 'middle'
yield 'end'
14.7. GENERATORS 231
It looks like a function, but uses a keyword yield that we havenβt met before.
Letβs see how it works after running this code
type(f)
[105]:
function
[105]:
gen = f()
[106]: gen
---------------------------------------------------------------------------
<ipython-input-110-6e72e47198db> in <module>
----> 1 next(gen)
StopIteration:
The generator function f() is used to create generator objects (in this case gen).
Generators are iterators, because they support a next method.
The first call to next(gen)
The second call to next(gen) starts executing from the next line
def f():
[111]: yield 'start'
yield 'middle' # This line!
yield 'end'
def g(x):
[112]: while x < 100:
yield x
x = x * x
generator
[114]:
next(gen)
[115]:
2
[115]:
next(gen)
[116]:
4
[116]:
next(gen)
[117]:
16
[117]:
next(gen)
[118]:
---------------------------------------------------------------------------
<ipython-input-118-6e72e47198db> in <module>
----> 1 next(gen)
StopIteration:
β’ The body of g() executes until the line yield x, and the value of x is returned.
def g(x):
[120]: while 1:
yield x
x = x * x
5000779
[121]:
But we are creating two huge lists here, range(n) and draws.
This uses lots of memory and is very slow.
If we make n even bigger then this happens
n = 100000000
[122]: draws = [random.uniform(0, 1) < 0.5 for i in range(n)]
In summary, iterables
This is not something that you will use every day, but it is still useful β you should learn it
at some stage.
234 CHAPTER 14. MORE LANGUAGE FEATURES
What happens here is that each successive call uses itβs own frame in the stack
β’ a frame is where the local variables of a given function call are held
β’ stack is memory used to process function calls
β a First In Last Out (FILO) queue
This example is somewhat contrived, since the first (iterative) solution would usually be pre-
ferred to the recursive solution.
Weβll meet less contrived applications of recursion later on.
14.9 Exercises
14.9.1 Exercise 1
The first few numbers in the sequence are 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55.
Write a function to recursively compute the π‘-th Fibonacci number for any π‘.
14.9.2 Exercise 2
Complete the following code, and test it using this csv file, which we assume that youβve put
in your current working directory
dates = column_iterator('test_table.csv', 1)
14.9.3 Exercise 3
prices
3
8
7
21
Using try β except, write a program to read in the contents of the file and sum the num-
bers, ignoring lines without numbers.
14.10 Solutions
14.10.1 Exercise 1
Letβs test it
print([x(i) for i in range(10)])
[129]:
14.10.2 Exercise 2
dates = column_iterator('test_table.csv', 1)
i = 1
for date in dates:
print(date)
if i == 10:
break
i += 1
Date
2009-05-21
2009-05-20
2009-05-19
2009-05-18
2009-05-15
2009-05-14
2009-05-13
2009-05-12
2009-05-11
14.10.3 Exercise 3
7
21
Writing numbers.txt
f = open('numbers.txt')
[132]:
total = 0.0
for line in f:
try:
total += float(line)
except ValueError:
pass
f.close()
print(total)
39.0
Chapter 15
Debugging
15.1 Contents
β’ Overview 15.2
β’ Debugging 15.3
β’ Other Useful Magics 15.4
βDebugging is twice as hard as writing the code in the first place. Therefore, if
you write the code as cleverly as possible, you are, by definition, not smart enough
to debug it.β β Brian Kernighan
15.2 Overview
Are you one of those programmers who fills their code with print statements when trying to
debug their programs?
Hey, we all used to do that.
(OK, sometimes we still do thatβ¦)
But once you start writing larger programs youβll need a better system.
Debugging tools for Python vary across platforms, IDEs and editors.
Here weβll focus on Jupyter and leave you to explore other settings.
Weβll need the following imports
import numpy as np
[1]: import matplotlib.pyplot as plt
%matplotlib inline
15.3 Debugging
237
238 CHAPTER 15. DEBUGGING
def plot_log():
[2]: fig, ax = plt.subplots(2, 1)
x = np.linspace(1, 2, 10)
ax.plot(x, np.log(x))
plt.show()
---------------------------------------------------------------------------
<ipython-input-2-c32a2280f47b> in <module>
5 plt.show()
6
----> 7 plot_log() # Call the function, generate plot
<ipython-input-2-c32a2280f47b> in plot_log()
2 fig, ax = plt.subplots(2, 1)
3 x = np.linspace(1, 2, 10)
----> 4 ax.plot(x, np.log(x))
5 plt.show()
6
This code is intended to plot the log function over the interval [1, 2].
But thereβs an error here: plt.subplots(2, 1) should be just plt.subplots().
(The call plt.subplots(2, 1) returns a NumPy array containing two axes objects, suit-
able for having two subplots on the same figure)
The traceback shows that the error occurs at the method call ax.plot(x, np.log(x)).
15.3. DEBUGGING 239
The error occurs because we have mistakenly made ax a NumPy array, and a NumPy array
has no plot method.
But letβs pretend that we donβt understand this for the moment.
We might suspect thereβs something wrong with ax but when we try to investigate this ob-
ject, we get the following exception:
ax
[3]:
---------------------------------------------------------------------------
<ipython-input-3-b00e77935981> in <module>
----> 1 ax
The problem is that ax was defined inside plot_log(), and the name is lost once that func-
tion terminates.
Letβs try doing it a different way.
We run the first cell block again, generating the same error
def plot_log():
[4]: fig, ax = plt.subplots(2, 1)
x = np.linspace(1, 2, 10)
ax.plot(x, np.log(x))
plt.show()
---------------------------------------------------------------------------
<ipython-input-4-c32a2280f47b> in <module>
5 plt.show()
6
----> 7 plot_log() # Call the function, generate plot
<ipython-input-4-c32a2280f47b> in plot_log()
2 fig, ax = plt.subplots(2, 1)
3 x = np.linspace(1, 2, 10)
----> 4 ax.plot(x, np.log(x))
5 plt.show()
6
%debug
You should be dropped into a new prompt that looks something like this
ipdb>
ipdb> ax
array([<matplotlib.axes.AxesSubplot object at 0x290f5d0>,
<matplotlib.axes.AxesSubplot object at 0x2930810>], dtype=object)
Itβs now very clear that ax is an array, which clarifies the source of the problem.
To find out what else you can do from inside ipdb (or pdb), use the online help
ipdb> h
Undocumented commands:
======================
retval rv
ipdb> h c
c(ont(inue))
Continue execution, only stop when a breakpoint is encountered.
plot_log()
Here the original problem is fixed, but weβve accidentally written np.logspace(1, 2,
10) instead of np.linspace(1, 2, 10).
242 CHAPTER 15. DEBUGGING
Now there wonβt be any exception, but the plot wonβt look right.
To investigate, it would be helpful if we could inspect variables like x during execution of the
function.
To this end, we add a βbreak pointβ by inserting breakpoint() inside the function code
block
def plot_log():
breakpoint()
fig, ax = plt.subplots()
x = np.logspace(1, 2, 10)
ax.plot(x, np.log(x))
plt.show()
plot_log()
Now letβs run the script, and investigate via the debugger
> <ipython-input-6-a188074383b7>(6)plot_log()
-> fig, ax = plt.subplots()
(Pdb) n
> <ipython-input-6-a188074383b7>(7)plot_log()
-> x = np.logspace(1, 2, 10)
(Pdb) n
> <ipython-input-6-a188074383b7>(8)plot_log()
-> ax.plot(x, np.log(x))
(Pdb) x
array([ 10. , 12.91549665, 16.68100537, 21.5443469 ,
27.82559402, 35.93813664, 46.41588834, 59.94842503,
77.42636827, 100. ])
We used n twice to step forward through the code (one line at a time).
Then we printed the value of x to see what was happening with that variable.
To exit from the debugger, use q.
243
Chapter 16
Pandas
16.1 Contents
β’ Overview 16.2
β’ Series 16.3
β’ DataFrames 16.4
β’ Exercises 16.6
β’ Solutions 16.7
In addition to whatβs in Anaconda, this lecture will need the following libraries:
!pip install --upgrade pandas-datareader
[1]:
245
246 CHAPTER 16. PANDAS
/home/ubuntu/anaconda3/lib/python3.7/site-packages (from
requests>=2.3.0->pandas-datareader) (3.0.4)
Requirement already satisfied, skipping upgrade: certifi>=2017.4.17 in
/home/ubuntu/anaconda3/lib/python3.7/site-packages (from
requests>=2.3.0->pandas-datareader) (2019.6.16)
Requirement already satisfied, skipping upgrade: six>=1.5 in
/home/ubuntu/anaconda3/lib/python3.7/site-packages (from python-
dateutil>=2.5.0->pandas>=0.21->pandas-datareader) (1.12.0)
16.2 Overview
Just as NumPy provides the basic array data type plus core array operations, pandas
β’ reading in data
β’ adjusting indices
β’ working with dates and time series
β’ sorting, grouping, re-ordering and general data munging 1
β’ dealing with missing values, etc., etc.
More sophisticated statistical functionality is left to other packages, such as statsmodels and
scikit-learn, which are built on top of pandas.
This lecture will provide a basic introduction to pandas.
Throughout the lecture, we will assume that the following imports have taken place
16.3. SERIES 247
import pandas as pd
[2]: import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import requests
16.3 Series
Two important data types defined by pandas are Series and DataFrame.
You can think of a Series as a βcolumnβ of data, such as a collection of observations on a
single variable.
A DataFrame is an object for storing related columns of data.
Letβs start with Series
s = pd.Series(np.random.randn(4), name='daily returns')
[3]: s
0 1.520681
[3]: 1 0.619485
2 -0.127222
3 0.163491
Name: daily returns, dtype: float64
Here you can imagine the indices 0, 1, 2, 3 as indexing four listed companies, and the
values being daily returns on their shares.
Pandas Series are built on top of NumPy arrays and support many similar operations
s * 100
[4]:
0 152.068086
[4]: 1 61.948511
2 -12.722188
3 16.349095
Name: daily returns, dtype: float64
np.abs(s)
[5]:
0 1.520681
[5]: 1 0.619485
2 0.127222
3 0.163491
Name: daily returns, dtype: float64
AMZN 1.520681
[7]: AAPL 0.619485
MSFT -0.127222
GOOG 0.163491
Name: daily returns, dtype: float64
Viewed in this way, Series are like fast, efficient Python dictionaries (with the restriction
that the items in the dictionary all have the same typeβin this case, floats).
In fact, you can use much of the same syntax as Python dictionaries
s['AMZN']
[8]:
1.520680864827497
[8]:
s['AMZN'] = 0
[9]: s
AMZN 0.000000
[9]: AAPL 0.619485
MSFT -0.127222
GOOG 0.163491
Name: daily returns, dtype: float64
'AAPL' in s
[10]:
True
[10]:
16.4 DataFrames
While a Series is a single column of data, a DataFrame is several columns, one for each
variable.
In essence, a DataFrame in pandas is analogous to a (highly optimized) Excel spreadsheet.
Thus, it is a powerful tool for representing and analyzing data that are naturally organized
into rows and columns, often with descriptive indexes for individual rows and individual
columns.
Letβs look at an example that reads data from the CSV file pandas/data/test_pwt.csv
that can be downloaded here.
Hereβs the content of test_pwt.csv
"country","country isocode","year","POP","XRAT","tcgdp","cc","cg"
"Argentina","ARG","2000","37335.653","0.9995","295072.21869","75.716805379","5.5
"Australia","AUS","2000","19053.186","1.72483","541804.6521","67.759025993","6.7
"India","IND","2000","1006300.297","44.9416","1728144.3748","64.575551328","14.0
"Israel","ISR","2000","6114.57","4.07733","129253.89423","64.436450847","10.2666
"Malawi","MWI","2000","11801.505","59.543808333","5026.2217836","74.707624181","
"South Africa","ZAF","2000","45064.098","6.93983","227242.36949","72.718710427",
"United States","USA","2000","282171.957","1","9898700","72.347054303","6.032453
"Uruguay","URY","2000","3219.793","12.099591667","25255.961693","78.978740282","
Supposing you have this data saved as test_pwt.csv in the present working directory (type
%pwd in Jupyter to see what this is), it can be read in as follows:
16.4. DATAFRAMES 249
df = pd.read_csv('https://github.com/QuantEcon/QuantEcon.lectures.code/raw/master/pandas/
[11]: βͺdata/test_pwt.csv')
type(df)
pandas.core.frame.DataFrame
[11]:
df
[12]:
country country isocode year POP XRAT tcgdp \
[12]: 0 Argentina ARG 2000 37335.653 0.999500 2.950722e+05
1 Australia AUS 2000 19053.186 1.724830 5.418047e+05
2 India IND 2000 1006300.297 44.941600 1.728144e+06
3 Israel ISR 2000 6114.570 4.077330 1.292539e+05
4 Malawi MWI 2000 11801.505 59.543808 5.026222e+03
5 South Africa ZAF 2000 45064.098 6.939830 2.272424e+05
6 United States USA 2000 282171.957 1.000000 9.898700e+06
7 Uruguay URY 2000 3219.793 12.099592 2.525596e+04
cc cg
0 75.716805 5.578804
1 67.759026 6.720098
2 64.575551 14.072206
3 64.436451 10.266688
4 74.707624 11.658954
5 72.718710 5.726546
6 72.347054 6.032454
7 78.978740 5.108068
We can select particular rows using standard Python array slicing notation
df[2:5]
[13]:
country country isocode year POP XRAT tcgdp \
[13]: 2 India IND 2000 1006300.297 44.941600 1.728144e+06
3 Israel ISR 2000 6114.570 4.077330 1.292539e+05
4 Malawi MWI 2000 11801.505 59.543808 5.026222e+03
cc cg
2 64.575551 14.072206
3 64.436451 10.266688
4 74.707624 11.658954
To select columns, we can pass a list containing the names of the desired columns represented
as strings
df[['country', 'tcgdp']]
[14]:
country tcgdp
[14]: 0 Argentina 2.950722e+05
1 Australia 5.418047e+05
2 India 1.728144e+06
3 Israel 1.292539e+05
4 Malawi 5.026222e+03
5 South Africa 2.272424e+05
6 United States 9.898700e+06
7 Uruguay 2.525596e+04
To select both rows and columns using integers, the iloc attribute should be used with the
format .iloc[rows, columns]
df.iloc[2:5, 0:4]
[15]:
country country isocode year POP
[15]: 2 India IND 2000 1006300.297
3 Israel ISR 2000 6114.570
4 Malawi MWI 2000 11801.505
To select rows and columns using a mixture of integers and labels, the loc attribute can be
250 CHAPTER 16. PANDAS
Letβs imagine that weβre only interested in population and total GDP (tcgdp).
One way to strip the data frame df down to only these variables is to overwrite the
dataframe using the selection method described above
df = df[['country', 'POP', 'tcgdp']]
[17]: df
Here the index 0, 1,..., 7 is redundant because we can use the country names as an in-
dex.
To do this, we set the index to be the country variable in the dataframe
df = df.set_index('country')
[18]: df
POP tcgdp
[18]: country
Argentina 37335.653 2.950722e+05
Australia 19053.186 5.418047e+05
India 1006300.297 1.728144e+06
Israel 6114.570 1.292539e+05
Malawi 11801.505 5.026222e+03
South Africa 45064.098 2.272424e+05
United States 282171.957 9.898700e+06
Uruguay 3219.793 2.525596e+04
Next, weβre going to add a column showing real GDP per capita, multiplying by 1,000,000 as
we go because total GDP is in millions
df['GDP percap'] = df['total GDP'] * 1e6 / df['population']
[21]: df
One of the nice things about pandas DataFrame and Series objects is that they have
methods for plotting and visualization that work through Matplotlib.
For example, we can easily generate a bar plot of GDP per capita
df['GDP percap'].plot(kind='bar')
[22]: plt.show()
252 CHAPTER 16. PANDAS
At the moment the data frame is ordered alphabetically on the countriesβletβs change it to
GDP per capita
df = df.sort_values(by='GDP percap', ascending=False)
[23]: df
https://research.stlouisfed.org/fred2/series/UNRATE/downloaddata/UNRATE.csv
One option is to use requests, a standard Python library for requesting data over the Inter-
net.
To begin, try the following code on your computer
r = requests.get('http://research.stlouisfed.org/fred2/series/UNRATE/downloaddata/UNRATE.
[25]: βͺcsv')
1. You are not connected to the Internet β hopefully, this isnβt the case.
2. Your machine is accessing the Internet through a proxy server, and Python isnβt aware
of this.
'DATE,VALUE\r'
[26]:
source[1]
[27]:
'1948-01-01,3.4\r'
[27]:
source[2]
[28]:
254 CHAPTER 16. PANDAS
'1948-02-01,3.8\r'
[28]:
We could now write some additional code to parse this text and store it as an array.
But this is unnecessary β pandasβ read_csv function can handle the task for us.
We use parse_dates=True so that pandas recognizes our dates column, allowing for simple
date filtering
data = pd.read_csv(url, index_col=0, parse_dates=True)
[29]:
The data has been read into a pandas DataFrame called data that we can now manipulate in
the usual way
type(data)
[30]:
pandas.core.frame.DataFrame
[30]:
data.head() # A useful method to get a quick look at a data frame
[31]:
VALUE
[31]: DATE
1948-01-01 3.4
1948-02-01 3.8
1948-03-01 4.0
1948-04-01 3.9
1948-05-01 3.5
pd.set_option('precision', 1)
[32]: data.describe() # Your output might differ slightly
VALUE
[32]: count 861.0
mean 5.7
std 1.6
min 2.5
25% 4.5
50% 5.6
75% 6.8
max 10.8
We can also plot the unemployment rate from 2006 to 2012 as follows
data['2006':'2012'].plot()
[33]: plt.show()
16.5. ON-LINE DATA SOURCES 255
The maker of pandas has also authored a library called pandas_datareader that gives pro-
grammatic access to many data sources straight from the Jupyter notebook.
While some sources require an access key, many of the most important (e.g., FRED, OECD,
EUROSTAT and the World Bank) are free to use.
For now letβs work through one example of downloading and plotting data β this time from
the World Bank.
The World Bank collects and organizes data on a huge range of indicators.
For example, hereβs some data on government debt as a ratio to GDP.
The next code example fetches the data for you and plots time series for the US and Aus-
tralia
from pandas_datareader import wb
[34]:
govt_debt = wb.download(indicator='GC.DOD.TOTL.GD.ZS', country=['US', 'AU'], start=2005,οΏ½
βͺend=2016).stack().unstack(0)
ind = govt_debt.index.droplevel(-1)
govt_debt.index = ind
ax = govt_debt.plot(lw=2)
plt.title("Government Debt to GDP (%)")
plt.show()
256 CHAPTER 16. PANDAS
The documentation provides more details on how to access various data sources.
16.6 Exercises
16.6.1 Exercise 1
Write a program to calculate the percentage price change over 2013 for the following shares
ticker_list = {'INTC': 'Intel',
[35]: 'MSFT': 'Microsoft',
'IBM': 'IBM',
'BHP': 'BHP',
'TM': 'Toyota',
'AAPL': 'Apple',
'AMZN': 'Amazon',
'BA': 'Boeing',
'QCOM': 'Qualcomm',
'KO': 'Coca-Cola',
'GOOG': 'Google',
'SNE': 'Sony',
'PTR': 'PetroChina'}
A dataset of daily closing prices for the above firms can be found in
pandas/data/ticker_data.csv and can be downloaded here.
Plot the result as a bar graph like this one
16.7. SOLUTIONS 257
16.7 Solutions
16.7.1 Exercise 1
ticker = pd.read_csv('https://github.com/QuantEcon/QuantEcon.lectures.code/raw/master/
[36]: βͺpandas/data/ticker_data.csv')
ticker.set_index('Date', inplace=True)
price_change = pd.Series()
price_change.sort_values(inplace=True)
fig, ax = plt.subplots(figsize=(10,8))
price_change.plot(kind='bar', ax=ax)
plt.show()
258 CHAPTER 16. PANDAS
Footnotes
[1] Wikipedia defines munging as cleaning data from one raw form into a structured, purged
one.
Chapter 17
17.1 Contents
β’ Overview 17.2
β’ Exercises 17.7
β’ Solutions 17.8
17.2 Overview
pandas (derived from βpanelβ and βdataβ) contains powerful and easy-to-use tools for solving
exactly these kinds of problems.
In what follows, we will use a panel data set of real minimum wages from the OECD to cre-
ate:
259
260 CHAPTER 17. PANDAS FOR PANEL DATA
We will begin by reading in our long format panel data from a CSV file and reshaping the
resulting DataFrame with pivot_table to build a MultiIndex.
Additional detail will be added to our DataFrame using pandasβ merge function, and data
will be summarized with the groupby function.
Most of this lecture was created by Natasha Watkins.
We will read in a dataset from the OECD of real minimum wages in 32 countries and assign
it to realwage.
The dataset pandas_panel/realwage.csv can be downloaded here.
Make sure the file is in your current working directory
import pandas as pd
[1]:
# Display 6 columns for viewing purposes
pd.set_option('display.max_columns', 6)
realwage = pd.read_csv('https://github.com/QuantEcon/QuantEcon.lectures.code/raw/master/
βͺpandas_panel/realwage.csv')
The data is currently in long format, which is difficult to analyze when there are several di-
mensions to the data.
We will use pivot_table to create a wide format panel, with a MultiIndex to handle
higher dimensional data.
pivot_table arguments should specify the data (values), the index, and the columns we
want in our resulting dataframe.
By passing a list in columns, we can create a MultiIndex in our column axis
realwage = realwage.pivot_table(values='value',
[3]: index='Time',
columns=['Country', 'Series', 'Pay period'])
realwage.head()
17.3. SLICING AND RESHAPING DATA 261
Country Australia \
[3]: Series In 2015 constant prices at 2015 USD PPPs
Pay period Annual Hourly
Time
2006-01-01 20,410.65 10.33
2007-01-01 21,087.57 10.67
2008-01-01 20,718.24 10.48
2009-01-01 20,984.77 10.62
2010-01-01 20,879.33 10.57
Country β¦ \
Series In 2015 constant prices at 2015 USD exchange rates β¦
Pay period Annual β¦
Time β¦
2006-01-01 23,826.64 β¦
2007-01-01 24,616.84 β¦
2008-01-01 24,185.70 β¦
2009-01-01 24,496.84 β¦
2010-01-01 24,373.76 β¦
Country
Series In 2015 constant prices at 2015 USD exchange rates
Pay period Annual Hourly
Time
2006-01-01 12,594.40 6.05
2007-01-01 12,974.40 6.24
2008-01-01 14,097.56 6.78
2009-01-01 15,756.42 7.58
2010-01-01 16,391.31 7.88
To more easily filter our time series data, later on, we will convert the index into a
DateTimeIndex
realwage.index = pd.to_datetime(realwage.index)
[4]: type(realwage.index)
pandas.core.indexes.datetimes.DatetimeIndex
[4]:
The columns contain multiple levels of indexing, known as a MultiIndex, with levels being
ordered hierarchically (Country > Series > Pay period).
A MultiIndex is the simplest and most flexible way to manage panel data in pandas
type(realwage.columns)
[5]:
pandas.core.indexes.multi.MultiIndex
[5]:
realwage.columns.names
[6]:
FrozenList(['Country', 'Series', 'Pay period'])
[6]:
Like before, we can select the country (the top level of our MultiIndex)
realwage['United States'].head()
[7]:
262 CHAPTER 17. PANDAS FOR PANEL DATA
Stacking and unstacking levels of the MultiIndex will be used throughout this lecture to
reshape our dataframe into a format we need.
.stack() rotates the lowest level of the column MultiIndex to the row index
(.unstack() works in the opposite direction - try it out)
realwage.stack().head()
[8]:
Country Australia \
[8]: Series In 2015 constant prices at 2015 USD PPPs
Time Pay period
2006-01-01 Annual 20,410.65
Hourly 10.33
2007-01-01 Annual 21,087.57
Hourly 10.67
2008-01-01 Annual 20,718.24
Country \
Series In 2015 constant prices at 2015 USD exchange rates
Time Pay period
2006-01-01 Annual 23,826.64
Hourly 12.06
2007-01-01 Annual 24,616.84
Hourly 12.46
2008-01-01 Annual 24,185.70
Country Belgium β¦ \
Series In 2015 constant prices at 2015 USD PPPs β¦
Time Pay period β¦
2006-01-01 Annual 21,042.28 β¦
Hourly 10.09 β¦
2007-01-01 Annual 21,310.05 β¦
Hourly 10.22 β¦
2008-01-01 Annual 21,416.96 β¦
Country
Series In 2015 constant prices at 2015 USD exchange rates
Time Pay period
2006-01-01 Annual 12,594.40
Hourly 6.05
2007-01-01 Annual 12,974.40
Hourly 6.24
2008-01-01 Annual 14,097.56
[5 rows x 64 columns]
We can also pass in an argument to select the level we would like to stack
realwage.stack(level='Country').head()
[9]:
Series In 2015 constant prices at 2015 USD PPPs \
[9]: Pay period Annual Hourly
Time Country
2006-01-01 Australia 20,410.65 10.33
Belgium 21,042.28 10.09
Brazil 3,310.51 1.41
Canada 13,649.69 6.56
Chile 5,201.65 2.22
Time
Series In 2015 constant prices at 2015 USD exchange rates
Pay period Annual Hourly
Country
Australia 25,349.90 12.83
Belgium 20,753.48 9.95
Brazil 2,842.28 1.21
Canada 17,367.24 8.35
Chile 4,251.49 1.81
For the rest of lecture, we will work with a dataframe of the hourly real minimum wages
across countries and time, measured in 2015 US dollars.
To create our filtered dataframe (realwage_f), we can use the xs method to select values
at lower levels in the multiindex, while keeping the higher levels (countries in this case)
realwage_f = realwage.xs(('Hourly', 'In 2015 constant prices at 2015 USD exchange rates'),
[11]: level=('Pay period', 'Series'), axis=1)
realwage_f.head()
264 CHAPTER 17. PANDAS FOR PANEL DATA
[5 rows x 32 columns]
Similar to relational databases like SQL, pandas has built in methods to merge datasets to-
gether.
Using country information from WorldData.info, weβll add the continent of each country to
realwage_f with the merge function.
The CSV file can be found in pandas_panel/countries.csv and can be downloaded
here.
worlddata = pd.read_csv('https://github.com/QuantEcon/QuantEcon.lectures.code/raw/master/
[12]: βͺpandas_panel/countries.csv', sep=';')
worlddata.head()
[5 rows x 17 columns]
First, weβll select just the country and continent variables from worlddata and rename the
column to βCountryβ
worlddata = worlddata[['Country (en)', 'Continent']]
[13]: worlddata = worlddata.rename(columns={'Country (en)': 'Country'})
worlddata.head()
Country Continent
[13]: 0 Afghanistan Asia
1 Egypt Africa
2 Γ
land Islands Europe
3 Albania Europe
4 Algeria Africa
Time 2016-01-01
Country
Australia 12.98
Belgium 9.76
Brazil 1.24
Canada 8.48
Chile 1.91
[5 rows x 11 columns]
We can use either left, right, inner, or outer join to merge our datasets:
We will also need to specify where the country name is located in each dataframe, which will
be the key that is used to merge the dataframes βonβ.
266 CHAPTER 17. PANDAS FOR PANEL DATA
[5 rows x 13 columns]
Countries that appeared in realwage_f but not in worlddata will have NaN in the Conti-
nent column.
To check whether this has occurred, we can use .isnull() on the continent column and
filter the merged dataframe
merged[merged['Continent'].isnull()]
[16]:
2006-01-01 00:00:00 2007-01-01 00:00:00 2008-01-01 00:00:00 β¦ \
[16]: 247 3.42 3.74 3.87 β¦
247 0.23 0.45 0.39 β¦
247 1.50 1.64 1.71 β¦
[3 rows x 13 columns]
merged['Country'].map(missing_continents)
17 NaN
[17]: 23 NaN
32 NaN
100 NaN
38 NaN
17.4. MERGING DATAFRAMES AND FILLING NANS 267
108 NaN
41 NaN
225 NaN
53 NaN
58 NaN
45 NaN
68 NaN
233 NaN
86 NaN
88 NaN
91 NaN
247 Asia
117 NaN
122 NaN
123 NaN
138 NaN
153 NaN
151 NaN
174 NaN
175 NaN
247 Europe
247 Europe
198 NaN
200 NaN
227 NaN
241 NaN
240 NaN
Name: Country, dtype: object
merged[merged['Country'] == 'Korea']
[1 rows x 13 columns]
We will also combine the Americas into a single continent - this will make our visualization
nicer later on.
To do this, we will use .replace() and loop through a list of the continent values we want
to replace
replace = ['Central America', 'North America', 'South America']
[19]:
for country in replace:
merged['Continent'].replace(to_replace=country,
value='America',
inplace=True)
Now that we have all the data we want in a single DataFrame, we will reshape it back into
panel form with a MultiIndex.
We should also ensure to sort the index using .sort_index() so that we can efficiently fil-
ter our dataframe later on.
By default, levels will be sorted top-down
268 CHAPTER 17. PANDAS FOR PANEL DATA
2015-01-01 2016-01-01
Continent Country
America Brazil 1.21 1.24
Canada 8.35 8.48
Chile 1.81 1.91
Colombia 1.13 1.12
Costa Rica 2.56 2.63
[5 rows x 11 columns]
While merging, we lost our DatetimeIndex, as we merged columns that were not in date-
time format
merged.columns
[21]:
Index([2006-01-01 00:00:00, 2007-01-01 00:00:00, 2008-01-01 00:00:00,
[21]: 2009-01-01 00:00:00, 2010-01-01 00:00:00, 2011-01-01 00:00:00,
2012-01-01 00:00:00, 2013-01-01 00:00:00, 2014-01-01 00:00:00,
2015-01-01 00:00:00, 2016-01-01 00:00:00],
dtype='object')
Now that we have set the merged columns as the index, we can recreate a DatetimeIndex
using .to_datetime()
merged.columns = pd.to_datetime(merged.columns)
[22]: merged.columns = merged.columns.rename('Time')
merged.columns
The DatetimeIndex tends to work more smoothly in the row axis, so we will go ahead and
transpose merged
merged = merged.transpose()
[23]: merged.head()
[5 rows x 32 columns]
Grouping and summarizing data can be particularly useful for understanding large panel
datasets.
17.5. GROUPING AND SUMMARIZING DATA 269
A simple way to summarize data is to call an aggregation method on the dataframe, such as
.mean() or .max().
For example, we can calculate the average real minimum wage for each country over the pe-
riod 2006 to 2016 (the default is to aggregate over rows)
merged.mean().head(10)
[24]:
Continent Country
[24]: America Brazil 1.09
Canada 7.82
Chile 1.62
Colombia 1.07
Costa Rica 2.53
Mexico 0.53
United States 7.15
Asia Israel 5.95
Japan 6.18
Korea 4.22
dtype: float64
Using this series, we can plot the average real minimum wage over the past decade for each
country in our data set
import matplotlib.pyplot as plt
[25]: %matplotlib inline
import matplotlib
matplotlib.style.use('seaborn')
plt.show()
270 CHAPTER 17. PANDAS FOR PANEL DATA
Passing in axis=1 to .mean() will aggregate over columns (giving the average minimum
wage for all countries over time)
merged.mean(axis=1).head()
[26]:
Time
[26]: 2006-01-01 4.69
2007-01-01 4.84
2008-01-01 4.90
2009-01-01 5.08
2010-01-01 5.11
dtype: float64
We can also specify a level of the MultiIndex (in the column axis) to aggregate over
merged.mean(level='Continent', axis=1).head()
[28]:
Continent America Asia Australia Europe
[28]: Time
2006-01-01 2.80 4.29 10.25 4.80
2007-01-01 2.85 4.44 10.73 4.94
2008-01-01 2.99 4.45 10.76 4.99
2009-01-01 3.23 4.53 10.97 5.16
2010-01-01 3.34 4.53 10.95 5.17
We can plot the average minimum wages in each continent as a time series
merged.mean(level='Continent', axis=1).plot()
[29]: plt.title('Average real minimum wage')
plt.ylabel('2015 USD')
plt.xlabel('Year')
plt.show()
272 CHAPTER 17. PANDAS FOR PANEL DATA
The groupby method achieves the first step of this process, creating a new
DataFrameGroupBy object with data split into groups.
Letβs split merged by continent again, this time using the groupby function, and name the
resulting object grouped
grouped = merged.groupby(level='Continent', axis=1)
[32]: grouped
Calling an aggregation method on the object applies the function to each group, the results of
which are combined in a new data structure.
For example, we can return the number of countries in our dataset for each continent using
.size().
In this case, our new data structure is a Series
grouped.size()
[33]:
Continent
[33]: America 7
Asia 4
Europe 19
dtype: int64
Calling .get_group() to return just the countries in a single group, we can create a kernel
density estimate of the distribution of real minimum wages in 2016 for each continent.
grouped.groups.keys() will return the keys from the groupby object
import seaborn as sns
[34]:
continents = grouped.groups.keys()
plt.show()
This lecture has provided an introduction to some of pandasβ more advanced features, includ-
ing multiindices, merging, grouping and plotting.
Other tools that may be useful in panel data analysis include xarray, a python package that
extends pandas to N-dimensional data structures.
17.7 Exercises
17.7.1 Exercise 1
In these exercises, youβll work with a dataset of employment rates in Europe by age and sex
from Eurostat.
The dataset pandas_panel/employ.csv can be downloaded here.
Reading in the CSV file returns a panel dataset in long format. Use .pivot_table() to
construct a wide format dataframe with a MultiIndex in the columns.
Start off by exploring the dataframe and the variables available in the MultiIndex levels.
Write a program that quickly returns all values in the MultiIndex.
17.8. SOLUTIONS 275
17.7.2 Exercise 2
Filter the above dataframe to only include employment as a percentage of βactive populationβ.
Create a grouped boxplot using seaborn of employment rates in 2015 by age group and sex.
Hint: GEO includes both areas and countries.
17.8 Solutions
17.8.1 Exercise 1
employ = pd.read_csv('https://github.com/QuantEcon/QuantEcon.lectures.code/raw/master/
[35]: βͺpandas_panel/employ.csv')
employ = employ.pivot_table(values='Value',
index=['DATE'],
columns=['UNIT','AGE', 'SEX', 'INDIC_EM', 'GEO'])
employ.index = pd.to_datetime(employ.index) # ensure that dates are datetime format
employ.head()
UNIT
AGE
SEX
INDIC_EM
GEO United Kingdom
DATE
2007-01-01 4,131.00
2008-01-01 4,204.00
2009-01-01 4,193.00
2010-01-01 4,186.00
2011-01-01 4,164.00
This is a large dataset so it is useful to explore the levels and variables available
employ.columns.names
[36]:
FrozenList(['UNIT', 'AGE', 'SEX', 'INDIC_EM', 'GEO'])
[36]:
17.8.2 Exercise 2
To easily filter by country, swap GEO to the top level and sort the MultiIndex
employ.columns = employ.columns.swaplevel(0,-1)
[38]: employ = employ.sort_index(axis=1)
We need to get rid of a few items in GEO which are not countries.
A fast way to get rid of the EU areas is to use a list comprehension to find the level values in
GEO that begin with βEuroβ
geo_list = employ.columns.get_level_values('GEO').unique().tolist()
[39]: countries = [x for x in geo_list if not x.startswith('Euro')]
employ = employ[countries]
employ.columns.get_level_values('GEO').unique()
Select only percentage employed in the active population from the dataframe
employ_f = employ.xs(('Percentage of total population', 'Active population'),
[40]: level=('UNIT', 'INDIC_EM'),
axis=1)
employ_f.head()
GEO
AGE
SEX Total
DATE
2007-01-01 59.30
2008-01-01 59.80
2009-01-01 60.30
2010-01-01 60.00
2011-01-01 59.70
18.1 Contents
β’ Overview 18.2
β’ Endogeneity 18.5
β’ Summary 18.6
β’ Exercises 18.7
β’ Solutions 18.8
In addition to whatβs in Anaconda, this lecture will need the following libraries:
!pip install linearmodels
[1]:
18.2 Overview
Linear regression is a standard tool for analyzing the relationship between two or more vari-
ables.
In this lecture, weβll use the Python package statsmodels to estimate, interpret, and visu-
alize linear regression models.
Along the way, weβll discuss a variety of topics, including
As an example, we will replicate results from Acemoglu, Johnson and Robinsonβs seminal pa-
per [3].
279
280 CHAPTER 18. LINEAR REGRESSION IN PYTHON
In the paper, the authors emphasize the importance of institutions in economic development.
The main contribution is the use of settler mortality rates as a source of exogenous variation
in institutional differences.
Such variation is needed to determine whether it is institutions that give rise to greater eco-
nomic growth, rather than the other way around.
Letβs start with some imports:
import numpy as np
[2]: import matplotlib.pyplot as plt
%matplotlib inline
import pandas as pd
import statsmodels.api as sm
from statsmodels.iolib.summary2 import summary_col
from linearmodels.iv import IV2SLS
18.2.1 Prerequisites
18.2.2 Comments
[3] wish to determine whether or not differences in institutions can help to explain observed
economic outcomes.
How do we measure institutional differences and economic outcomes?
In this paper,
β’ economic outcomes are proxied by log GDP per capita in 1995, adjusted for exchange
rates.
β’ institutional differences are proxied by an index of protection against expropriation on
average over 1985-95, constructed by the Political Risk Services Group.
These variables and other data used in the paper are available for download on Daron Ace-
mogluβs webpage.
We will use pandasβ .read_stata() function to read in data contained in the .dta files to
dataframes
df1 = pd.read_stata('https://github.com/QuantEcon/QuantEcon.lectures.code/raw/master/ols/
[3]: βͺmaketable1.dta')
df1.head()
18.3. SIMPLE LINEAR REGRESSION 281
Letβs use a scatterplot to see whether any obvious relationship exists between GDP per capita
and the protection against expropriation index
plt.style.use('seaborn')
[4]:
df1.plot(x='avexpr', y='logpgp95', kind='scatter')
plt.show()
The plot shows a fairly strong positive relationship between protection against expropriation
and log GDP per capita.
Specifically, if higher protection against expropriation is a measure of institutional quality,
then better institutions appear to be positively correlated with better economic outcomes
(higher GDP per capita).
Given the plot, choosing a linear model to describe this relationship seems like a reasonable
assumption.
We can write our model as
282 CHAPTER 18. LINEAR REGRESSION IN PYTHON
ππππππ95π = π½0 + π½1 ππ£ππ₯πππ + π’π
where:
Visually, this linear model involves choosing a straight line that best fits the data, as in the
following plot (Figure 2 in [3])
# Dropping NA's is required to use numpy's polyfit
[5]: df1_subset = df1.dropna(subset=['logpgp95', 'avexpr'])
X = df1_subset['avexpr']
y = df1_subset['logpgp95']
labels = df1_subset['shortnam']
ax.set_xlim([3.3,10.5])
ax.set_ylim([4,10.5])
ax.set_xlabel('Average Expropriation Risk 1985-95')
ax.set_ylabel('Log GDP per capita, PPP, 1995')
ax.set_title('Figure 2: OLS relationship between expropriation \
risk and income')
plt.show()
18.3. SIMPLE LINEAR REGRESSION 283
The most common technique to estimate the parameters (π½βs) of the linear model is Ordinary
Least Squares (OLS).
As the name implies, an OLS model is solved by finding the parameters that minimize the
sum of squared residuals, i.e.
π
min β π’Μ2π
π½Μ π=1
where π’Μπ is the difference between the observation and the predicted value of the dependent
variable.
To estimate the constant term π½0 , we need to add a column of 1βs to our dataset (consider
the equation if π½0 was replaced with π½0 π₯π and π₯π = 1)
df1['const'] = 1
[6]:
Now we can construct our model in statsmodels using the OLS function.
We will use pandas dataframes with statsmodels, however standard arrays can also be
used as arguments
reg1 = sm.OLS(endog=df1['logpgp95'], exog=df1[['const', 'avexpr']], \
[7]: missing='drop')
type(reg1)
statsmodels.regression.linear_model.OLS
[7]:
results = reg1.fit()
[8]: type(results)
statsmodels.regression.linear_model.RegressionResultsWrapper
[8]:
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
Using our parameter estimates, we can now write our estimated relationship as
Μ
ππππππ95 π = 4.63 + 0.53 ππ£ππ₯πππ
This equation describes the line that best fits our data, as shown in Figure 2.
We can use this equation to predict the level of log GDP per capita for a value of the index of
expropriation protection.
18.3. SIMPLE LINEAR REGRESSION 285
For example, for a country with an index value of 7.07 (the average for the dataset), we find
that their predicted level of log GDP per capita in 1995 is 8.38.
mean_expr = np.mean(df1_subset['avexpr'])
[10]: mean_expr
6.515625
[10]:
predicted_logpdp95 = 4.63 + 0.53 * 7.07
[11]: predicted_logpdp95
8.3771
[11]:
An easier (and more accurate) way to obtain this result is to use .predict() and set
ππππ π‘πππ‘ = 1 and ππ£ππ₯πππ = ππππ_ππ₯ππ
results.predict(exog=[1, mean_expr])
[12]:
array([8.09156367])
[12]:
We can obtain an array of predicted ππππππ95π for every value of ππ£ππ₯πππ in our dataset by
calling .predict() on our results.
Plotting the predicted values against ππ£ππ₯πππ shows that the predicted values lie along the
linear line that we fitted above.
The observed values of ππππππ95π are also plotted for comparison purposes
# Drop missing observations from whole sample
[13]:
df1_plot = df1.dropna(subset=['logpgp95', 'avexpr'])
fix, ax = plt.subplots()
ax.scatter(df1_plot['avexpr'], results.predict(), alpha=0.5,
label='predicted')
ax.legend()
ax.set_title('OLS predicted values')
ax.set_xlabel('avexpr')
ax.set_ylabel('logpgp95')
plt.show()
286 CHAPTER 18. LINEAR REGRESSION IN PYTHON
So far we have only accounted for institutions affecting economic performance - almost cer-
tainly there are numerous other factors affecting GDP that are not included in our model.
Leaving out variables that affect ππππππ95π will result in omitted variable bias, yielding
biased and inconsistent parameter estimates.
We can extend our bivariate regression model to a multivariate regression model by
adding in other factors that may affect ππππππ95π .
[3] consider other factors such as:
Letβs estimate some of the extended models considered in the paper (Table 2) using data from
maketable2.dta
df2 = pd.read_stata('https://github.com/QuantEcon/QuantEcon.lectures.code/raw/master/ols/
[14]: βͺmaketable2.dta')
Now that we have fitted our model, we will use summary_col to display the results in a sin-
gle table (model numbers correspond to those in the paper)
info_dict={'R-squared' : lambda x: f"{x.rsquared:.2f}",
[15]: 'No. observations' : lambda x: f"{int(x.nobs):d}"}
results_table = summary_col(results=[reg1,reg2,reg3],
float_format='%0.2f',
stars = True,
model_names=['Model 1',
'Model 3',
'Model 4'],
info_dict=info_dict,
regressor_order=['const',
'avexpr',
'lat_abst',
'asia',
'africa'])
print(results_table)
18.5 Endogeneity
As [3] discuss, the OLS models likely suffer from endogeneity issues, resulting in biased and
inconsistent model estimates.
Namely, there is likely a two-way relationship between institutions and economic outcomes:
To deal with endogeneity, we can use two-stage least squares (2SLS) regression, which
is an extension of OLS regression.
This method requires replacing the endogenous variable ππ£ππ₯πππ with a variable that is:
The new set of regressors is called an instrument, which aims to remove endogeneity in our
proxy of institutional differences.
The main contribution of [3] is the use of settler mortality rates to instrument for institu-
tional differences.
They hypothesize that higher mortality rates of colonizers led to the establishment of insti-
tutions that were more extractive in nature (less protection against expropriation), and these
institutions still persist today.
Using a scatterplot (Figure 3 in [3]), we can see protection against expropriation is negatively
correlated with settler mortality rates, coinciding with the authorsβ hypothesis and satisfying
the first condition of a valid instrument.
# Dropping NA's is required to use numpy's polyfit
[16]: df1_subset2 = df1.dropna(subset=['logem4', 'avexpr'])
X = df1_subset2['logem4']
y = df1_subset2['avexpr']
labels = df1_subset2['shortnam']
ax.set_xlim([1.8,8.4])
ax.set_ylim([3.3,10.4])
ax.set_xlabel('Log of Settler Mortality')
ax.set_ylabel('Average Expropriation Risk 1985-95')
ax.set_title('Figure 3: First-stage relationship between settler mortality \
and expropriation risk')
plt.show()
18.5. ENDOGENEITY 289
The second condition may not be satisfied if settler mortality rates in the 17th to 19th cen-
turies have a direct effect on current GDP (in addition to their indirect effect through institu-
tions).
For example, settler mortality rates may be related to the current disease environment in a
country, which could affect current economic performance.
[3] argue this is unlikely because:
β’ The majority of settler deaths were due to malaria and yellow fever and had a limited
effect on local people.
β’ The disease burden on local people in Africa or India, for example, did not appear to
be higher than average, supported by relatively high population densities in these areas
before colonization.
As we appear to have a valid instrument, we can use 2SLS regression to obtain consistent and
unbiased parameter estimates.
First stage
The first stage involves regressing the endogenous variable (ππ£ππ₯πππ ) on the instrument.
The instrument is the set of all exogenous variables in our model (and not just the variable
we have replaced).
Using model 1 as an example, our instrument is simply a constant and settler mortality rates
πππππ4π .
Therefore, we will estimate the first-stage regression as
ππ£ππ₯πππ = πΏ0 + πΏ1 πππππ4π + π£π
290 CHAPTER 18. LINEAR REGRESSION IN PYTHON
The data we need to estimate this equation is located in maketable4.dta (only complete
data, indicated by baseco = 1, is used for estimation)
# Import and select the data
[17]: df4 = pd.read_stata('https://github.com/QuantEcon/QuantEcon.lectures.code/raw/master/ols/
βͺmaketable4.dta')
df4 = df4[df4['baseco'] == 1]
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
Second stage
We need to retrieve the predicted values of ππ£ππ₯πππ using .predict().
We then replace the endogenous variable ππ£ππ₯πππ with the predicted values ππ£ππ₯ππ
Μ π in the
original linear model.
Our second stage regression is thus
ππππππ95π = π½0 + π½1 ππ£ππ₯ππ
Μ π + π’π
df4['predicted_avexpr'] = results_fs.predict()
[18]:
results_ss = sm.OLS(df4['logpgp95'],
df4[['const', 'predicted_avexpr']]).fit()
print(results_ss.summary())
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
The second-stage regression results give us an unbiased and consistent estimate of the effect
of institutions on economic outcomes.
The result suggests a stronger positive relationship than what the OLS results indicated.
Note that while our parameter estimates are correct, our standard errors are not and for this
reason, computing 2SLS βmanuallyβ (in stages with OLS) is not recommended.
We can correctly estimate a 2SLS regression in one step using the linearmodels package, an
extension of statsmodels
Note that when using IV2SLS, the exogenous and instrument variables are split up in the
function arguments (whereas before the instrument included exogenous variables)
iv = IV2SLS(dependent=df4['logpgp95'],
[19]: exog=df4['const'],
endog=df4['avexpr'],
instruments=df4['logem4']).fit(cov_type='unadjusted')
print(iv.summary)
Parameter Estimates
==============================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
------------------------------------------------------------------------------
const 1.9097 1.0106 1.8897 0.0588 -0.0710 3.8903
avexpr 0.9443 0.1541 6.1293 0.0000 0.6423 1.2462
==============================================================================
Endogenous: avexpr
Instruments: logem4
292 CHAPTER 18. LINEAR REGRESSION IN PYTHON
Given that we now have consistent and unbiased estimates, we can infer from the model we
have estimated that institutional differences (stemming from institutions set up during colo-
nization) can help to explain differences in income levels across countries today.
[3] use a marginal effect of 0.94 to calculate that the difference in the index between Chile
and Nigeria (ie. institutional quality) implies up to a 7-fold difference in income, emphasizing
the significance of institutions in economic development.
18.6 Summary
18.7 Exercises
18.7.1 Exercise 1
In the lecture, we think the original model suffers from endogeneity bias due to the likely ef-
fect income has on institutional development.
Although endogeneity is often best identified by thinking about the data and model, we can
formally test for endogeneity using the Hausman test.
We want to test for correlation between the endogenous variable, ππ£ππ₯πππ , and the errors, π’π
ππ£ππ₯πππ = π0 + π1 πππππ4π + ππ
Second, we retrieve the residuals ππΜ and include them in the original equation
If πΌ is statistically significant (with a p-value < 0.05), then we reject the null hypothesis and
conclude that ππ£ππ₯πππ is endogenous.
Using the above information, estimate a Hausman test and interpret your results.
18.8. SOLUTIONS 293
18.7.2 Exercise 2
The OLS parameter π½ can also be estimated using matrix algebra and numpy (you may need
to review the numpy lecture to complete this exercise).
The linear equation we want to estimate is (written in matrix form)
π¦ = ππ½ + π’
To solve for the unknown parameter π½, we want to minimize the sum of squared residuals
minπ’Μβ² π’Μ
π½Μ
Rearranging the first equation and substituting into the second equation, we can write
Solving this optimization problem gives the solution for the π½ Μ coefficients
π½ Μ = (π β² π)β1 π β² π¦
Using the above information, compute π½ Μ from model 1 using numpy - your results should be
the same as those in the statsmodels output from earlier in the lecture.
18.8 Solutions
18.8.1 Exercise 1
# Load in data
[20]: df4 = pd.read_stata('https://github.com/QuantEcon/QuantEcon.lectures.code/raw/master/ols/
βͺmaketable4.dta')
print(reg2.summary())
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
The output shows that the coefficient on the residuals is statistically significant, indicating
ππ£ππ₯πππ is endogenous.
18.8.2 Exercise 2
# Load in data
[21]: df1 = pd.read_stata('https://github.com/QuantEcon/QuantEcon.lectures.code/raw/master/ols/
βͺmaketable1.dta')
df1 = df1.dropna(subset=['logpgp95', 'avexpr'])
# Compute Ξ²_hat
Ξ²_hat = np.linalg.solve(X.T @ X, X.T @ y)
Ξ²_0 = 4.6
Ξ²_1 = 0.53
19.1 Contents
β’ Overview 19.2
β’ Summary 19.8
β’ Exercises 19.9
β’ Solutions 19.10
19.2 Overview
In a previous lecture, we estimated the relationship between dependent and explanatory vari-
ables using linear regression.
But what if a linear relationship is not an appropriate assumption for our model?
One widely used alternative is maximum likelihood estimation, which involves specifying a
class of distributions, indexed by unknown parameters, and then using the data to pin down
these parameter values.
The benefit relative to linear regression is that it allows more flexibility in the probabilistic
relationships between variables.
Here we illustrate maximum likelihood by replicating Daniel Treismanβs (2016) paper, Rus-
siaβs Billionaires, which connects the number of billionaires in a country to its economic char-
acteristics.
The paper concludes that Russia has a higher number of billionaires than economic factors
such as market size and tax rate predict.
295
296 CHAPTER 19. MAXIMUM LIKELIHOOD ESTIMATION
19.2.1 Prerequisites
19.2.2 Comments
Letβs consider the steps we need to go through in maximum likelihood estimation and how
they pertain to this study.
The first step with maximum likelihood estimation is to choose the probability distribution
believed to be generating the data.
More precisely, we need to make an assumption as to which parametric class of distributions
is generating the data.
β’ e.g., the class of all normal distributions, or the class of all gamma distributions.
β’ e.g., the class of normal distributions is a family of distributions indexed by its mean
π β (ββ, β) and standard deviation π β (0, β).
Weβll let the data pick out a particular element of the class by pinning down the parameters.
The parameter estimates so produced will be called maximum likelihood estimates.
Hence we consider distributions that take values only in the nonnegative integers.
(This is one reason least squares regression is not the best tool for the present problem, since
the dependent variable in linear regression is not restricted to integer values)
One integer distribution is the Poisson distribution, the probability mass function (pmf) of
which is
ππ¦ βπ
π(π¦) = π , π¦ = 0, 1, 2, β¦ , β
π¦!
We can plot the Poisson distribution over π¦ for different values of π as follows
poisson_pmf = lambda y, ΞΌ: ΞΌ**y / factorial(y) * exp(-ΞΌ)
[2]: y_values = range(0, 25)
ax.grid()
ax.set_xlabel('$y$', fontsize=14)
ax.set_ylabel('$f(y \mid \mu)$', fontsize=14)
ax.axis(xmin=0, ymin=0)
ax.legend(fontsize=14)
plt.show()
298 CHAPTER 19. MAXIMUM LIKELIHOOD ESTIMATION
Notice that the Poisson distribution begins to resemble a normal distribution as the mean of
π¦ increases.
Letβs have a look at the distribution of the data weβll be working with in this lecture.
Treismanβs main source of data is Forbesβ annual rankings of billionaires and their estimated
net worth.
The dataset mle/fp.dta can be downloaded here or from its AER page.
pd.options.display.max_columns = 10
[3]:
# Load in data and view
df = pd.read_stata('https://github.com/QuantEcon/QuantEcon.lectures.code/raw/master/mle/fp.
βͺdta')
df.head()
[5 rows x 36 columns]
Using a histogram, we can view the distribution of the number of billionaires per country,
numbil0, in 2008 (the United States is dropped for plotting purposes)
numbil0_2008 = df[(df['year'] == 2008) & (
[4]: df['country'] != 'United States')].loc[:, 'numbil0']
plt.subplots(figsize=(12, 8))
plt.hist(numbil0_2008, bins=30)
plt.xlim(left=0)
plt.grid()
plt.xlabel('Number of billionaires in 2008')
plt.ylabel('Count')
plt.show()
19.4. CONDITIONAL DISTRIBUTIONS 299
From the histogram, it appears that the Poisson assumption is not unreasonable (albeit with
a very low π and some outliers).
π¦
ππ π βππ
π(π¦π β£ xπ ) = π ; π¦π = 0, 1, 2, β¦ , β. (1)
π¦π !
To illustrate the idea that the distribution of π¦π depends on xπ letβs run a simple simulation.
We use our poisson_pmf function from above and arbitrary values for π½ and xπ
y_values = range(0, 20)
[5]:
# Define a parameter vector with estimates
Ξ² = np.array([0.26, 0.18, 0.25, -0.1, -0.22])
np.array([6, 5, 4, 4, 7])]
for X in datasets:
ΞΌ = exp(X @ Ξ²)
distribution = []
for y_i in y_values:
distribution.append(poisson_pmf(y_i, ΞΌ))
ax.plot(y_values,
distribution,
label=f'$\mu_i$={ΞΌ:.1}',
marker='o',
markersize=8,
alpha=0.5)
ax.grid()
ax.legend()
ax.set_xlabel('$y \mid x_i$')
ax.set_ylabel(r'$f(y \mid x_i; \beta )$')
ax.axis(xmin=0, ymin=0)
plt.show()
In our model for number of billionaires, the conditional distribution contains 4 (π = 4) pa-
rameters that we need to estimate.
We will label our entire parameter vector as π½ where
19.5. MAXIMUM LIKELIHOOD ESTIMATION 301
π½0
β‘π½ β€
π½ = β’ 1β₯
β’π½2 β₯
β£π½3 β¦
To estimate the model using MLE, we want to maximize the likelihood that our estimate π½Μ is
the true parameter π½.
Intuitively, we want to find the π½Μ that best fits our data.
First, we need to construct the likelihood function β(π½), which is similar to a joint probabil-
ity density function.
Assume we have some data π¦π = {π¦1 , π¦2 } and π¦π βΌ π(π¦π ).
If π¦1 and π¦2 are independent, the joint pmf of these data is π(π¦1 , π¦2 ) = π(π¦1 ) β
π(π¦2 ).
If π¦π follows a Poisson distribution with π = 7, we can visualize the joint pmf like so
def plot_joint_poisson(ΞΌ=7, y_n=20):
[6]: yi_values = np.arange(0, y_n, 1)
plot_joint_poisson(ΞΌ=7, y_n=20)
302 CHAPTER 19. MAXIMUM LIKELIHOOD ESTIMATION
Similarly, the joint pmf of our data (which is distributed as a conditional Poisson distribu-
tion) can be written as
π π¦
ππ π βππ
π(π¦1 , π¦2 , β¦ , π¦π β£ x1 , x2 , β¦ , xπ ; π½) = β π
π=1
π¦π !
π π¦
ππ π βππ
β(π½ β£ π¦1 , π¦2 , β¦ , π¦π ; x1 , x2 , β¦ , xπ ) = β π
π=1
π¦π !
=π(π¦1 , π¦2 , β¦ , π¦π β£ x1 , x2 , β¦ , xπ ; π½)
Now that we have our likelihood function, we want to find the π½Μ that yields the maximum
likelihood value
maxβ(π½)
π½
The MLE of the Poisson to the Poisson for π½ Μ can be obtained by solving
π π π
max( β π¦π log ππ β β ππ β β log π¦!)
π½
π=1 π=1 π=1
However, no analytical solution exists to the above problem β to find the MLE we need to use
numerical methods.
19.6. MLE WITH NUMERICAL METHODS 303
Many distributions do not have nice, analytical solutions and therefore require numerical
methods to solve for parameter estimates.
One such numerical method is the Newton-Raphson algorithm.
Our goal is to find the maximum likelihood estimate π½.Μ
At π½,Μ the first derivative of the log-likelihood function will be equal to 0.
Letβs illustrate this by supposing
Ξ² = np.linspace(1, 20)
[7]: logL = -(Ξ² - 10) ** 2 - 10
dlogL = -2 * Ξ² + 20
ax1.set_ylabel(r'$log \mathcal{L(\beta)}$',
rotation=0,
labelpad=35,
fontsize=15)
ax2.set_ylabel(r'$\frac{dlog \mathcal{L(\beta)}}{d \beta}$ ',
rotation=0,
labelpad=35,
fontsize=19)
ax2.set_xlabel(r'$\beta$', fontsize=15)
ax1.grid(), ax2.grid()
plt.axhline(c='black')
plt.show()
304 CHAPTER 19. MAXIMUM LIKELIHOOD ESTIMATION
π log β(π½)
The plot shows that the maximum likelihood value (the top plot) occurs when ππ½ = 0
(the bottom plot).
Therefore, the likelihood is maximized when π½ = 10.
We can also ensure that this value is a maximum (as opposed to a minimum) by checking
that the second derivative (slope of the bottom plot) is negative.
The Newton-Raphson algorithm finds a point where the first derivative is 0.
To use the algorithm, we take an initial guess at the maximum value, π½0 (the OLS parameter
estimates might be a reasonable guess), then
As can be seen from the updating equation, π½ (π+1) = π½ (π) only when πΊ(π½ (π) ) = 0 ie. where the
first derivative is equal to 0.
(In practice, we stop iterating when the difference is below a small tolerance threshold)
Letβs have a go at implementing the Newton-Raphson algorithm.
First, weβll create a class called PoissonRegression so we can easily recompute the values
of the log likelihood, gradient and Hessian for every iteration
class PoissonRegression:
[8]:
def __init__(self, y, X, Ξ²):
self.X = X
self.n, self.k = X.shape
# Reshape y as a n_by_1 column vector
self.y = y.reshape(self.n,1)
# Reshape Ξ² as a k_by_1 column vector
self.Ξ² = Ξ².reshape(self.k,1)
def ΞΌ(self):
return np.exp(self.X @ self.Ξ²)
def logL(self):
y = self.y
ΞΌ = self.ΞΌ()
return np.sum(y * np.log(ΞΌ) - ΞΌ - np.log(factorial(y)))
def G(self):
y = self.y
ΞΌ = self.ΞΌ()
return X.T @ (y - ΞΌ)
19.6. MLE WITH NUMERICAL METHODS 305
def H(self):
X = self.X
ΞΌ = self.ΞΌ()
return -(X.T @ (ΞΌ * X))
Our function newton_raphson will take a PoissonRegression object that has an initial
guess of the parameter vector π½ 0 .
The algorithm will update the parameter vector according to the updating rule, and recalcu-
late the gradient and Hessian matrices at the new parameter estimates.
Iteration will end when either:
β’ The difference between the parameter and the updated parameter is below a tolerance
level.
β’ The maximum number of iterations has been achieved (meaning convergence is not
achieved).
So we can get an idea of whatβs going on while the algorithm is running, an option
display=True is added to print out values at each iteration.
def newton_raphson(model, tol=1e-3, max_iter=1000, display=True):
[9]:
i = 0
error = 100 # Initial error value
# Print iterations
if display:
Ξ²_list = [f'{t:.3}' for t in list(model.Ξ².flatten())]
update = f'{i:<13}{model.logL():<16.8}{Ξ²_list}'
print(update)
i += 1
Letβs try out our algorithm with a small dataset of 5 observations and 3 variables in X.
X = np.array([[1, 2, 5],
[10]: [1, 1, 3],
[1, 4, 2],
[1, 5, 2],
[1, 3, 1]])
y = np.array([1, 0, 1, 1, 0])
Iteration_k Log-likelihood ΞΈ
--------------------------------------------------------------------------------
---------
0 -4.3447622 ['-1.49', '0.265', '0.244']
1 -3.5742413 ['-3.38', '0.528', '0.474']
2 -3.3999526 ['-5.06', '0.782', '0.702']
3 -3.3788646 ['-5.92', '0.909', '0.82']
4 -3.3783559 ['-6.07', '0.933', '0.843']
5 -3.3783555 ['-6.08', '0.933', '0.843']
Number of iterations: 6
Ξ²_hat = [-6.07848205 0.93340226 0.84329625]
As this was a simple model with few observations, the algorithm achieved convergence in only
6 iterations.
You can see that with each iteration, the log-likelihood value increased.
Remember, our objective was to maximize the log-likelihood function, which the algorithm
has worked to achieve.
Also, note that the increase in log β(π½ (π) ) becomes smaller with each iteration.
This is because the gradient is approaching 0 as we reach the maximum, and therefore the
numerator in our updating equation is becoming smaller.
The gradient vector should be close to 0 at π½Μ
poi.G()
[11]:
array([[-3.95169226e-07],
[11]: [-1.00114804e-06],
[-7.73114559e-07]])
The iterative process can be visualized in the following diagram, where the maximum is found
at π½ = 10
logL = lambda x: -(x - 10) ** 2 - 10
[12]:
def find_tangent(Ξ², a=0.01):
y1 = logL(Ξ²)
y2 = logL(Ξ²+a)
x = np.array([[Ξ², 1], [Ξ²+a, 1]])
m, c = np.linalg.lstsq(x, np.array([y1, y2]), rcond=None)[0]
return m, c
Ξ² = np.linspace(2, 18)
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(Ξ², logL(Ξ²), lw=2, c='black')
labelpad=25,
fontsize=15)
ax.grid(alpha=0.3)
plt.show()
Note that our implementation of the Newton-Raphson algorithm is rather basic β for more
robust implementations see, for example, scipy.optimize.
Now that we know whatβs going on under the hood, we can apply MLE to an interesting ap-
plication.
Weβll use the Poisson regression model in statsmodels to obtain a richer output with stan-
dard errors, test values, and more.
statsmodels uses the same algorithm as above to find the maximum likelihood estimates.
Before we begin, letβs re-estimate our simple model with statsmodels to confirm we obtain
the same coefficients and log-likelihood value.
X = np.array([[1, 2, 5],
[13]: [1, 1, 3],
[1, 4, 2],
[1, 5, 2],
[1, 3, 1]])
y = np.array([1, 0, 1, 1, 0])
Iterations 7
Poisson Regression Results
==============================================================================
Dep. Variable: y No. Observations: 5
Model: Poisson Df Residuals: 2
Method: MLE Df Model: 2
Date: Sun, 20 Oct 2019 Pseudo R-squ.: 0.2546
Time: 17:06:51 Log-Likelihood: -3.3784
converged: True LL-Null: -4.5325
Covariance Type: nonrobust LLR p-value: 0.3153
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const -6.0785 5.279 -1.151 0.250 -16.425 4.268
x1 0.9334 0.829 1.126 0.260 -0.691 2.558
x2 0.8433 0.798 1.057 0.291 -0.720 2.407
==============================================================================
Now letβs replicate results from Daniel Treismanβs paper, Russiaβs Billionaires, mentioned ear-
lier in the lecture.
Treisman starts by estimating equation Eq. (1), where:
β’ π¦π is ππ’ππππ ππ ππππππππππππ π
β’ π₯π1 is log πΊπ·π πππ πππππ‘ππ
β’ π₯π2 is log ππππ’πππ‘ππππ
β’ π₯π3 is π¦ππππ ππ πΊπ΄π π π β years membership in GATT and WTO (to proxy access to in-
ternational markets)
# Add a constant
df['const'] = 1
# Variable sets
reg1 = ['const', 'lngdppc', 'lnpop', 'gattwto08']
reg2 = ['const', 'lngdppc', 'lnpop',
'gattwto08', 'lnmcap08', 'rintr', 'topint08']
reg3 = ['const', 'lngdppc', 'lnpop', 'gattwto08', 'lnmcap08',
'rintr', 'topint08', 'nrrents', 'roflaw']
Then we can use the Poisson function from statsmodels to fit the model.
Weβll use robust standard errors as in the authorβs paper
# Specify model
[15]: poisson_reg = sm.Poisson(df[['numbil0']], df[reg1],
missing='drop').fit(cov_type='HC0')
print(poisson_reg.summary())
results_table = summary_col(results=results,
float_format='%0.3f',
stars=True,
model_names=reg_names,
info_dict=info_dict,
regressor_order=regressor_order)
results_table.add_title('Table 1 - Explaining the Number of Billionaires \
in 2008')
print(results_table)
nrrents -0.005
(0.010)
roflaw 0.203
(0.372)
Pseudo R-squared 0.86 0.90 0.90
No. observations 197 131 131
=================================================
Standard errors in parentheses.
* p<.1, ** p<.05, ***p<.01
The output suggests that the frequency of billionaires is positively correlated with GDP
per capita, population size, stock market capitalization, and negatively correlated with top
marginal income tax rate.
To analyze our results by country, we can plot the difference between the predicted an actual
values, then sort from highest to lowest and plot the first 15
data = ['const', 'lngdppc', 'lnpop', 'gattwto08', 'lnmcap08', 'rintr',
[17]: 'topint08', 'nrrents', 'roflaw', 'numbil0', 'country']
results_df = df[data].dropna()
# Calculate difference
results_df['difference'] = results_df['numbil0'] - results_df['prediction']
As we can see, Russia has by far the highest number of billionaires in excess of what is pre-
dicted by the model (around 50 more than expected).
Treisman uses this empirical result to discuss possible reasons for Russiaβs excess of billion-
aires, including the origination of wealth in Russia, the political climate, and the history of
privatization in the years after the USSR.
19.8 Summary
19.9 Exercises
19.9.1 Exercise 1
Suppose we wanted to estimate the probability of an event π¦π occurring, given some observa-
tions.
312 CHAPTER 19. MAXIMUM LIKELIHOOD ESTIMATION
π¦
π(π¦π ; π½) = ππ π (1 β ππ )1βπ¦π , π¦π = 0, 1
where ππ = Ξ¦(xβ²π π½)
Ξ¦ represents the cumulative normal distribution and constrains the predicted π¦π to be be-
tween 0 and 1 (as required for a probability).
π½ is a vector of coefficients.
Following the example in the lecture, write a class to represent the Probit model.
To begin, find the log-likelihood function and derive the gradient and Hessian.
The scipy module stats.norm contains the functions needed to compute the cmf and pmf
of the normal distribution.
19.9.2 Exercise 2
Use the following dataset and initial values of π½ to estimate the MLE with the Newton-
Raphson algorithm developed earlier in the lecture
1 2 4 1
β‘1 1 1β€ β‘0β€ 0.1
β’ β₯ β’ β₯
X = β’1 4 3β₯ π¦ = β’1β₯ π½ (0) = β‘
β’0.1β₯
β€
β’1 5 6β₯ β’1β₯ β£0.1β¦
β£1 3 5β¦ β£0β¦
Verify your results with statsmodels - you can import the Probit function with the follow-
ing import statement
from statsmodels.discrete.discrete_model import Probit
[18]:
Note that the simple Newton-Raphson algorithm developed in this lecture is very sensitive to
initial values, and therefore you may fail to achieve convergence with different starting values.
19.10 Solutions
19.10.1 Exercise 1
π
log β = β [π¦π log Ξ¦(xβ²π π½) + (1 β π¦π ) log(1 β Ξ¦(xβ²π π½))]
π=1
π
Ξ¦(π ) = π(π )
ππ
where π is the marginal normal distribution.
19.10. SOLUTIONS 313
π
π log β π(xβ²π π½) π(xβ²π π½)
= β [π¦π β (1 β π¦π ) ]x
ππ½ π=1
β²
Ξ¦(xπ π½) 1 β Ξ¦(xβ²π π½) π
π
π 2 log β β² π(xβ²π π½) + xβ²π π½Ξ¦(xβ²π π½) ππ (xβ²π π½) β xβ²π π½(1 β Ξ¦(xβ²π π½))
β² = β β π(x π π½)[π¦ π β² 2
+ (1 β π¦ π ) β² 2
]xπ xβ²π
ππ½ππ½ π=1
[Ξ¦(x π π½)] [1 β Ξ¦(x π π½)]
Using these results, we can write a class for the Probit model as follows
class ProbitRegression:
[19]:
def __init__(self, y, X, Ξ²):
self.X, self.y, self.Ξ² = X, y, Ξ²
self.n, self.k = X.shape
def ΞΌ(self):
return norm.cdf(self.X @ self.Ξ².T)
def οΏ½(self):
return norm.pdf(self.X @ self.Ξ².T)
def logL(self):
ΞΌ = self.ΞΌ()
return np.sum(y * np.log(ΞΌ) + (1 - y) * np.log(1 - ΞΌ))
def G(self):
ΞΌ = self.ΞΌ()
οΏ½ = self.οΏ½()
return np.sum((X.T * y * οΏ½ / ΞΌ - X.T * (1 - y) * οΏ½ / (1 - ΞΌ)),
axis=1)
def H(self):
X = self.X
Ξ² = self.Ξ²
ΞΌ = self.ΞΌ()
οΏ½ = self.οΏ½()
a = (οΏ½ + (X @ Ξ².T) * ΞΌ) / ΞΌ**2
b = (οΏ½ - (X @ Ξ².T) * (1 - ΞΌ)) / (1 - ΞΌ)**2
return -(οΏ½ * (y * a + (1 - y) * b) * X.T) @ X
19.10.2 Exercise 2
X = np.array([[1, 2, 4],
[20]: [1, 1, 1],
[1, 4, 3],
[1, 5, 6],
[1, 3, 5]])
y = np.array([1, 0, 1, 1, 0])
Iteration_k Log-likelihood ΞΈ
--------------------------------------------------------------------------------
---------
0 -2.3796884 ['-1.34', '0.775', '-0.157']
1 -2.3687526 ['-1.53', '0.775', '-0.0981']
2 -2.3687294 ['-1.55', '0.778', '-0.0971']
3 -2.3687294 ['-1.55', '0.778', '-0.0971']
Number of iterations: 4
Ξ²_hat = [-1.54625858 0.77778952 -0.09709757]