Basic Econometrics with Python

Fabián Alejandro Triana Alarcón

National University of Colombia

Faculty of Economic Sciences
2019

This work is licensed under the Creative Commons Attribution-NonCommercial 4.0 International
 Preface
This is a manual for conducting basic econometric analysis using the language of
Python programming; does not understand, as part of its content, a solid theoretical body that can be taken
as an adequate basis for teaching Econometrics and is merely limited to being a practical reference for
application of the techniques and models specific to the subject. Given the applied approach of this work, the reader
you will not find in this a text that serves as a theoretical guide and allows you to strengthen your knowledge about the
Econometrics, as it is not the objective proposed here; for that purpose, it is suggested to refer to the manuals.
generally specialized employees used as reference material in undergraduate Economics courses.
The content presented in this work is particularly aimed at students of Economics, Statistics-
and related fields, to professionals whose work is associated with the field of Econometrics and to all
a person, with some prior knowledge, who feels affinity for the subject and is interested in improving their skills
skills in the use of useful computer tools in the development of applied econometrics exercises.
Given the minimum level of depth with which the theoretical foundations on which it is structured are addressed
Given the econometric analysis in this work, and its marked practical approach, it is expected that the reader possesses some.
good prior knowledge (at an introductory level) about the subject, so that it does not seem complex to him
follow the proposed exercises and be able to achieve your learning objectives.
In addition to basic prior knowledge of Econometrics, it is expected that the reader is familiar with
with the use of computer tools, in such a way that understands, without major problems, the instructions
provided for the development of the proposed exercises and the instructions given for installation and use
of the programs used. Likewise, it is necessary for the reader to have some knowledge of the language
English, as most functions, methods, and attributes used are represented through
words in English, which, given the very high level of the programming language that will be used, result
very informative and greatly facilitate the understanding of the tasks being carried out.
I wish to highlight my interest in having this work serve as a non-theoretical reference and facilitate learning.
I give Econometrics due to the practical treatment it aims to achieve; likewise, that it allows for the approach
of all persons interested in the subject, with a certain level of knowledge, regardless of the field of their training
formation, and promotes the use of technological tools by everyone. Although true
Previous knowledge or experience with programming languages is welcome and will be very useful for
the understanding of the work by the reader is essential because of the treatment given to the
information in this work, I aim for the presentations to be as illustrative and clear as possible,
so that every reader feels comfortable with the study of this work and achieves their learning goals
satisfactorily.

The Author,
Fabián Alejandro Triana Alarcón
General index

Installation of the programs and preparation of the environment 1

Python. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
Anaconda. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
Jupyter Notebook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
Installation and first steps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2. Exploration from the data 13

Preparation of the environment13
Data Import 15
Descriptive statistics16

3. Regression linear by Ordinary Least Squares 21

Simple linear regression ............................... 23
Multiple linear regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

Verification of assumptions 39
Assumptions about the model structure .............................. 40
Preparation of the environment................................................. 40
Data import 40
Number of observations greater than number of parameters41
Variation in the explanatory variables 42
No multicollinearity 43
No specification bias 48
Assumptions about the term of error 52
Mean value equal to zero. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
Homoscedasticity53
No autocorrelation 57
Normality. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

Minimums Squares in 2 Stages 66

Preparation of the environment67
Data import .......................................................... 68
Manual process 71
First stage 71
Second stage74
Automatic process 76
MCO vs. MC2E . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
Comparison of models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6. Variable discrete dependent - Binary case 82
Linear probability model 83
Preparation of the environment.................................................................................83
Data import 84
Data preparation ...
Logit Model 88
Logit Model in Python90
Probit Model 92
Probit Model in Python 93
Comparison of models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

7. Variable discrete dependent - Non-binary case 99

Multinomial Logit Model in Python 103
Environment preparation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
Data import ................................................................................... 103
Data preparation 105
Construction and estimation 105

References 111
Chapter 1

Installation of the programs and preparation of the

environment

Throughout this work, it is intended to develop practical exercises, with illustrative purposes, that serve
as a reference for the use of the Python programming language as a useful tool for economic analysis
three and that allow the reader to become familiar with it and strengthen their skills in handling it.
This work is based on the use of Python 3, with Anaconda distribution and working on Jupyter Notebook;
although these names may seem strange to the reader, especially the novice, it is appropriate to point out that this
you shouldn't worry, since these are concepts without much complexity (for our purposes) and that,
Next, it is intended to clarify, so that you feel complete confidence in the study of this work.

Python
Python is an object-oriented, interpreted, high-level programming language with dynamic typing.
dynamic (Python Software Foundation, n.d.), created by Guido van Rossum, a computer scientist and
Dutch programmer, and which appeared in its initial version in 1991. Python is of very high level, which suggests
that its syntax is relatively easy to understand and this is reflected in the fact that its process of
learning is rather simple compared to that of other programming languages.
It has been noted that we will work with Python 3, however, the reader has not been informed about the difference.
between Python, plain and simple, and Python3. This should not cause you any unease, as the distinction is quite simple:
Python is the programming language, while Python3 is a version of that language.
is available in version Python 2). In this work, Python 3 will be used as it is the version that
there is evidence of a growing trend toward widespread adoption and is characterized by active development.
The use of Python throughout this work is applied and specific, that is, it is employed for purposes
very clearly delineated, performing useful tasks to address the intended econometric analysis
carry out; at no time does this work aim to be a programming manual or anything that resembles it,
therefore the reader should not consider that with the study of this work they will learn to program, as such an objective
exceeds, by far, the scope of this guide. To the reader interested in programming as such, it is suggested
to consult other types of resources that may be useful for this purpose, as this work uses Python
only as a tool and not as the main topic to be addressed.
As has just been pointed out, Python is a tool that we will use to develop analysis.
econometric; however, such an approach may not be entirely clear to the reader, who is natural
that it questions how Python will be used to carry out econometric analysis? The answer to such
the question does not imply greater complexity: it will only make use of the functions, methods, and attributes that
they offer us Python and its specialized libraries (which we will address later) to handle the

1
data, extract the required information, build the models, make the estimates and run the tests
necessary statistics. Although such tasks do involve operations with a certain level of difficulty and imply
elaborate calculations, the user will not have to make an extraordinary effort, in fact, they will not even have to
to define a significant difficulty function on its own or to write 'complex' code, in terms of
Python and its specialized libraries will provide the tools we will need throughout our work.

Anaconda
Although it is possible to install Python using the official distribution from the Python Software Foundation website.
tion, [Link] such a means will not be resorted to, as the Anaconda distribution will be used.
The explanation for such a decision lies in the fact that the Anaconda distribution installs automatically-
specifically multiple libraries (very useful for the tasks that will be undertaken), it includes a package manager
proprio, conda, the standard package management system of Python (pip) and, in addition, gives us access to Anacon-
from Navigator, a graphical interface with which we can access Jupyter Notebook, which is where, ultimately
Finally, we will write the Python code with which we will carry out our econometric analysis.
Anaconda Distribution is a product of the American company Anaconda, Inc., (formerly known as)
known as ContinuumAnalytics) and consists of a free and open-source Python distribution widely
used in the field of data science and machine learning. For more information, it is suggested to consult the
official documentation website [Link] although, for our purposes-
This will not even be necessary, as the required instructions will be given here by the reader.
to correctly follow the development of the topics discussed.

Jupyter Notebook
Jupyter Notebook is an open-source web application that allows you to create and share documents that
contains code, graphs, equations, and text (Project Jupyter, 2019). It supports various programming languages
Among them, Julia, Python, and R have an attractive, easy-to-use, and very user-friendly interface.
user, and is part of Project Jupyter, an initiative that emerged in 2014 from IPython and directed by Fernan-
Pérez, Colombian physicist from the University of Antioquia, assistant professor in the Department of Statistics.
from UC Berkeley and researcher at the Berkeley Institute for Data Science.
The choice of Jupyter Notebook as the application in which our economic analysis will be developed
Trico is based on the pleasantness of its interface and its extremely simple and intuitive handling, which will allow us to
write code, visualize results, create graphs and make comments in the same document, enrich
vastly enhancing our analysis and facilitating the understanding of the content created (in fact, this book
It has been written in Jupyter Notebook, about
A LTXy markdown, with later editing in Overleaf). For further
E
Information about Project Jupyter, you are invited to visit the official project website: [Link]

Installation and first steps

At the beginning of this chapter, it was pointed out that "this work is based on the use of Python 3, with distribution
Anaconda and working on Jupyter Notebook”; the reader has already been given a brief description of each one.
of these terms, however, a specific set of instructions has not yet been provided to be able to
access the 'material' we will be working with. Therefore, in this section, we will proceed to present a
detailed description of the steps to follow to install the required programs and prepare the environment,
so that the main content of this basic econometrics guide can effectively begin with
Python.

2
 The first thing we need to do is have the Anaconda distribution. In case you don't have it installed.
previously (which is very likely, especially for an unfamiliar user on the subject), we must
access the official site [Link] Upon entering such an address, we find ourselves
We start with a page like the following:

We can observe that at the bottom there is a section that says Windows | macOS
Linux. At this point, we must click on the name corresponding to the operating system of our
computer; this, with the purpose of accessing the Anaconda Installer for our operating system. In is-
In this particular case, we will select Windows; however, the reader must be careful, as they must
select the operating system of your computer (which may not be Windows). When selecting the system
corresponding operation, we will be able to see the following:

We observe that there are two installers: one for Python 3 (version 3.71specifically) and another one for Python
2 (version 2.72specifically); the one we must select, as previously noted, is that of
1 The most recent version at the time of publication of this work.
2 The latest version at the time of publication of this work.

3
Python 3. Regarding the specific Graphical Installer, it depends on the processor of our computer: if
it is a 64-bit processor, we will naturally download the 64-Bit Graphical Installer, and if it is about
For a 32-bit processor, our choice will be the 32-Bit Graphical Installer. For the reader who is unfamiliar with the
computer system, this can be checked on Windows in the properties of the computer or at:

There, you will find, among other data, the information corresponding to the system; information similar to what
is presented below:

Once the corresponding Graphical Installer is selected, you just need to save the file and, once downloaded-
To do it completely, it must be executed. This process is very simple (as it just consists of following instructions.
very concise), so it will not be presented step by step, as doing so would be superfluous. Once removed -
loaded and installed the program, we will access Anaconda Navigator; for this purpose, we will resort to a
search from the Windows start button, as indicated below:

4
 We will click on Anaconda Navigator (Desktop application) and wait for it to load.
full (this may take some time, so the user should not get anxious). We will receive
so the following message:

5
 About this message, we just need to click on Oko in Ok, and don’t show again. Once rea-
Once this step is completed, we will have access to Anaconda Navigator, where we will find a home screen similar to this.
to which is presented below:

The reader will remember that we will use Python 3, with Anaconda distribution, working on Jupyter Notebook.
At this point, if the given instructions have been followed, we already have the Anaconda distribution for Python 3.
on our computer; the only thing needed is to have access to Jupyter Notebook, which is, to tell the truth,
the simplest step of all: it just consists of clicking.
The reader can observe that Jupyter Notebook is one of the applications present on the home screen.
from Anaconda Navigator; to access it, simply click on Launch:

6
 Once pressed, we will be automatically redirected to localhost. The redirection process may take a while.
a little, for which the user is advised to be patient; once this process is completed, we will be
Jupyter Notebook, where we will finally be able to start our econometric analysis and get in touch
cash with Python.
The description provided so far about the process to access Jupyter Notebook has been pre-
sitting considering the immediately prior installation of Anaconda; however, once it
With this, it is not necessary to access Jupyter Notebook through Anaconda Navigator. One more way.
brief, and with a faster execution, is the access from Anaconda Prompt; it is enough to perform the search and
click on the corresponding application:

7
 Once clicked, you must type jupyter notebook in Anaconda Prompt and press the Enter key.
(Intro), after which Jupyter Notebook will open and additional information will appear in Anaconda Prompt.
about the executed process.
Finally, the simplest way (once Anaconda is installed) and, of course, the most obvious
to access Jupyter Notebook, you need to search for the application and click on the corresponding icon.
tooth

8
 Any of the 3 access methods mentioned above is completely valid and should allow
to the user to enter Jupyter Notebook; while in this application, you must click on the New tab.
find in the upper right section, and then click on Python 3, as indicated below:

Once we have clicked on Python 3, a new tab will open in the browser, which will have a
appearance like the following:

9
 It is at this point where our real work begins, for it is now that we can start to
write Python code and carry out the econometric analysis that we have planned.
The reader will notice that the appearance of the Jupyter notebook is quite pleasant, with some similarity
a traditional text editor or a common office program. In reality, the feeling of security that
This 'familiar' appearance is effectively reflected in ease of handling: we see an environment
not overloaded and most commands are self-explanatory; to the user not completely unfamiliar with
Those who are knowledgeable in computing will find the use of Jupyter notebooks truly simple. Perhaps the reader has noticed
that the Python logo is located at the top right of the notebook, and that near it is written
Python 3: this indicates that the code written in the notebook corresponds to Python code. Likewise,
we observe that, at the top of the notebook, next to the Jupyter Notebook logo, there is the
Untitled; this is the name that has been assigned by default to our notebook, which, if it
we wish, we can change to the one of our preference by clicking on Untitled, or by resorting to the
Save as options and Rename from the File menu:

We will change the name of our notebook and we will call it, for illustrative purposes, 'My first note-'
book”:

Now, we will proceed to write our first line of code; for this, we will use the prede- function.
findaprint(...). The reader can observe that in the notebook there is a cell without any content; it is here
where we should write our code. For this specific example, the argument of the function print(...)
This is the first line of code of my notebook written in Python.

At this point, the user has written their first line of Python code in the notebook, but has this code...
Has there been any effect? Has any action been taken? The answer is simply no. Our line of code

10
it has not generated any results as we have only written it; we have not executed it yet. To execute
To run the code, we must press the Run button, located at the top, below the Cell tab, or, in another way
alternative and faster, use the key combination Ctrl + Enter.
When executing the code, we will obtain the following:

We can see that the message contained in the code has been printed below our cell and that
in the brackets next to the code cell, which were empty, has appeared the
number one. This number indicates that it is the first execution of code we have performed in the notebook; if
we enter the cell again and run the code again (with the combination Ctrl + Enter),
we will observe that number one is no longer present but number two is (the reader is invited to verify this for themselves)
cuenta).
We have already created a Jupyter notebook, we have assigned it a name, and we have written and executed something.
of code within this. These are some of the basic instructions that the user, without any excuse, must
to know; finally, we will save the changes made in the notebook and proceed to stop and close it.
Knowing this set of operations, the user will have the fundamental preparation required to tackle
successfully the true content in which this work specializes.
The process used to save the changes made in the notebook is very simple and practically
identical to that used with any common office program, so it will not pose any difficulty to
any user with basic computer knowledge. Just use the key combination Ctrl + S,
or click on the Save and Checkpoint tab (identified by the shape of a floppy disk) on the ribbon
options.
Although the notebook can be closed by closing the tab in the browser where it is located, such a process-
the result is not advisable. To close the notebook correctly, you must click on Close and
Open the File menu, as indicated below:

Having completed such a process, the browser tab will automatically close and we will be redirected.
to the Jupyter Notebook home screen, where we should find the notebook we just created,
saved and closed
We can observe that our notebook 'My first notebook' has a size of 814 B and has been saved.
given with extension .ipynb; this is because Jupyter notebooks are files that are stored auto-
Mathematically with that extension. To access the notebook again, simply click on its name;
In this way, a new tab will open in the browser where the notebook will be found.

11
 It is important to emphasize that every time the notebook is opened, a "new session" is being started.
no line of code has been executed and, although it is possible to view the output of the lines of code, the
the output displayed is not the result of executions from the current 'session'. This is important because, if
It intends to execute a particular line of code that requires the prior execution of another specific line and this
the last one has not been executed, an error will occur.
The process to enter Anaconda Navigator and Jupyter Notebook has already been outlined (mentioning 3
alternatives), now it is time to address the reverse process: how to exit. Although it is possible to exit Jupyter
Notebook with just closing the browser in which it is located, such conduct is inadequate. The method
To properly exit Jupyter Notebook, you need to close the active notebooks and then close it.
Jupyter Notebook by clicking on the tab located in the upper right section:

Once pressed, we will receive the following message:

Now we can close the browser with complete peace of mind. As for Anaconda Navigator, after
If the process just described has been completed, if it is open, simply close it as if it were a program.
Any automatic, just by pressing the tab with laxen at the top right.
To access Jupyter Notebook again, it is only necessary to repeat the necessary instructions that are
have been indicated throughout this chapter, which are actually just accessing Anaconda Navigator and,
from this application, enter Jupyter Notebook (or simply access from Anaconda Prompt or Jupyter
Notebook directly) being redirected to localhost.
With this chapter, the discussion of the programs to be used comes to a close; it is now time to delve into the true
subject that addresses this work, so from the following chapter, the purpose of this will be developed
practical guide: econometric analysis with Python.

12
Chapter 2

Data exploration

It is impossible to carry out an analysis without having data to analyze; this is nothing short of absurd. Thus,
if we want to conduct an econometric analysis exercise, we must have the required information to
to carry it out; the reader should not worry about this comment, as, in addition to the practical guide
what is offered to him/her, information is also provided with which the models are built and structured
analysis.
For each of the exercises that will be developed in this work, the source from which it has been obtained will be indicated.
Having the information used in the corresponding econometric analysis. Likewise, it is possible to obtain the
datasets used, in the same formats that are employed in this work, by contacting the author of it
(just send an email to the address fatrianaa@[Link]), who will give you access to the requested information.
Rida, or alternatively, access the website of the Computer and Communications Unit of the Faculty of
Economic Sciences of the National University of Colombia.
This work, given its practical nature and the utility it can signify for undergraduate students
The degree in Economics mainly relies on Gujarati and Porter's Econometrics (2010) and Introduction to the
Econometrics. A modern approach by Wooldridge (2010), very simple texts widely used such as
guide in the basic courses of Econometrics, which this work refers to multiple times,
to develop the illustrative exercises posed and with which it is intended that the reader can compare
results, in such a way that it shows that the procedures carried out are correct and lead to the solutions
expected.
The dataset that will be used for simple linear regression by Ordinary Least Squares that will be addressed.
In the next chapter is the employee in Example 7.1, Infant mortality in relation to GDP per capita and
the literacy rate of women, by Gujarati and Porter (2010). The reader is invited to consult the example, of
so that I can verify the results obtained.
The source of the information used by Gujarati and Porter (2010) is “Chandan Mukherjee, Howard White and
Marc Whyte, Econometrics and Data Analysis for Developing Countries, Routledge, London, 1998, p. 456.
data used for the exercise in Python that takes place in this work were obtained from the data file of
Gujarati is a language spoken in the Indian state of Gujarat.

At this point, our work with Python begins. The first thing the user must do is access
Jupyter Notebook (through Anaconda Navigator or by your preferred method) and create a new notebook.
In this notebook is where the required code will be written and the proposed econometric analysis will be carried out.

Environment preparation
The first cell in which the user will write and execute code will have the following content:

13
In [1]: We import the required libraries:
import numpy as np
import pandas as pd
import [Link] as sm
import [Link] as plt
importseabornassns
%matplotlibinline
[Link]("seaborn-white")

What are these lines of code about? The reader will notice that the first thing that has been written, as a co-
We import the required libraries; this is what, effectively-
actively, is achieved by executing this code cell. Basically, what we achieve is to prepare our
work environment when gaining access to the toolkit that we will need for our analysis; such
tools are part of libraries.
A library is simply a set of functions and methods designed to perform specific tasks.
cases, grouped in a single space (the library), and which is added to basic Python with the aim of achieving
results that cannot be easily achieved with their integrated (native) functions. To the familiar reader
With the programming language R, it suffices to say that a Python library is essentially equivalent to a
R package.
The libraries we import are NumPy, pandas, Statsmodels, Matplotlib, and Seaborn; each of these has
un rol particular a desempeñar en nuestro trabajo, el cual se presentará, junto a breves descripciones, en un
moment. One thing the reader should notice is that importing libraries is not just about 'importing
the libraries", without specifying any type of additional instruction: the NumPy library is simply imported
how to import numpy but it is imported as import numpy as np. What is the reason for this? What effect
generate?
The import process with as, which is a reserved word in Python, is carried out for the purpose
to assign an alias to the imported libraries. The reader might be wondering: if the library already has a
name, why is it wanted to assign a nickname? The answer to the question is very simple: just for
practicality. When importing the NumPy library with the alias np, that library is identified by this name
Therefore, every time you need to use a function or method from NumPy, it is enough to refer to it.
shy not shy.
The utility may be more apparent in the case of the Matplotlib library: instead of referencing
In [Link], whenever a function or method from it is required, it is sufficient to
refer to aplt. Thus, instead of using 17 characters, only 3 are required; although to the reader this
It may seem trivial at this moment, but over time you will see that the importation will turn out to be tremendously practical.
and greatly facilitates the work.
The imported libraries are NumPy, pandas, Statsmodels, Matplotlib, and Seaborn. Why is this group imported?
specific to libraries and not others? To understand this, a brief description of each of these is offered,
so that the reader has clarity about their purpose.

NumPy is a Python library for scientific computing that provides a set of routines for the
execution of procedures on objects, such as vectors and matrices, which include, among others,
mathematical and logical operations, basic linear algebra, random simulation, and fundamental statistics
(The SciPy community, 2019).
Pandas is an open-source library that offers high-performance and easy-to-use tools.
for data management and analysis in Python. It is a library that allows for the execution of a process
optimal information analysis directly in Python, without having to resort to a more advanced language
specificity, such as [Link] used in both academic and commercial fields, in areas

14
 various that include finance, neuroscience, statistics, and economics (Pandas, n.d.). The creator
Benevolent Dictator for Life depends on the mathematician and American developer Wes McKinney.

Statsmodels is a library for econometric and statistical analysis in Python that provides functions for
the estimation of multiple statistical models, as well as the conduct of tests and analysis of information
statistical information. The results generated by the library functions are contrasted with the
of recognized statistical packages, so that their accuracy can be guaranteed. The library has
its origin in the module models from [Link] in 2009, after a process of correction, testing and
Improvement, Statsmodels was launched independently (Perktold, Seabold, and Taylor, 2018).

Matplotlib is a 2D graphics library used in Python for generating high-quality images.

Dad. With Matplotlib it is possible to generate histograms, bar graphs, scatter plots, among others.
others, by using just a few lines of code. The pyplot module provides an interface
similar to MATLAB, useful for simple graphics, allowing complete control over various properties
Matplotlib is a contribution from John Hunter, an American neurobiologist, and his multiple
collaborators (The Matplotlib development team, 2018).

Seaborn is a library for creating statistical graphics in Python. It is based on Matplotlib.

and its graph-building functions are oriented towards datasets, so it is closely
grade with the data structures of pandas. The purpose of Seaborn is to constitute visualization in a
key element in the exploration and understanding of information (Waskom, 2018).

The reader, considering the brief descriptions just presented, should already have an idea
of the role played by each of the libraries mentioned in the analysis to be carried out; without
embargo, if this is not the case and the previously exposed information seems a bit confusing to you, in words
the utility of each of the libraries is as follows:

NumPy: Numerical Analysis.

pandas: Data management.
Statsmodels: Construction of models and making estimates.
Matplotlib: Creation of graphs (general).
Seaborn: Additional functions for creating specialized graphs.

Once the lines of code from the first cell are executed, the required libraries will have been imported and
we will have at our disposal the set of tools that we will use to carry out the analysis
econometric. Now, we will proceed to load the information we will work with.

Data import
The import of the data is carried out using functions from the pandas library. Depending on the type of
file in which the dataset information is stored, we will need to resort to a partitioning function
cellular and a set of specific parameters.
The dataset we will use corresponds to a .txt file, so we will resort to the function
read_csv(...) dependencies will indicate the values corresponding to the delimiter parameter.

Note: The required parameter of the function read_csv(...) is the filepath_or_buffer.

file path, that is, the location of the file containing the dataset, expressed as a string.

For our specific case, the following code is the one that imports the dataset:

15
In [2]:data=pd.read_csv("C:/Users/FCE/Documents/Econometrics/[Link]",
sep=" " , delimiter="\t")

The reader should note that it is not pd.read_csv(...) but read_csv(...). This is because
read_csv(...) is not a 'native' function of Python but belongs to pandas, which is why it is necessary-
Specify where it comes from. Since the import was made with import pandas as pd, the code
pd.read_csv(...) informs that the function read_csv(...) belongs to pandas. If the import is not h-
it should have been carried out using an alias, pandas.read_csv(...) should have been used for this purpose.
follow the result.

Note: The user must take into account that the value assigned to the parameter filepath_or_buffer will be the
the same as the code in the previous example, as it depends on the location of the file
[Link] (or the file in which you have saved the dataset) on your computer.

Once the line of code in this cell is executed, the dataset should have been imported correctly.
To check it, the .head() method is used on the object that contains the dataset (in this case, it
has assigned the name data):

In [3]:[Link]()

Out[3]: CM FLR PGNP TFR

0 128 37 1870 6.66
1 204 22 130 6.15
2 202 16 310 7.00
3 197 65 570 6.25
4 96 76 2050 3.81

We can observe that the import process has been successful and the dataset information has been soul-
He correctly seated in a Pandas DataFrame. Now, we can continue calmly with the development.
from our analysis, as we have successfully completed the first step.

Descriptive statistics
Before starting to examine the values of the dataset variables and identify possible relationships that
it's advisable to conduct a recognition of one's own dataset, that is, to obtain information
about its dimensions, the type of data it contains, etc. To know the dimensions of the dataset you can
refer to the attribute .shape:

In [4]: [Link]

(64, 4)

The result of applying the [Link] is a Python tuple (...) that reports the number of
rows and columns (in that order) that the DataFrame has; thus, our dataset contains 64 rows and 4 columns.
Now, we know there are 4 columns, but what are they called? What is the name of each of these 4 co-
columns? When the method .head() has been applied, the first observations of the dataset have been visualized.
including the header where the column names are located; however, if specifically
If you want to know the names of the columns, you can use the attribute .columns:

16
In [5]: [Link]

["CM","FLR","PGNP","TFR"]

With all the obtained information, it is known that the dataset to be worked on has 64 rows and 4 columns.
and that the column names are 'CM', 'FLR', 'PGNP', and 'TFR'. However, it is still unknown what type of
the data contain those columns, are they words? are they numbers? only whole numbers? To find out,
the attribute can be used: .dtypes

In [6]:[Link]

Out[6]:CM int64
FLR int64
PGNP int64
TFR float64
dtype: object

The information generated by the [Link] indicates that three variables (columns) contain data from
tipoint64 (integer number) and one of type float64 (decimal number).
Finally, a way to obtain the information that is found with the application of the three attributes just
described is the [Link](), which allows them to be visualized in a single space:

In [7]: [Link]()

<class '[Link]'>
RangeIndex: 64 entries, 0 to 63
Data columns (total 4 columns):
CM 64 non-null int64
FLR 64 non-null int64
PGNP 64 non-null int64
TFR 64 non-null float64
dtypes: float64(1), int64(3)
memory usage: 2.1 KB

El resultado de la aplicación de este método es tan solo un resumen de la información que hemos obtenido
previously. Additionally, the data on the associated memory usage is presented, which in this case is 2.1
Kb; the type of object to which the dataset corresponds, which is a pandas DataFrame; and the range of the index,
which corresponds to the interval[0, 63] It is very important to keep in mind that indexing in Python starts
from 0 and from 1).
To obtain descriptive statistics for the numerical variables of the dataset, one can resort to the method
.describe():

In [8]:[Link]()

Out[8]: CM FLR PGNP TFR

count 64.000000 64.000000 64.000000 64.000000
mean 141.500000 51.187500 1401.25 5.549687
std 75.978067 26.007859 2725.695775 1.508993

17
 min 12.000000 9.000000 120.000000 1.690000
25 % 82.000000 29.000000 300.000000 4.607500
50% 138.500000 48.000000 620.000000 6.040000
75% 192.500000 77.250000 1317.500000 6.615000
max 312.000000 95.000000 19830.000000 8.49

As we can see, the method .describe() allows us to obtain basic statistical information about
the numerical variables, which include, for example, the average, the maximum and minimum values, and the standard deviation
standard. Thus, it can be evidenced, for example, that the minimum value of the variable 'CM' is 12, the maximum value
de la variable ’PGNP’ es 19,830 y el promedio de la variable ’TFR’ corresponde a 5.55.
By applying the method .describe() on the DataFrame, statistical information is obtained by default.
basic of all the numerical variables present in this. If one is only interested in the information of one
a specific variable, this can be selected using square brackets [ ], indicating its name
(in quotes) within these.
Thus, if we are only interested in the 'CM' variable, the code to be used will be the following:

In [9]:data["CM"].describe()

count 64.000000
mean 141.500000
std 75.978067
min 12.000000
25% 82.000000
50% 138.500000
75 % 192.500000
max 312.000000
Name: CM, dtype: float64

Yes, on the contrary, if what interests us is a group of variables and not just one, we can select it.
using the brackets [ ] and including the names of the variables of interest (in quotation marks, separated
by commas) within a Python list [...]:

In [10]:data[["CM","PGNP","FLR"]].describe()

Out[10]: CM PGNP FLR

count 64.000000 64.000000 64.000000
mean 141.500000 1401.250000 51.187500
std 75.978067 2725.695775 26.007859
min 12.000000 120.000000 9.000000
25% 82.000000 300.000000 29.000000
50% 138.500000 620.000000 48.000000
75% 192.500000 1317.500000 77.250000
max 312.000000 19830.000000 95.000000

By this point, it is known that the dataset contains information about 4 numerical variables (3 of integer values.
and 1 of non-integer values), counting with data of 641individuals (entities) for each of these. It has already
1 Inthis case, the information of each individual (entity) corresponds to a single row, so, since the dataset has 64 rows,
we know that there is information on 64 individuals (entities). However, in other cases, such as data panels, the information
the same individual (entity) does not correspond to a single row, so one must be careful in such situations.

18
A very brief review of the data has been carried out, which allows for a greater level of familiarity with
the same and it is possible to identify potential relationships to examine.
So far, a data recognition process has been carried out through which
it has obtained purely numerical information: the number of rows and columns, the mean and the standard deviation
of each variable, etc. This gives us some idea about the data; however, often, 'a picture is worth a thousand words.'
a thousand words2 it is possible for a graphic to have immense communicative power and contribute enormously
to gain a better understanding of the data that will be worked with.
The pairplot(...) function from Seaborn is particularly useful for the context in which we find ourselves.
This function creates a grid where we can visualize the relationship that exists between pairs of variables.
(through scatter plots) and the distribution of each of these (through histograms):

In [11]: g=[Link](data, plot_kws={"color":"darkblue"},

diag_kws={"color":"darkblue"})
Relationship between pairs of variables
fontsize= 25,
fontweight="bold"
plt.subplots_adjust(top=0.9)
[Link](1,-.02,
Preparation:
fontsize= 13, fontweight="bold"
ha="right")
[Link](1, -0.05,
Triana, F.\n(2019)
fontsize= 12, ha="right"
[Link]()
2 Oh, instead of words, rather numbers in this case.

19
 In the diagonal, the distribution of each of the variables can be observed, while in the diagrams of
dispersion finds the relationship between pairs of these; thus, it can be observed that, for example, income
it focuses on values below 5000 and there is a negative relationship between the literacy rate of the
women and infant mortality.
At this time, we already have basic information about the variables of the dataset and we have made progress.
a brief visual inspection of the relationships presented among these; now, this can be used
knowledge to obtain more concrete expressions of such relationships, examining them from a point of
A somewhat more technical view: the next chapter will address the topic of linear regression by the method of
Ordinary Least Squares.

20
Chapter 3

Linear Regression by Least Squares

Ordinaries

Linear regression is the introductory topic usually covered in Econometrics courses.

undergraduate degree in Economics; according to Greene (2003), 'the linear regression model is the most useful tool of
"toolkit of econometricians" (p. 7). This topic is often the first contact that the student has.
with matter of purely econometric content and tends to be approached in some detail, such that
this has enough clarity about it. It is common for courses to present linear regression, along with
to the applicable theoretical foundations, supported by a guide text from a recognized author and
that there is noticeable emphasis on the associated assumptions, so that the student acquires a level of
adequate knowledge about this.
In this work, considering its practical nature and taking into account the warning made in the preface
that the theoretical foundations will be addressed with a minimum level of depth, it will be limited to exposing
the general idea of linear regression in very superficial terms; this is with the purpose of focusing on the
stated objective and not to make the reader from other fields of study or with a
not very solid econometric training, as it does not delve into the technical foundations and mathematical formalities on
the foundations of linear regression and its estimation methods.
Regression analysis deals with the study of the dependence of one variable (dependent variable)
regarding one or more variables (explanatory variables) with the aim of estimating or predicting the
average the population average of the first in terms of the known or fixed values (in
repeated samples) of the second. (Gujarati and Porter, 2010, p.15)
In simple terms, regression analysis studies the relationship that exists between the dependent variable and
the explanatory variable(s). For the case of linear regression, such relationship can be expressed as a
linear function in the parameters, but what parameters are being discussed?
To solve the previous question, the first thing to do is to consider a bivariate regression, in the
which is the dependent variable (identified as y) can be expressed as a function1(with linearity
in the parameters) of the explanatory variable (identified with x). Specifically:

y i= β 0 + β 1 x i+ u I
This expression corresponds to the Population Regression Function (PRF) (Gujarati and Porter, 2010) and the
the parameters that have been referred to are simply theβas we can see, theseβ they cannot be seen
affected by some type of transformation that modifies their linear relationship with the dependent variable: to this
1 It is not exactly a function, since, as Fahrmeier, Kneib, and Lang (2007) point out, a fundamental characteristic
from regression analysis is that the relationship between the dependent variable and the explanatory variables does not correspond to a function in
strict sense, as it is affected by random disturbances (p. 19).

21
it refers to collinearity in the parameters, which does not necessarily have to apply to the variables
explanatory (Greene, 2003; Gujarati and Porter, 2010). Regarding this, it corresponds to the term of error, which
It captures the effect of the variables not explicitly included in the model and that affect it.
Now, in practice, the Population Regression Function is unknown, so we resort to the Function
of Sample Regression (FRM) as an approximation to this. To the reader interested in the issue and who wishes
to understand an explanation regarding this, it is suggested to consult Chapter 2, Regression Analysis with two variables:
Some basic ideas from 'Gujarati, D.N. and Porter, D.C. (2010). Econometrics'.
The Sample Regression Function is:

yi= β 0ˆ+ β 1 xˆi+ uI ˆ

This expression is very similar to that of the Population Regression Function, but it has a difference.
fundamental: the terms of the Sample Regression Function have a direct correspondence with the terms
included in the expression of the Population Regression Function, but they are not exactly the same, in
As far as the estimators are concerned. Thus, the estimator of the parameter k is identified as k. ˆ
As noted by Gujarati and Porter (2010), "the main objective of regression analysis is to estimate the
Population Regression Function based on the Sample Regression Function" (p. 44). In this sense, it
what is sought is to find the values of theβˆ that are closer to the true values of theβ.
To carry out the estimation, the use of two methods is common: Ordinary Least Squares and Maximum
Likelihood. The Ordinary Least Squares method, or OLS, is the most common and tends to be addressed
as a fundamental theme in introductory econometrics courses.
Highlighting, once again, the level of depth with which the theoretical foundations are addressed in this work.
it is minimal, it will only be presented (with the specialists' pardon, not addressing the issue with the right detail)
the general idea of the OLS method: the objective is, basically, to minimize the sum of the residuals (term) ˆ
in the Sample Regression Function) squares.
For the case of a bivariate regression, what is sought with the OLS method is the lowest possible value of
n
ˆ 2i considering that:
∑ i=1you,
n n
∑ 2
∑ ˆ ˆ 2
uˆ=
i ( y i− β 0 − β 1 x i )
i =1 i =1

Thus, the estimatorsβˆ0y βˆ1with which the smallest sum of the squared residuals is found are the estimates.
MCO pains, which, upon fulfilling a set of specific assumptions, demonstrate properties
very attractive statistics.
ˆ β and become familiar
To the reader interested in understanding the process by which the values are obtained from the
with the corresponding mathematical bases, you are invited to consult Chapter 3, Regression model with two
variables: problema de estimación, de “Gujarati, D.N. y Porter, D.C. (2010).Econometría” y el Capítulo 2,El modelo
of simple regression, from 'Wooldridge, J.M. (2010). Introduction to econometrics. A modern approach.'
A very brief general idea of linear regression and the OLS method has already been presented. Now, it is
It's time to effectively start the practical work; namely, econometric analysis using Python.
The first step will be to carry out an estimation, using Ordinary Least Squares, of a model.
of simple linear regression. The simple linear regression model is nothing more than a bivariate linear regression,
that is, in which only two explicit variables are involved: the dependent variable as a linear function of the
parameters of the explanatory variable. That is, simply, a model to which a FRM corresponds in the form
yI= β 0ˆ+ β 1 xˆI+ ui which
ˆ should be familiar to the reader, as it has been addressed previously.

22
 Simple linear regression
To carry out simple linear regression, we will use the same dataset as in the previous chapter and require
we will laugh at the same tools used in this, so we will proceed to import the respective ones
libraries and the data:

In [1]: import numpy as np

import pandas as pd
[Link]
import [Link] as plt
importseabornassns
%matplotlibinline
[Link]("seaborn-white")

In [2]: data=pd.read_csv("[Link]",
sep=" ", delimiter="\t")

In the previous chapter, a brief overview of the data was conducted. Now, we will proceed to the actual realization.
Estimation of a Simple Linear Regression using the Ordinary Least Squares method, taking into account
ration the following:

The DataFrame contains information about the following variables:

'CM'. This variable refers to Child Mortality and corresponds to infant mortality. It is about
the number of deaths, in a year, of children under 5 years old per 1000 live births
alive.
‘FLR’. This variable refers to Female Literacy Rate and corresponds to the literacy rate of
the women.
'PGNP' refers to Per Capita Gross National Product and corresponds to GDP per capita in 1980.
'TFR' refers to Total Fertility Rate and corresponds to the total fertility rate.

In the simple linear regression that will be presented, the dependent variable will be 'CM' (infant mortality) and
the explanatory variable will be 'PGNP' (GDP per capita).

A visual examination will be conducted to assess the relationship between the variables of interest. To do this, we
It will build a scatter plot using the method .scatter(), applied to an Axes.
Matplotlib. The code to be used is the following:

In [3]: fig, ax = [Link](figsize=(10,5))

GDP per capita vs. Infant mortality
fontweight="bold"
[Link](x=data["PGNP"], y=data["CM"], s= 50, color="darkblue")
ax.set_xlabel('GDP per capita', fontsize=15)
ax.set_ylabel("Child Mortality")\n(per 1000 live births)
fontsize= 15)
plt.subplots_adjust(top=0.9)
plt.tick_params(labelsize= 15)
[Link](.9,-.02,
Preparation:

23
 fontsize= 13, fontweight="bold",
ha="right")
[Link](.9,-.08,
Triana, F.\n(2019)
fontsize= 12, ha="right"
[Link]()

Although it is a somewhat lengthy block of code, the reader should focus their attention only on the first,
especially the third, and the last line of code in the cell, as they contain the 'essence' of the
graph; the others only refer to details that can be dispensed with. There will be no emphasis.
in the explanation of this code, as it is not the main issue being addressed and most lines 'se
"explain themselves." To the reader who is, however, interested in this matter, it is recommended to study
detailed of Matplotlib.
From the graph, it can be observed that the relationship between the variables is apparently negative. Without
embargo, this is just an impression that is generated from observation and can be erroneous, reason for
which will quantitatively examine this relationship seeking to verify if the conjecture is correct; for such
The purpose will be to carry out a linear regression using the OLS method.
To create the model, the OLS(...) function from Statsmodels is used. The arguments that are included co-
correspond, respectively, to the dependent variable and the explanatory variable. The result will be assigned to
an object that will be called MySimpleModel:

In [4]: MySimpleModel = [Link](data['CM'], data['PGNP'])

The model has already been created, that is, its structure has already been defined; however, it has not yet been implemented.
no estimation will be made from this. To carry out the estimation, the [Link]() will be used.
they will store the results generated in an object that will be called SimpleResults:

24
In [5]: SimpleResults=[Link]()
When executing this line of code, no visible result is obtained. This is because what has been
achieved with the execution is to store the results generated in an object, without indicating the execution of
any particular action with said object. Now, if you want to visualize the results of the model, you can
simply use the integradaprint(...) function, using the application of the method as an argument
.summary() about the object that stores the results of that model. Thus:
In [6]:print([Link]())
OLS Regression Results
==============================================================================
Dependent Variable: CM R-squared: 0.056
Model: OLS Adj. R-squared: 0.041
Method: Least Squares F-statistic: 3.710
Date: xxx, xx xxx xxxx Prob (F-statistic): 0.0586
Time: xx:xx:xx Log-Likelihood: -413.92
No. Observations: 64 AIC: 829.8
Df Residuals: 63 BIC: 832.0
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef standard error t P>|t| 0.025 0.975
------------------------------------------------------------------------------
PGNP 0.0124 0.006 1.926 0.059 -0.000 0.025
==============================================================================
Omnibus: 7.668 Durbin-Watson: 0.755
Prob(Omnibus): 0.022 Jarque-Bera (JB): 7.941
Skew: -0.561 Prob(JB): 0.0189
Kurtosis: 4.312 Cond. No. 1.00
==============================================================================

Warnings:
Standard Errors assume that the covariance matrix of the errors is correctly
specified.

What is obtained from the execution of the code is a summary of the results instance in which it
presents the most relevant information about these: at the top, the estimation method is included
employee, the number of observations, the degrees of freedom (of the residuals and the model) and the R2(standard
and adjusted), among other relevant data. In the middle section, the coefficients can be found, along with their
respective standard errors, statistics confidence intervals. In the last section are the warnings,
which are often very useful and can help identify potential problems.
It can be observed that the coefficient of the explanatory variable is positive, which contradicts the conclusion.
It is visually evident that the relationship, apparently, was negative; something that is also somewhat illogical, as one would expect
that countries with the highest per capita income had a lower infant mortality rate. In light of this situation-
Strange action is worth questioning its cause: how can such inconsistency be explained? The answer
it is extremely simple: the reader must notice that for this model only the coefficient has been calculated
corresponding to the variable 'PGNP' and the existence of the intercept parameter has not been considered, that is, it
has conducted a regression through the origin.

25
 Instead of consideringi= β 0+ β 1 xi+ ui it has been assumed that the FRP is of the formi= β 1 xi+ ui .
The OLS(...) function from Statsmodels does not assume by default that the model to be estimated includes a constant.
as an explanatory variable and, therefore, does not estimate the coefficient corresponding to the parameter of
intercept. In order to obtain the estimate of that parameter, it is necessary to specify that a constant must
add to the set of explanatory variables of the model; this is achieved, for example, by using the function
add_constant(...) from Statsmodels, taking as argument the corresponding object of the independent variables
applications included. Thus:

In [7]: MySimpleModel2 = [Link](data['CM'], sm.add_constant(data['PGNP']))

To make the estimation and visualize the corresponding results, it is enough to repeat the same structure.
update of the codes previously used (in the case of the model without an intercept term). Thus:

In [8]: SimpleResults2=[Link]()
print([Link]())

OLS Regression Results

==============================================================================
Dep. Variable: CM R-squared: 0.166
Model: OLS Adj. R-squared: 0.153
Method: Least Squares F-statistic: 12.36
Date: xxx, xx xxx xxxx Prob (F-statistic): 0.000826
Time: xx:xx:xx Log-Likelihood: -361.64
No. Observations: 64 AIC: 727.3
Df Residuals: 62 BIC: 731.6
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef standard error t P>|t| 0.025 0.975
------------------------------------------------------------------------------
const 157.4244 9.846 15.989 0.000 137.743 177.105
PGNP -0.0114 0.003 -3.516 0.001 -0.018 -0.005
==============================================================================
Omnibus: 3.321 Durbin-Watson: 1.931
Prob(Omnibus): 0.190 Jarque-Bera (JB): 2.545
Skew: 0.345 Prob(JB): 0.280
Kurtosis: 2.309 Cond. No. 3430
==============================================================================

Warnings:
Standard Errors assume that the covariance matrix of the errors is correctly
specified.
The condition number is large, 3.43e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

Now, it can be observed that the coefficient corresponding to the explanatory variable is negative, just as
It had been presumed in the visual examination, and this is statistically significant. What will be done next
to graph the original values and the values estimated by the model, so that they can be contrasted and

26
A visual examination of the goodness of fit of this is possible. For this, the method will be used again.
.scatter(...) applied on an Axes of Matplotlib.
In [9]:fig, ax=[Link](figsize=(10,5))
Regression lineal simple (OLS)
fontweight="bold"
[Link](data["PGNP"], data["CM"], s= 50,
label="Valores originales",
color="darkblue")
[Link](data["PGNP"], [Link](), s= 50,
label="Valores estimados",
color="red")
GDP per capita
ax.set_ylabel("Infant mortality")\n(per 1000 live births)
fontsize= 15)
[Link](frameon=True, loc="lower left")
plt.subplots_adjust(top=0.9)
plt.tick_params(labelsize= 15)
[Link](.9,-.02,
Preparation:
font size= 13, font weight="bold"
ha="right")
[Link](.9,-.08,
Triana, F.\n(2019)
fontsize= 12, ha="right"
[Link]()

27
 It is observed that the goodness of fit of the model is not the best. What is the reason for this discouraging result?
If the reader examines the previous graph, they will notice that the part corresponding to the original data presents
a very slight resemblance to the graph of a function of the form y= f ( x ) = , cona>xa 0,x> 0y
representation of the value dex on the x-axis and the value dey on the y-axis. This resemblance
suggests that the functional formi= β 0+ β 1 xi+ uimaybe it is not the most suitable for modeling a relationship
a
With this data, how much can it resemble a function of the form f?( x ) = (a= 0) that a x
1
one of the forms( x ) = ax( a= 0)Considering this, perhaps the wayi= β 0+ β 1 x + uigenerate
better results.
Note: A model with a functional shape.i= β 0+ β 1 1
+ uIis called a reciprocal model.
x
I invite the reader to consult Chapter 6, Extensions of the linear regression model with two variables, by Gujarati.
and Porter (2010), for more information on this.
To work with a reciprocal model, it should be taken into account that the estimation will not be made from the
1
values of the values of. Evidently, this implies that the model is nonlinear in the variable
x
x, but does this mean that what will be carried out is a nonlinear regression? The answer is no. Although the
the model is no longer linear in the variable x, the reader must remember that the linearity considered refers to
the parameters and not the variables, therefore, despite the non-linearity of the variable, it continues to be a
linear regression model, regarding the linearity in the parameters has not been affected.
The model has been created using 'PGNP' as an explanatory variable; now, it will not be included, but rather it will
will resort to its reciprocal: '1/PGNP'. The values of 'PGNP' are included in the dataset that has been imported, without
embargo, the values of '1/PGNP' are unknown. How can a variable whose values are not used?
Do you have information?
Although the values of '1/PGNP' are unknown, there is still the possibility of obtaining both the values.
The 'PGNP' are well known and only a simple mathematical operation needs to be performed. The calculation of these
values, despite being simple, can be tedious given the amount of them and can lead to making mistakes
errors if executed 'manually'; fortunately, pandas is very useful in this context, as it avoids the
carrying out such a 'manual' task for the calculation of each of the values, by allowing operations to be applied
basic Python mathematics that allow for easy results.
The variable corresponding to the reciprocal of the original variable will be created, since it is not part of the
dataset. To create this variable, you simply assign it to a pandas Series, indicating
the name of the corresponding DataFrame and then specifying the name you want to give to the
new variable (written in quotes) inside brackets [ ]. Thus:
In [10]:data["1/PGNP"]= 1/data["PGNP"]
The correct execution of the process is verified by using the method .head():

In [11]:[Link]()
Out[11]: CM FLR PGNP TFR 1/PGNP
0 128 37 1870 6.66 0.000535
1 204 22 130 6.15 0.007692
2 202 16 310 7.00 0.003226
3 197 65 570 6.25 0.001754
4 96 76 2050 3.81 0.000488
As the reader can see, the last column is named after the new variable that has been assigned.
("1/PGNP") whose values are the results of taking, row by row, the number one and dividing it by the corresponding value.
dependent on the column 'PGNP'. To verify the accuracy of the operation performed, feel free
to take the information from any row and execute the corresponding calculation.

28
 Now, we will estimate the model andi= β 0+ β 1 1
xi + ui taking into account that1xcorresponds to the variable
newly created '1/PGNP', and we will observe if an improvement in the fit is generated compared to the one obtained with the estimate.
decision of youi= β 0 + β 1 xi+ uI For the construction and estimation of the model and visualization of the results,
we will use again the OLS(...) function from Statsmodels and the methods .fit() and .summary():

In [12]: MySimpleModel3 = [Link](data["CM"], sm.add_constant(data[["1/PGNP"]]))

ResultadosSimple3=[Link]()
print([Link]())

OLS Regression Results

==============================================================================
Dep. Variable: CM R-squared: 0.459
Model: OLS Adj. R-squared: 0.450
Method: Least Squares F-statistic: 52.61
Date: xxx, xx xxx xxxx Prob (F-statistic): 7.82e-10
Time: xx:xx:xx Log-Likelihood: -347.79
No. Observations: 64 AIC: 699.6
Df Residuals: 62 BIC: 703.9
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef standard error t P>|t| 0.025 0.975
------------------------------------------------------------------------------
const 81.7944 10.832 7.551 0.000 60.141 103.447
1/PGNP 27270 3759.999 7.254 0.000 19800 3.48e+04
==============================================================================
Omnibus: 0.147 Durbin-Watson: 1.959
Prob(Omnibus): 0.929 Jarque-Bera (JB): 0.334
Skew: 0.065 Prob(JB): 0.846
Kurtosis: 2.671 Cond. No. 534.
==============================================================================

Warnings:
Standard Errors assume that the covariance matrix of the errors is correctly
specified.

Now, we proceed to graph the original values and the estimated values, using the method
.scatter(...) of a Matplotlib Axes, so that we can observe if it is indeed generated a
improvement in the fit:

In [13]:fig, ax=[Link](figsize=(10,5))
Simple linear regression (OLS)
fontweight="bold")
[Link](data['PGNP'], data['CM'], s=50,
label="Valores originales",
color="darkblue"
[Link](data['PGNP'], [Link](), s=50,

29
 label="Valores estimados",
color="red"
ax.set_xlabel('GDP per capita', fontsize= 15)
ax.set_ylabel("Infant mortality")\n(per 1000 live births)
fontsize= 15)
[Link](frameon=True, loc="lower right"
plt.subplots_adjust(top=0.9)
plt.tick_params(labelsize= 15)
[Link](.9,-.02,
"Elaboración:",
fontsize= 13, fontweight="bold"
ha="right")
[Link](.9,-.08,
Triana, F.\n(2019)
fontsize= 12, ha="right"
[Link]()

As we can observe, there is a clear improvement, as the adjustment seems to be more suitable; the me-
Jora is visually evident; however, there has also been a noticeable increase in the coefficient.
of determination, or2 , which has increased from 0.166 to 0.459, as the reader can verify in the tables of
summary of each regression. Thus, for this particular case, the functional form of the reciprocal model seems to be
more appropriate than the traditional functional form.

30
 Multiple linear regression
We have already performed a linear regression using a single explanatory variable; now we will use several
explanatory. In a multiple linear regression model, the FRP is of the form andi= β 0+ β 1 x1I+
... + β k xwhat+ ui so that explanatory variables are included (excluding the constant ofβ 0 ) and must be
estimark+ 1 parameters (those of the variables and the intercept term).
In the multiple linear regression model that we will use, the dependent variable will be 'CM' (mortality
infantile) and the explanatory variables will be 'PGNP' (GDP per capita) and 'FLR' (female literacy rate)
The first thing we will do is graphically examine the relationship that exists between the explanatory variables.
and the dependent variable; this time, we will resort to the pairplot(...) function of Seaborn, as it allows us to
it allows graphing pair relationships in data sets.
The function pairplot(...) of Seaborn includes multiple arguments; however, for our purposes...
In addition to (obviously) the mandatory parameter data, the important arguments are x_vars and y_vars;
first allow specifying the variables to be graphed on the x-axis and the second ones
wants to graph on the y-axis. For our particular case, y_vars = 'CM' and x_vars = ['PGNP',
The reader should note that, since there are several variables, they must be included in a list.
hon [...]).

In [14]: g = [Link](data, x_vars=["PGNP", "FLR"],

y_vars="CM", height= 5,
plot_kws={"color":"darkblue","s":50})
[Link]('Correlation between variables', fontsize=18)
fontweight="bold", y= 1.08
[Link][0].set_xlabel('PGNP')\n(GDP per capita)
CM\n(Infant mortality)
FLR\n(Women's literacy rate)
fontsize= 15)
[Link]:
ax.tick_params(labelsize= 15)
[Link](1,-.09,
Preparation:
fontsize= 13, fontweight="bold"
ha="right")
[Link](1, -0.15,
Triana, F.\n(2019)
fontsize= 12, ha="right"
[Link]()

31
 From our graph, we can observe the relationship that exists between the explanatory variables.
and the dependent variable is, apparently, negative for both. Now, it is time to quantify such
relationships. We must assign the variables to objects, so that the construction becomes less tedious.
model and we have a more 'clean' code. To select variables from a DataFrame, you must use
the brackets '[ ]', taking into account the following:

When it comes to a single variable, it is enough to write its name (in quotes) inside the braces.
test.

When it comes to multiple variables, their names should be written (in quotes, separated by commas)
inside a Python list [...] and use such a list within the brackets.

The dependent variable will be assigned to object Y and the explanatory variables will be assigned to object X.
like this:

In [15]: Y=data["CM"]
X=data[["PGNP","FLR"]]

Now we have assigned the variables to the corresponding objects, we can now use such objects.
as arguments of the OLS(...) function of Statsmodels for building the model, remembering that its
the arguments are, respectively, the dependent variable and the explanatory variables:

In [16]: MyModel=[Link](Y, X)

Llevamos a cabo la estimación, con el mé[Link](), y guardamos los resultados generados en un objeto
what will be called Results:

In [17]: Results = [Link]()

32
 Finally, we visualize the results using the built-in (native) print(...) function, using co-
I argue the application of the [Link]():

In [18]: print([Link]())

OLS Regression Results

==============================================================================
Dep. Variable: CM R-squared: 0.387
Model: OLS Adj. R-squared: 0.368
Method: Least Squares F-statistic: 19.61
Date: xxx, xx xxx xxxx Prob (F-statistic): 2.52e-07
Time: xx:xx:xx Log-Likelihood: -400.07
No. Observations: 64 AIC: 804.1
Df Residuals: 62 BIC: 808.5
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef standard error t P>|t| 0.025 0.975
------------------------------------------------------------------------------
PGNP -0.0060 0.006 -0.974 0.334 -0.018 0.006
FLR 1.8838 0.325 5.796 0.000 1.234 2.534
==============================================================================
Omnibus: 7.742 Durbin-Watson: 1.601
Prob(Omnibus): 0.021 Jarque-Bera (JB): 2.919
Skew: 0.161 Prob(JB): 0.232
Kurtosis: 2.005 Cond. No. 62.1
==============================================================================

Warnings:
Standard Errors assume that the covariance matrix of the errors is correctly
specified.

It can be evidenced that the intercept term has not been estimated, that is, a regression has been carried out.
through the origin. This, as mentioned earlier, is because the OLS(...) function of Statsmodels
do not assume by default that a constant should be included in the set of explanatory variables. In order for the
the estimation is carried out with an intercept, it is necessary to specify that a constant must be used as a re-
aggressor; for this purpose, the function add_constant(...) from Statsmodels can be used, taking as
argument is the object that contains the other explanatory variables.

In [19]: MyModel2=[Link](Y, sm.add_constant(X))

The estimation and visualization of results are carried out using the methods that have been used.
previously (.fit() and .summary()):

In [20]: Results2 = [Link]()

print([Link]())

33
 OLS Regression Results
==============================================================================
Dependent Variable: CM R-squared: 0.708
Model: OLS Adj. R-squared: 0.698
Method: Least Squares F-statistic: 73.83
Date: xxx, xx xxx xxxx Prob (F-statistic): 5.12e-17
Time: xx:xx:xx Log-Likelihood: -328.10
No. Observations: 64 AIC: 662.2
Df Residuals: 61 BIC: 668.7
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef standard error t P>|t| 0.025 0.975
------------------------------------------------------------------------------
const 263.6416 11.593 22.741 0.000 240.460 286.824
PGNP -0.0056 0.002 -2.819 0.006 -0.010 -0.002
FLR -2.2316 0.210 -10.629 0.000 -2.651 -1.812
==============================================================================
Omnibus: 0.732 Durbin-Watson: 2.186
Probability (Omnibus): 0.693 Jarque-Bera (JB): 0.559
Skew: 0.228 Prob(JB): 0.756
Kurtosis: 2.949 Cond. No. 6770
==============================================================================

Warnings:
Standard Errors assume that the covariance matrix of the errors is correctly
specified.
The condition number is large, 6.77e+03. This might indicate that there
strong multicollinearity or other numerical problems.

The reader can observe that the results we have achieved have been obtained with the use of
methods, which are essentially functions associated with a specific object, directly applicable
about this.
.fit() y .summary() are methods applicable to objects generated by the OLS(...) function, however,
these are not the only methods specific to such types of objects. To know the methods and attributes of a
object in Python, the built-in (native) function dir(...) is used, which takes an object as an argument
from Python and returns its attributes and methods within a Python list '[...]'.
To know the methods and attributes of the objectResultados2, the code to use is the following:

In [21]: dir(Resultados2)

Out[21]:['HC0_se',
HC1_se
HC2_se
HC3_se
_HCCM
'__class__',

34
__delattr__
'__dict__',
'__dir__',
'__doc__',
'__eq__',
'__format__',
__ge__
__getattribute__
'__gt__',
'__hash__',
__init__
'__init_subclass__',
__le__
'__lt__',
'__module__',
__ne__
__new__
'__reduce__',
'__reduce_ex__',
'__repr__',
'__setattr__',
'__sizeof__',
__str__
'__subclasshook__',
'__weakref__',
'_cache',
'_data_attr',
_get_robustcov_results
'_is_nested',
'_wexog_singular_values',
aic
bic
bse
centered_tss
compare_f_test
compare_lm_test
compare_lr_test
'condition_number',
'conf_int',
'conf_int_el',
cov_HC0
cov_HC1
cov_HC2
cov_HC3
cov_kwds
'cov_params',
'cov_type',
df_model

35
 df_resid
diagn
eigenvals
the_test
is
f_pvalue
f_test
fitted values
fvalue
get_influence
get_prediction
get_robustcov_results
initialize
k_constant
'llf',
load
model
mse_model
'mse_resid'
'mse_total',
nobs
normalized_cov_params
outlier_test
params
predict
pvalues
remove_data
resid
resid_pearson
rsquared
adjusted_r_squared
save
scale
'ssr'
'summary',
summary2
't_test'
t_test_pairwise
tvalues
uncentered_tss
use_t
wald_test
wald_test_terms
wresid

Each of the objects in this Python list [...] is a method or attribute of the model we have built.
Thus, if we want to know the multiple determination coefficient oR2 , which is identified with the attribute
.rsquared, podemos emplear el siguiente código:

36
In [22]:print("The R2 of the model is:", [Link])

The model's R2 is: 0.7076654981900712

If we are interested in the estimated coefficients, we can resort to the attribute .params, as follows:

The coefficients of the model are:\n, [Link])

The coefficients of the model are:

const 263.641586
PGNP -0.005647
FLR -2.231586
dtype: float64

Likewise, the methods allow for the execution of specific actions on the object. For example, the method-
do.get_robustcov_results() allows you to obtain a new instance of results assuming by default
robust covariance.

In [24]: RobustResults2=Results2.get_robustcov_results()
print(Resultados2_robustos.summary())

OLS Regression Results

==============================================================================
Dependent Variable: CM R-squared: 0.708
Model: OLS Adj. R-squared: 0.698
Method: Least Squares F-statistic: 93.06
Date: xxx, xx xxx xxxx Prob (F-statistic): 2.94e-19
Time: 17:34:52 Log-Likelihood: -328.10
No. Observations: 64 AIC: 662.2
Df Residuals: 61 BIC: 668.7
Df Model: 2
Covariance Type: HC1
==============================================================================
coef standard error t P>|t| 0.025 0.975
------------------------------------------------------------------------------
const 263.6416 11.983 22.002 0.000 239.681 287.602
PGNP -0.0056 0.001 -4.788 0.000 -0.008 -0.003
FLR -2.2316 0.194 -11.516 0.000 -2.619 -1.844
==============================================================================
Omnibus: 0.732 Durbin-Watson: 2.186
Prob(Omnibus): 0.693 Jarque-Bera (JB): 0.559
Skew: 0.228 Prob(JB): 0.756
Kurtosis: 2.949 Cond. No. 6770
==============================================================================

Warnings:
Standard Errors are heteroscedasticity robust (HC1)

37
The condition number is large, 6.77e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

The reader is invited to explore the other available methods and attributes, so that they can know everything
the additional information with which you can deepen your analysis.
At this point, the third chapter of this work comes to a close; the reader has already had an introduction to regression.
simple and multiple linear regression in Python with Ordinary Least Squares estimation; likewise, it has observed
how to carry out a minimal exploratory analysis and how to create basic graphs. Now, just as
as mentioned in this chapter, under certain assumptions, OLS estimators exhibit statistical properties
attractive; such assumptions correspond to the topic addressed in the following chapter.

38
Chapter 4

Verification of assumptions

The OLS estimators have a number of attractive statistical properties given certain assumptions;
This statement is based on the Gauss-Markov Theorem, which, following Gujarati and Porter (2010), can...
to express as: "Given the assumptions of the classical linear regression model, the ordinary least squares estimators
squares, within the class of unbiased linear estimators, have minimum variance, that is, they are MVUE.
(p. 72).
MELI (BLUE in English) refers to 'Best Linear Unbiased Estimators' and is the qualification that
the OLS estimators are acquired when the assumptions of the classical linear regression model are met,
also known as the Gaussian model or standard linear regression model.
The classical linear regression model proposes, according to Gujarati and Porter (2010), the following assumptions:
tos

It is a linear regression model, that is, linear in the parameters.

2. Fixed values of independent values from the error term. This means that it requires cova-
launch 0 enteriand each variablexyes . It is decircov(ui ,xyes ) = 0∀ i,j
3. Average value of the disturbanceiequal to 0.E(ui ) = 0∀ I

4. Homoscedasticity or constant variancei This is, var(ui ) = σ2∀ i

5. No autocorrelation, or serial correlation, between the disturbances. That is, cov(ui ,uj ) = 0∀ i= j
6. The number of observations (n) must be greater than the number of parameters to be estimated. If there are variables
explanatory (excluding the intercept), then k must be estimated+ 1parameters (for a regression
with the intercept term), therefore, in such a case, n> k+ 1.
7. There must be variation in the values of the variables.
8. There should be no exact collinearity among the variables.

9. There is no specification bias, that is, the model is correctly specified.

In this chapter, the fulfillment of some of these assumptions is verified using various tools.
something that the reader must take into account is that some of the previously mentioned assumptions
specifically corresponds to the structure (equation) of the model (for example, assumption 1) while
others make explicit reference to the characteristics of the error term (for example, assumption 4 and the su-
position 5). Thus, for greater clarity, we will first examine the assumptions regarding the structure of the model
and we will subsequently address the assumptions related to the term error, so it will not continue the
same order in which the assumptions have been stated.

39
 Assumptions about the structure of the model
The classic linear regression model assumes the fulfillment of a series of assumptions about the structure.
(equation) of the model. It is possible to verify the compliance with some of these assumptions in Python using
functions specific to some specialized libraries like Statsmodels.
We will estimate a multiple linear regression model and verify the compliance of the budgets about
the structure of it. The dataset we will use will be, just like in the previous chapter, the one from the Example
7.1, Infant mortality in relation to GDP per capita and the literacy rate of women, by Gujarati and Porter
(2010).

Environment Preparation
The first cell in which the user will write and execute code will have the following content:

In [1]:import numpy as np
importpandasaspd
[Link]
[Link]
import [Link] as plt
importseabornassns
%matplotlibinline
[Link]("seaborn-white")

The reader should already know that the previous block of code allows loading the necessary tools for
carry out the required actions in the proposed econometric analysis; in case there is no clarity
In this regard, you are asked to consult the 'Environment Preparation' section of the previous chapter, where it is provided
an explanation about such an issue.
The code in this cell is exactly the same as that used in the environment setup of the previous chapter,
except for the following: we will additionally use the stats module from the Statsmodels library, which
it will be imported with the alias desms.
The import of [Link] is fundamental to our work, as it is where they are conte-
We will use most of the functions, without which it would not be possible to easily verify compliance.
falsification of the assumptions of the classical linear regression model. The user must be careful in case of
reuse the code block for setting up the environment from the previous chapter, correctly include the line of
import [Link].

Data import
The data importation is carried out with functions from the pandas library. Depending on the type of
file in which the information of the dataset is stored, we will have to resort to a partition function
and a set of specific parameters.
For our specific case, the following code is the one that imports the dataset:

In [2]:data=pd.read_csv("[Link]",
delimiter="\t")

To the reader who does not understand this line of code, it is suggested, again, to refer to Chapter 2, where
you will find an explanation that will surely clarify your doubts, or to Chapter 1 of Triana and Galindo (2019).

40
 Once the line of code in this cell has been executed, the dataset should have been imported correctly.
To verify it, we resort to the method .head(), applied to the object that contains the dataset (in this case it
we have assigned the name data):

In [3]:[Link]()

Output[3]: CM FLR PGNP TFR

0 128 37 1870 6.66
1 204 22 130 6.15
2 202 16 310 7.00
3 197 65 570 6.25
4 96 76 2050 3.81

We can observe that the import process has been successful and the information from the dataset has been soul-
You can correctly fill in a DataFrame in pandas. Now, we can continue peacefully with the development.
from our analysis, as we have satisfactorily completed the first step. We will proceed to verify the
compliance with some of the assumptions regarding the structure of the model.

Assumption of number of observations greater than number of parameters (Assumption #6)

One of the assumptions about the structure of the model is that the number of observations with which it is performed
the estimation must exceed the number of parameters (βto be estimated. If we intend to work with all
The observations of the DataFrame, it is enough to examine its dimensions to determine the maximum amount.
of parameters that could be estimated. The dimensions of a DataFrame correspond to the number of rows and
columns by which it is formed and can be known very easily (remember that it has already been done
same operation in the second chapter): just use the attribute .shape, like this:

In [4]:[Link]

(64, 4)

The result we obtain is a Python tuple (...) that informs us that our DataFrame is of size
64× 4, so it has 64 rows and 4 columns. Considering that each row corresponds to an observation and
each column to a variable, the maximum number of parameters to estimate is 63.
What happens if the number of observations is less than the number of parameters to be estimated? 'If there are fewer
observations thenXcannot be of full rank” (Greene, 2003, p. 14). Because of this, when
the number of parameters to estimate exceeds the number of available observations to make the estimation,
it is obtained that the variance is infinite, therefore the OLS method cannot be used (James, Witten, Hastie, and
Tibshirani, 2013).
As an example, a test regression will be performed, for purely illustrative purposes, in which
They will use only 3 observations to estimate 4 parameters. For this, only the first three will be selected.
Observations of the DataFrame using brackets [ ]:

In [5]: Test_Model=[Link](data["CM"][0:3],
sm.add_constant(data[["FLR", "PGNP", "TFR"]][0:3])
Resultados_de_Prueba=Modelo_de_Prueba.fit()
print(Test_Results.summary())

41
 OLS Regression Results
==============================================================================
Dependent Variable: CM R-squared: 1.000
Model: OLS Adj. R-squared: nan
Method: Least Squares F-statistic: 0.000
Date: xxx, xx xxx xxxx Prob (F-statistic): nan
Time: xx:xx:xx Log-Likelihood: 66.160
No. Observations: 3 AIC: -126.3
Df Residuals: 0 BIC: -129.0
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef standard error t P>|t| 0.025 0.975
------------------------------------------------------------------------------
const 3.2453 inf 0 nan nan nan
FLR 2.0501 inf 0 nan nan nan
PGNP -0.0692 inf -0 nan nan nan
TFR 26.7721 inf 0 nan nan nan
==============================================================================
Omnibus: nan Durbin-Watson: 0.854
Prob(Omnibus): nan Jarque-Bera (JB): 0.511
Skew: -0.678 Prob(JB): 0.774
Kurtosis: 1,500 Cond. No. 753.
==============================================================================

Warnings:
Standard Errors assume that the covariance matrix of the errors is correctly
specified.
The input rank is higher than the number of observations.

As the reader can observe, the standard errors are infinite, and the p-values are not confidence intervals.
are defined, so it is not useful to use OLS estimators when the number of observations
is less than the number of parameters that are intended to be estimated.
Certainly, the example that has just been presented is very simple; however, it is possible that the dataset
to employ a large number of explanatory variables, even far greater than the number of observations
available, which is why Assumption #6, despite its innocent appearance, should not be taken lightly.
When there are many more explanatory variables than observations, various techniques can be used.
and variable selection processes can be carried out; however, such procedures exceed the
scope of this work and will not be addressed therein.

Assumption of variation in the values of the explanatory variables (Assumption #7)

El supuesto #7 indica que las variables explicativas deben ser “variables”; es decir, que los valores de las
explanatory variables must change and cannot correspond to the same constant. Intuitively, for
an explanatory variable be 'explanatory' its changes must be associated1to changes in the dependent variable;
1 This does not mean that there is a causal relationship. The fact that the movement of one variable is associated with the movement of another

does not mean that this is causing the change of the other.

42
if the value of an explanatory variable does not change, it becomes difficult to identify the way it is associated with the
dependent variable.
Identifying compliance with this assumption is quite simple: 'if the sample standard deviation of the
xiif it is zero, then the assumption is not satisfied; if not, this assumption is satisfied” (Wooldridge, 2010, p.
48).
To verify the assumption of variability in the variables in Python, one can use,
for example, the method .apply(), using the function std(...) from NumPy:

In [6]: [Link]([Link]([Link], axis= 0),

columns=["Desviación estándar"])

Out[6]: Standard deviation

CM 75.382151
FLR 25.803873
PGNP 2704.317439
TFR 1.497158

Note: The function DataFrame(...) is used in this code for purely aesthetic purposes, to
obtain a more pleasant presentation of the result. The really important part, which generates the
content of interest, [Link]([Link], axis = 0).

As evidenced, all standard deviations are different from 0, so there is variability in

the variables (if the standard deviation of any variable were 0, it would not be a variable, but rather a
constant).

Assumption of No Multicollinearity (Assumption #8)

One of the assumptions about the structure of the model is that there is no exact collinearity among the
explanatory variables. To the reader who is unaware of the definition of collinearity, it is noted that this, according to
Dormann et al. (2013), 'describe the situation in which there are two or more predictor variables in a statistical model'
are linearly related” (p. 28).
One method to detect collinearity between variables is to calculate the simple correlation coefficient between
pairs of these; if any of these coefficients is greater, in absolute value, than 0.9, then there is a symptom of
very strong collinearity.
To examine the correlation between the numerical variables of a pandas DataFrame, it is possible to use
the [Link](). Like this:

In [7]: [Link]()

Out[7]: CM FLR PGNP TFR

CM 1.000000 -0.818285 -0.407697 0.671135
FLR -0.818285 1.000000 0.268530 -0.625954
PGNP -0.407697 0.268530 1.000000 -0.185718
TFR 0.671135 -0.625954 -0.185718 1.000000

If the user prefers a somewhat more visually 'rich' representation, they can resort to the function
heatmap(...) from Seaborn, which draws a color-coded matrix. The argument to be passed to the function,
for this specific case, it is the correlation matrix, which, as we have just seen, is obtained with
the application of the method .corr() on the DataFrame.
The code to be used is the following:

43
In [8]: fig, ax = [Link](figsize=(8,6))
Simple correlation\nbetween explanatory variables
fontsize= 18
fontweight="bold", x= 0.43
[Link]([Link](), ax=ax)
plt.subplots_adjust(top= 0.85)
[Link](.9,-.02,
Preparation:
fontsize= 13, fontweight="bold"
ha="right")
[Link](.9,-.07,
Triana, F.\n(2019)
fontsize= 12, ha="right"
[Link]()

44
 Note: The reader should note that, since Seaborn is a library based on Matplotlib and that the function
heatmap(...) de Seaborn returns an Axes of a Matplotlib figure, the code used in the block-
The previous code is common for this latest library. Since Matplotlib plays an auxiliary role in
in this work, no efforts will be devoted to thoroughly examining the structure of its code, so
To the interested reader in this, it is suggested to consult additional resources.

We can observe that the default color scale used by the heatmap(...) function of Seaborn does not
seems to be the most suitable. It is possible to modify this by including a specific value for the cmap parameter,
assigning the scale that seems most appropriate to us. As an example, the following code will use the
scale 'RdYlGn' (reds, yellows, greens):

In [9]:fig, ax=[Link](figsize=(8,6))
[Link]("Simple correlation\namong explanatory variables
fontsize= 18,
fontweight="bold", x= 0.43
[Link]([Link](), ax=ax, cmap="RdYlGn")
plt.subplots_adjust(top= 0.85)
[Link](.9,-.02,
"Elaboración:",
fontsize= 13, fontweight="bold"
ha="right")
[Link](.9,-.07,
Triana, F.\n(2019)
fontsize= 12, ha="right"
[Link]()

45
 If we consider that our graph is still not informative enough, we can add the coef-
correlation coefficient corresponding to each of the cells. To achieve this goal, we must
include the parameter annotate in the function heatmap(...) of Seaborn and assign it the value of True (annotate a
boolean type parameter). Thus:

In [10]:fig, ax=[Link](figsize=(8,6))
Simple correlation\namong explanatory variables
fontsize= 18,
fontweight="bold", x= 0.43
[Link]([Link](), ax=ax, cmap="RdYlGn"}
annot=True)
plt.subplots_adjust(top= 0.85)
[Link](.9,-.02,
Preparation:

46
 fontsize= 13, fontweight="bold"
ha="right")
[Link](.9,-.07,
Triana, F.\n(2019)
fontsize= 12, ha="right"
[Link]()

Thus, we note that the strongest correlation occurs between the variables 'CM' and 'FLR', with a coefficient
greater than 0.8 (in absolute value); however, such a magnitude of the coefficient should not be a concern
however, as for Hypothesis #8 regarding collinearity between explanatory variables and the variable 'CM', it is the
dependent variable, while 'FLR' is an explanatory variable.
Assumption #8 refers to exact collinearity; this is when an explanatory variable has a re-
exact linear relationship with another. When exact collinearity is present, the regression coefficients are indeterminate.
mined and the standard errors are infinite (Gujarati and Porter, 2010).

47
 In our specific case, none of the explanatory variables shows a very strong correlation with
any other explanatory variable different from itself, therefore, perfect multicollinearity is not present and
It is possible to estimate the parameters.
Another useful technique for detecting multicollinearity is through the Variance Inflation Factor (VIF,
high VIFs reflect an increase in the variances of the regression coefficients
dear due to collinearity between predictor variables, compared to those that would be obtained when
the predictors are orthogonal” (Murray, Nguyen, Lee, Remmenga, and Smith, 2012).
Thus, in case there is no multicollinearity, or it is mild, it is expected that the magnitude of the VIFs will be
smaller and the stronger the collinearity, the higher the VIF will be; Gujarati and Porter (2010), based on the
Kleinbaum, Kupper, and Muller's work (1988) suggests a practical rule: "if the VIF of a variable is
greater than 10 (this happens if2j if it exceeds 0.90), it is said that the variable is very collinear" (p.340).
What are the VIFs of the variables in the dataset being worked on? To obtain these values,
the function variance_inflation_factor(...) from the module stats.outliers_influence will be used
from the Statsmodels library.
The function variance_inflation_factor(...) receives two arguments: the matrix that contains the values-
explanatory variables (exog) and an index (exog_idx) that indicates which variable corresponds to the calculated VIF.
This function calculates the VIF for the specified variable; however, it is of interest to know the VIF of each one of
the explanatory variables.
It is possible to apply the function variance_inflation_factor(...) manually as many times as
interest variables there; for our case, there are 3 explanatory variables, so the function should only apply-
do it 3 times, an operation that is not excessively tedious. However, if one had a considerable number
Of explanatory variables, it would be annoying to repeat the operation multiple times, so it is useful
to "automate" such a task. One option to achieve this is to use a for loop:

In [11]: [Link].outliers_influenceimportvariance_inflation_factor
vif=[]
exog=sm.add_constant(data[[Link][1:4]])
range=sm.add_constant(data[[Link][1:4]]).shape[1]
foriinrange
[Link](variance_inflation_factor([Link], i))
[Link]({'VIF': vif[1:]}, index=[Link][1:4])

Out[11]: VIF
FLR 1.711845
PGNP 1.078306
TFR 1.645150

Thus, it can be observed that none of the explanatory variables exhibit a high VIF, so it does not
suggests that some of these variables are highly collinear. This result is in line with what was obtained in the
correlation matrix, in which it was found that none of the simple correlation coefficients were very
high.

Assumption of No Specification Bias (Assumption #9)

One of the assumptions about the model is that it must be correctly specified, that is, that the
the functional form it adopts is the correct one and that the included variables are appropriate.
Visually, it is possible to examine if the adopted functional form is appropriate for graphing, through
a scatter plot, the estimated values against the residuals; it is expected that the behavior is

48
relatively stable, with points distributed symmetrically around zero, without showing patterns
specific.
Since the model has not yet been created nor the estimation carried out, there is no instance of
results to be used. Therefore, the first thing to be done is to carry out the construction and corresponding estimation.
dientes:

In [12]: Y=data["CM"]
X=data[["PGNP","FLR"]]
Model=[Link](Y, sm.add_constant(X))
Results = [Link]()
print([Link]())

OLS Regression Results

==============================================================================
Dep. Variable: CM R-squared: 0.708
Model: OLS Adj. R-squared: 0.698
Method: Least Squares F-statistic: 73.83
Date: xxx, xx xxx xxxx Prob (F-statistic): 5.12e-17
Time: xx:xx:xx Log-Likelihood: -328.10
No. Observations: 64 AIC: 662.2
Df Residuals: 61 BIC: 668.7
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef standard error t P>|t| 0.025 0.975
------------------------------------------------------------------------------
const 263.6416 11.593 22.741 0.000 240.460 286.824
PGNP -0.0056 0.002 -2.819 0.006 -0.010 -0.002
FLR -2.2316 0.210 -10.629 0.000 -2.651 -1.812
==============================================================================
Omnibus: 0.732 Durbin-Watson: 2.186
Prob(Omnibus): 0.693 Jarque-Bera (JB): 0.559
Skew: 0.228 Prob(JB): 0.756
Kurtosis: 2.949 Condo. No. 6770
==============================================================================

Warnings:
Standard Errors assume that the covariance matrix of the errors is correctly
specified.
The condition number is large, 6.77e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

The model and its results are already available, but we have not yet assigned the residuals to an object, so
we will carry out such assignment, so that the code used is simplified. We will assign the residues,
which are an attribute of the results instance (the user can verify it with the funcióndir(...)),
to the objectWaste, like this:

In [13]:Residuos=[Link]

49
 Similarly, it is possible to assign the estimated values obtained from the application of the method.
.predict() to a specific object:

In [14]:Valores_estimados=[Link]()

In [15]:fig, ax=[Link](figsize=(8,6))
Estimated values vs.\nWaste
fontsize= 18,
fontweight="bold"
[Link](Estimated_values, Residuals,
color="darkblue"
ax.set_xlabel("Estimated values", fontsize= 15)
ax.set_ylabel("Residuals", fontsize= 15)
ax.tick_params(labelsize= 13)
[Link](.9,-.02,
Preparation:
fontsize= 13, fontweight="bold"
ha="right")
[Link](.9,-.07,
Triana, F.\n(2019)
fontsize= 12, ha="right"
[Link]()

50
 As observed, there does not appear to be a specific pattern and the points are symmetrically distributed.
builds around zero, so the visual analysis suggests that the adopted functional form is not incorrect.
Another aspect to consider regarding the structure of the model is that it should be maintained throughout the esti-
information; that is, not being unstable. To verify the stability in the parameters, or no structural change,
The CUSUM test based on OLS residuals can be used, which takes no change as the null hypothesis.
structural.
The function used in Python for the CUSUM test based on OLS residuals is
breaks_olsresid(...)[Link]; its arguments include the residuals of the estima- model
made by MCO and the ddof parameter, with a default value equal to 0 and referring to the number of degrees of freedom
employees in the correction of error variance.
Now, we will use the objectResiduos as an argument of the function breaks_olsresid(...); the
results obtained from the execution of the function will be stored in an object that we will call
TestResults, like this:

In [16]: TestResults=sms.breaks_cusumolsresid(Residues)
print(ResultadosTest)

51
(0.5191974341185455, 0.9503227705917948, [(1, 1.63), (5, 1.36), (10, 1.22)])

We can see that the result we get from our code is just a Python tuple (...) that
it contains a dataset; this is because [Link] does not have a built-in presentation mode
more "elaborate" of the result.
To improve the presentation of the Test results, a pandas Series can be used.
so that we can properly organize the available information. We will proceed to create said series
assigning as values the components of the test result and as a Python unist index [...] with its
respective names. Like this:

In [17]:Nombres=["Estadístico","p-value del estadístico","Valores críticos"]

[Link](TestResults, index=Names)

Statistic 0.519197
p-value of the statistic 0.950323
Critical values [(1, 1.63), (5, 1.36), (10, 1.22)]
dtype: object

Assuming a significance level of 5%, how− value> α then the null hypothesis is not rejected
it is concluded that there is no structural change (at a 95% confidence level).
Some of the assumptions of the classical linear regression model have been briefly evaluated regarding
the structure of the model. Now, the compliance with the assumptions regarding the error term will be examined.

Assumptions about the term of error

The model used in the section 'Assumptions about the structure of the model' will be continued, for which
that the 'Preparation of the environment' and 'Data Import' will not take place, as these have already been completed
successfully carried out previously. The reader who has not followed the information presented is advised
go to the beginning of this chapter and run the first two code cells; once this is done, you should be able to
continue with the presentation without any issues.

Average value of the disturbance equal to zero (Assumption #3)

The residuals of the model can be positive or negative, that is, it is possible that the estimated value by
the model may be inferior to the true one or superior to it, but the distribution must be relatively symmetrical.
It is not expected that all residues are positive or that all are negative, as a ...
systematic underestimation or overestimation; it is appropriate that there is a balance between positive and negative residuals.
negatives, so that their average value is zero.
This assumption is relatively easy to evaluate, as we only need to calculate the average value of the
waste, for which the function mean(...) of NumPy is suitable:

In [18]: [Link](Residuals)

1.254552017826427e-13

The argument that the function mean(...) of NumPy receives is the vector of residues, which, conveniently-
We have previously assigned to the waste object. The result obtained is practically zero, because

52
It is a number, expressed in scientific notation, that is very small; this suggests that Assumption #3
this holds true for this example.
Although we have used a function, it may be helpful to point out to the reader that there is also a method.
with which it is possible to calculate the average, which, as expected, is .mean():

In [19]: [Link]()

Out[19]: 1.254552017826427e-13

Homoscedasticity (Assumption #4)

The assumption of homoscedasticity in the error term essentially indicates that it must exhibit variance.
constant variance. There are formal statistical tests to verify compliance with this assumption, without
embargo, the first thing we will do is use informal methods, so that the reader has an idea of
all the clearer about the issue to be addressed.
We will examine the presence of heteroscedasticity (that is, non-constant variance or 'non-
homoscedasticity") informally through visual analysis. What we will do is graph the
residues of the model, with the aim of observing if their behavior exhibits a particular pattern that may
indicate the presence of heteroscedasticity. (See Section 11.5, Chapter 11, Heteroscedasticity: what
What happens if the variance of the error is not constant? from Gujarati and Porter (2010)
To create the corresponding graph we will use the method .plot() on an Axes created with the function
subplots(...) from Matplotlib, using the model's residuals as an argument; these residuals have been
previously assigned to an object, which is why this object will be used as an argument. To the reader who
If you have not followed the exercise carried out in the section 'Assumptions about the model structure', you are invited to
consult it and execute the corresponding code cells, with which you can generate the required objects
next.
The content of the code cell for generating the graph is as follows:

In [20]:fig, ax=[Link](figsize=(10,5))
Figure title: 'Residuals'
fontweight="bold"
[Link](Residuals, color="darkblue")
ax.set_xlabel("Observations", fontsize= 15)
ax.set_ylabel("Residuals", fontsize= 15)
ax.tick_params(labelsize= 14)
[Link](.9,-.02,
"Elaboración:",
fontsize= 13, fontweight="bold"
ha="right")
.9,-.08,
Triana, F.\n(2019)
fontsize= 12, ha="right"
[Link]()

53
 As can be seen, the residues do not exhibit a particular pattern and the range of their changes does not stop.
It may vary significantly, so at first glance, there is no evidence of heteroscedasticity. However,
That conclusion may be hasty.
A similar analysis can be carried out using the estimated values and the square of the residuals.
that allows controlling by the nature of the waste (positive or negative) based on its magnitude and not
in its sign. It is expected that no particular pattern will emerge, so the variance of the residuals will be
constant and these do not depend on the magnitude of the returned variable.
To execute this analysis we will use the method .scatter(...) using as arguments the vec-
of the estimated values and that of the squared residues. We already have the former, but we only have information on
formation of the waste and not of the square waste; however, it is sufficient to apply the operator ** to
obtain these magnitudes, since this is the one that allows in Python to perform exponentiation following the struct-
turabase**exponent:

In [21]:fig, ax=[Link](figsize=(10,5))
[Link]("Estimated values vs."\nSquare residues
fontsize= 18,
fontweight="bold"
[Link](Estimated_Values, Residuals**2,
color="darkblue"
ax.set_xlabel("Estimated values")\n(dependent variable)
fontsize= 15)
ax.set_ylabel("Squared residues", fontsize= 15)
ax.tick_params(labelsize= 14)
[Link](.9,-.02,
Preparation:

54
 fontsize= 13, fontweight="bold"
ha="right")
[Link](.9,-.08,
Triana, F.\n(2019)
fontsize= 12, ha="right"
plt.subplots_adjust(top= 0.85)
[Link]()

According to the result obtained, it does not seem that the residues are heteroscedastic, although there is a li-
It gives the impression that they are growing with the magnitude of the dependent variable; to reach a conclusion.
"valid", visual analysis is not sufficient: formal statistical tests are required.
One of the statistical tests commonly used to verify the presence of heteroscedasticity
it is the Breusch-Pagan Lagrange Multiplier Test (Greene, 2003), which basically acts as
null hypothesis homoscedasticity.
The function used in Python for the Breusch-Pagan Lagrange Multiplier Test for heteros-
Cauchy ishet_breuschpagan(...) from [Link]; its arguments include the residuals of the mo-
work and the set of variables that can cause heteroscedasticity.
A container for the waste is already available, however, the variables have not yet been assigned.
explanatory variables of the model to a particular object. As the explanatory variables consist of an attribute of the
instance of model results, we can assign them to an object by applying this attribute, as follows:

Explanatory = [Link]

Since we have the objects to use as arguments, we can make use of the function
the_breuschpagan(...)[Link]. Following the structure we have used for the tests
From the previous section, the code to execute is as follows:

55
In [23]: TestResults = sms.het_breuschpagan(Residues, Explanatory)
Nombres=["Estadístico LM","p-value del estadístico LM","Estadístico F",
p-value of the F statistic
[Link](TestResults, index=Names)

LM Statistic 1.591582
p-value of the LM statistic 0.451224
F Statistic 0.777832
p-value of the F statistic 0.463904
dtype: float64

Assuming a significance level of 5%, how− value> α then the null hypothesis is not rejected
it is concluded that homoscedasticity is present in the error term (working with aα of 0.05).

Note 1: This test is equivalent to the one generated by the function bptest(...) in R (package lmtest).

Note 2: This test is equivalent to the one generated by the code estat hettest, rhs fstatdeStata.

Another commonly used statistical test to verify the presence of heteroscedasticity

Lagrange multiplier test, which makes use of auxiliary regressions, essentially takes
as the null hypothesis of homoscedasticity, and is able to identify more general forms of heteroscedasticity
the Breusch-Pagan Test (Verbeek, 2004).
The function used in Python for the White Test of Lagrange Multipliers for heteroscedasticity
eshet_white(...) [Link]; your arguments include the model's residuals and the set of
variables to be used in the auxiliary regressions.
The arguments used in the function het_white(...) from [Link] correspond to the same.
used in the functionhet_breuschpagan(...)of that module and, since these have been assigned
previously to specific objects, it is enough to reuse such objects.
Following the structure that has been used for the other statistical tests conducted, the code to
execute is the following:

In [24]: TestResults = sms.het_white(Residues, Explanatory)

Nombres=["Estadístico LM","p-value del estadístico LM",
"Estadístico F","p-value del estadístico F"]
[Link](TestResults, index=Names)

LM Statistic 2.356718
p-value of the LM statistic 0.797902
F Statistic 0.443486
p-value of the F statistic 0.816252
dtype: float64

The results of the White Test for Lagrange Multipliers for heteroscedasticity reaffirm the
conclusion of the Breusch-Pagan Test conducted previously. Assuming a significance level of 5%, as
p− value> α then the null hypothesis is not rejected and it is concluded that homoscedasticity is present in
the error term (at a 95% confidence level).

Note: This test is equivalent to the one generated by the code estat imtest, white in Stata.

56
No autocorrelation, or serial correlation, between the disturbances (Assumption #5)
The assumption of no autocorrelation among the disturbances indicates that they do not follow systematic patterns.
and they are not correlated with each other; such an assumption can be justified in the case of cross-sectional data, but
tend to be violated when working with time series, as successive observations usually
are strongly correlated (See Section 3.2, Chapter 3, Regression models with two variables:
estimation problem, by Gujarati and Porter (2010)).
This assumption is usually verified with statistical tests, some of which are addressed as continuous.
A statistical test commonly used to verify the presence of first-order autocorrelation.
it is the Durbin-Watson Test, which essentially takes as the null hypothesis no autocorrelation (of order 1) and
generate a statistic whose value is contrasted with the corresponding critical values (Chatterjee and Simonoff,
2013) to establish a conclusion.
The function used in Python for the Durbin-Watson Test is durbin_watson(...) from [Link];
his arguments include the model's waste.
Following the same structure as the other statistical tests performed, and taking into account that the
wastes have already been assigned to a specific object, the code to execute is as follows:

In [25]:ResultadosTest=sms.durbin_watson(Residuos)
Nombres=["Estadístico DW"]
[Link](TestResults, index=Names)

DW Statistic 2.186159
dtype: float64

The obtained value is 2.1861, then, is there first-order autocorrelation? The reader should note that, at
Unlike other statistical tests, in the case of the Durbin-Watson Test, a p-value is not reported with which
decide whether to reject the null hypothesis or not, then how is the conclusion determined?
The Durbin-Watson statistic is compared with critical values from which zones are established.
of positive autocorrelation, no autocorrelation, indeterminacy, and negative autocorrelation; depending on
The area in which the value of the statistic is located leads to a specific conclusion.
For a specific case, assuming a significant design level of 5%, the values of dLydU for2variables
explanatory (excluding the constant term) and 65 observations are, respectively, 1.536 and 1.662. The
The value of 2.1861 belongs to the interval(dU 4− dU ), which corresponds to the No Autocorrelation zone, for
it is concluded that there is no first-order autocorrelation (working with aα from 0.05).
Note: This test is equivalent to the one generated by the function dwtest(...) of R (package lmtest). The user must
keep in mind that, in the R function, unlike that of Python, the argument to be used is not the
waste unless the model itself.

The Durbin-Watson Test is a statistical test that allows the identification of first-order autocorrelation.
however, if one wants to examine higher-order serial correlation, this can no longer be employed and gene-
Typically, the Breusch-Godfrey Test is used, which essentially takes no autocorrelation as a hypothesis.
null and is based on the execution of an auxiliary regression in which lags of the error term are included as
explanatory variables (Kleiber and Zeileis, 2008).
The function used in Python for the Breusch-Godfrey Test of Lagrange Multipliers for autocorrelation
relationship of the residuals isacorr_breusch_godfrey(...)[Link]; its arguments include the
instance of results (not the residuals!) of the model and the number of lags, nlags, to include in the regression
auxiliary.
Following the same structure as the other statistical tests conducted, and taking into account that the
The instance of the model results has already been assigned to a specific object, the code to execute is the following:

57
In [26]:TestResults=sms.acorr_breusch_godfrey(Results, nlags= 1)
Nombres=["Estadístico LM","p-value del estadístico LM","Estadístico F",
p-value of the F statistic
[Link](TestResults, index=Names)

LM Statistic 0.728605
p-value of the LM statistic 0.393336
F statistic 0.690933
p-value of the F statistic 0.409143
dtype: float64

The results of the Breusch-Godfrey Lagrange Multiplier Test for autocorrelation of the residuals
two reaffirm the conclusion of the Durbin-Watson test conducted earlier. Assuming a level of signi-
financing of 5%, how− value> α then the null hypothesis is not rejected and it is concluded that it does not occur
first order serial correlation in the error term (working with a significance level of 5%).
It has already been indicated that the Breusch-Godfrey test is a general autocorrelation test, in that it allows
identify higher-order serial correlation greater than 1. In the previous code, we assigned the value of 1 to the parameter
metronlags of the function acorr_breusch_godfrey(...) from [Link], which we have examined
the presence of first-order autocorrelation, just like in the Durbin-Watson Test. Now, we will proceed
to examine the presence of autocorrelation up to order 4, for which it will take the value of 4. It is enough to
reuse the last code cell we have used and modify the value of the parameter nlags as follows:

In [27]: TestResults = sms.acorr_breusch_godfrey(Results, nlags=4)

Nombres=["Estadístico LM","p-value del estadístico LM","Estadístico F",
p-value of the F statistic
[Link](TestResults, index=Names)

LM Statistic 1.211387
p-value of the LM statistic 0.876220
F Statistic 0.274927
p-value of the F statistic 0.892978
dtype: float64

For the case of autocorrelation up to order 4, assuming a significance level of 5%, as−
value> α then the null hypothesis is not rejected and it is concluded that there is no serial correlation up to order
4 (working with aα of 0.05).
Note: This test is equivalent to the one generated by the bgtest(...) function of R (lmtest package).

Other statistical tests commonly used to examine the presence of autocorrelation in the
The error term is the TestBox-Piercey and the TestLjung-Box, which are closely related.
as a null hypothesis, basically, no autocorrelation and, like the Breusch-Godfrey Test, are useful for
determine the existence of autocorrelation up to order p.
The function used in Python for the Ljung-Box Test is the same as for the Box-Pierce Test. It is about
deacorr_ljungbox(...) [Link]; its arguments include the model's residuals, the number
of desired lags and a boolean variable to indicate whether only the results are wanted
of the Ljung-Box Test (performed by default) or if you also want to view the results of the Box-Pierce Test
(boxpierce = True o False).
Following the same structure as the other statistical tests performed, and taking into account that the
Waste has already been assigned to a specific object, the code to be executed is as follows:

58
In [28]: TestResults=sms.acorr_ljungbox(Residuals, lags= 1, boxpierce=True)
Nombres=["Estadístico LB","p-value del estadístico LB","Estadístico BP",
p-value of the BP statistic
[Link](TestResults, index=Names)

Statistic LB 0.7549153999366602
p-value of the LB statistic [0.3849244343516546]
BP statistician 0.7206010635759029
p-value of the BP statistic 0.3959468201677663
dtype: object

Both the p-value corresponding to the Ljung-Box Test statistic and the one corresponding to the Box test
The p-values are greater than 0.05, so the null hypothesis is not rejected and it is concluded that there is no correlation.
First-order serial relationship (note that lags = 1).
Just like in the Breusch-Godfrey Test, the Ljung-Box Test and the Box-Pierce Test can be used for
examine the existence of serial correlation up to order p. To examine up to order of autocorrelation,
Just assign the value of the parameter lags of the function acorr_ljungbox(...) from [Link].
Following the same structure as the other statistical tests conducted, and considering that the
wastes have already been assigned to a specific object, the code to execute to evaluate serial correlation up to
order 2 (for example) is as follows:

In [29]: TestResults = sms.acorr_ljungbox(Residuals, lags= 2, boxpierce=True)

Nombres=["Estadístico LB","p-value del estadístico LB","Estadístico BP",
p-value of the BP statistic
[Link](TestResults, index=Names)

Statistic LB [0.7549153999366602, 1.080970233422208]

p-value of the LB statistic [0.3849244343516546,0.5824656200187786]
BP Statistic [0.7206010635759029, 1.0268949980623265]
p-value of the BP statistic [0.3959468201677663, 0.5984289353087118]
dtype: object

The reader should note that the function corr_ljungbox(...) from [Link] generates p-values for the
statistical p-value corresponding, which it presents within a Python list "[...]", as it executes
the respective test for each of the included delays.

Note: This test is equivalent to the one generated by the function [Link](...) of R (package stats). Unlike
In R, the default test executed is the Box-Pierce test, so if you want to take it to
The Ljung-Box test is completed, such preference must be explicitly indicated in the 'type' parameter of the function.

Assumption of Normality
One of the assumptions about the error term is that it has a normal probability distribution. It is possible
that the careful reader has noticed that this assumption was not included in the list of assumptions of the model
classic linear regression, presented in the first part of this chapter. This is because the assumption of
normality in the distribution of the error term is not part of the classical linear regression model, because
what is not necessary to guarantee the properties of OLS estimators of linearity, unbiasedness and
minimum variance. So, why is it mentioned?

59
 The assumption of normality in the distribution of the error term is part of the classical regression model.
linear normal, and its existence is justified by the need to carry out statistical inference. While the
the other assumptions presented at the beginning of this chapter are required for the OLS estimators to be
MELI (or BLUE in English), such assumptions are not sufficient to guarantee that statistical inference is
valid.

Since the OLS method makes no assumptions about the probabilistic nature of the
Tica gavei it is of little help for the purpose of making inferences about the FRP through the
FRM, despite the Gauss-Markov theorem. This gap can be filled if it is assumed that the following
a certain probability distribution. For reasons that we will mention shortly, in the
contexto de regresión se supone, por lo general, que lasutienen la distribución de probabilidad
normal. (Gujarati and Porter, 2010, p. 98)

The assumption indicates that the error term has a normal distribution with a mean of zero and constant variance.
j) it is zero. Thus, the assumption can be summarized asI∼
and that, in addition, the covariance betweeniyouj(i=
NID (0, σ2 ), in other words: the error term is normally and independently distributed.
One of the most popular statistical tests used to verify the compliance of normality in
the distribution of the error term is the Jarque-Bera Test, which is based on skewness and kurtosis (Würtz and
Katzgraber, 2009) and takes as a null hypothesis, basically, normal distribution.
The function used in Python for the Jarque-Bera test for normality is jarque_bera(...)
[Link]; its arguments include the model's residuals.
Following the same structure as the other statistical tests conducted, and taking into account that the
residues have already been assigned to a specific object, the code to execute to evaluate normality in the distribution
the meaning of the error is as follows:

In [30]: TestResults = sms.jarque_bera(Residuals)

Nombres=["Jarque-Bera","p-value","Asimetría","Curtosis"]
[Link](TestResults, index=Names)

Jarque-Bera 0.559405
p-value 0.756009
Asymmetry 0.227575
Kurtosis 2.948855
dtype: float64

Assuming a significance level of 5%, how− value> α then the null hypothesis is not rejected
and it is concluded that normality is present in the distribution of the error term (working with aα of 0.05).
Note: This test is equivalent to the one generated by the function [Link](...) from R (package tseries).
the user must take into account that the argument used in the R function, just like in Python,
corresponds to the residues of the model.

Other widely used tests to assess normality are the Anderson-Darling test, which belongs to...
to the quadratic class of statistics based on the Empirical Distribution Function, and that of Kolmogorov-
Smirnov, also built from that function (Razaliy Wah, 2011). The first test is implemented
in Python with the function normal_ad(...) from [Link]:

In [31]: TestResults = sms.normal_ad(Residuals)

Nombres=["Anderson-Darling","p-value"]
[Link](TestResults, index=Names)

60
C:\Users\FCE\Anaconda3\lib\site-packages\statsmodels\stats\_adnorm.py:66:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`.
In the future this will be interpreted as an array index, `arr[[Link](seq)]`.
which will result either in an error or a different result.
S = [Link]((2*i[sl1]-1.0)/N*([Link](z)+[Link](1-z[sl2])), axis=axis)

Anderson-Darling 0.329856
p-value 0.508773
dtype: float64

Reinforcing the conclusion of the Jarque Bera Test, the result of Anderson-Darling suggests that the distribution
It is normal for waste.
Regarding the other test, the Kolmogorov-Smirnov test, this is implemented with the function
kstest_normal(...) of the module diagnostic of [Link]. The structure to be used for presenting -
The method for obtaining the results is the same as that used in the other statistical tests that have been carried out:

In [32]: TestResults=[Link].kstest_normal(Residues)
Nombres=["Estadístico KS","p-value"]
[Link](TestResults, index=Names)

KS Statistic 0.081275
p-value 0.200000
dtype: float64

Like with Jarque-Bera and Anderson-Darling, the results of the Kolmogorov-Smirnov test indicate
that the null hypothesis is not rejected and, therefore, it is not rejected that the distribution of the error term is
normal.
The last two tests performed are based on the Empirical Distribution Function, which, according to
Vaart (1998) is the natural estimator of the underlying distribution, F, if it is not known. To visualize
In Python, we can use Matplotlib functions to obtain the empirical distribution values.
corresponding to that function, it is enough to perform a few simple operations.
Vaart (1998) defines, for a random sample X1 , ...,X n the Empirical Distribution Function of the si-
next way:

1 n
n i∑
Fn (t) = 1 { X i≤ t }
=1
From this definition, it is very simple to determine the points with which the Function will be calculated.
Empirical Distribution. First, we must organize the values of the variable of interest from lowest to
mayor, for what the function sort(...) of NumPy is used:

In [33]: x = [Link](Residuos)

Each of these values corresponds to a 'probability' of being less than or equal; that is, if we have the
set {1,2,2,3,5,7,8,9,10,10} the probability that an element is less than or equal to 2 is 30% (since three
of the ten elements are less than or equal to two), the probability that an element is less than or equal to 7
It is 60% (since six of the ten elements are less than or equal to seven). This logic is very easy to implement.
In the code: it is only necessary to determine the number of elements, n, and divide a sequence 1, ..., n by this.
number, so that the corresponding probabilities are obtained:

61
In [34]: n=[Link]
y=[Link](1, n+1)/n

The pairs are already obtained.( x,y ) to graph the Empirical Distribution Function, but this only results in
useful if we can compare it to some reference: a theoretical distribution. The normal distribution is constructed
with two parameters: the mean and the standard deviation. If data is generated from a normal distribution
with a mean equal to the average of the residuals and a standard deviation equal to that of these, such data serves to
graph the theoretical distribution (as it is known that it is indeed a normal distribution), which
it can be used to compare with the empirical distribution; if the behavior is very similar, visually it
suggests that the distribution of the residues is normal.
Let's see what was previously described in practice. To generate random values from a distribution...
Normal distribution can use the function normal(...) from the random module of NumPy. This function must be
specify the mean (loc), the standard deviation (scale), and the number of values to generate (size):

In [35]: Theoretical_residuals=[Link](loc=[Link](Residuals),
scale=[Link](Residuos),
size= 10000)

From these theoretical residues, the 'theoretical' distribution function can be defined, for which it is enough
employ the same operations used in the empirical distribution:

In [36]: x_teor = [Link](Residuos_teoricos)

n_teor = x_teor.size
y_teor=[Link](1, n_teor+1)/n_teor

Finally, with the function subplots(...) from [Link] and the methods .scatter() and .plot(),
one can visualize the relationship between the previously defined functions:

In [37]: fig, ax = [Link](figsize=(8,6))

empirical distribution vs.\n"
"Distribución teórica",
fontsize= 18,
fontweight="bold"
[Link](x, y, color="darkblue",
label="Distribución empírica")
[Link](x_teor, y_teor, color="red",
label="Distribución teórica")
ax.set_xlabel(r"Values $x$")
fontsize= 15)
ax.set_ylabel("$P(X \leq x)$", fontsize= 15)
ax.tick_params(labelsize= 14)
[Link]()
[Link](.9,-.02,
Preparation:
fontsize= 13, fontweight="bold"
ha="right")
[Link](.9,-.08,
Triana, F.\n(2019)
fontsize= 12, ha="right"

62
 plt.subplots_adjust(top= 0.85)
[Link]()

The reader may observe that the distribution of the residues does not differ much from the appearance that corresponds.
to the normal distribution. The points corresponding to the empirical distribution function (the blue ones) are situated
very close to the reference line (the theoretical distribution).
Finally, another widely used graphic method, certainly more common than that of the distribution
Empirical, to assess normality, is the quantile-quantile plot, or Q-Q Plot. This plot can be constructed
with the function qqplot(...) from Statsmodels, passing as an argument the object that contains the data from
from which the quantiles are to be established. By default, this function works with the normal distribution,
reason why no modification should be introduced:

In [38]:fig, ax=[Link](figsize=(8,6))
Q-Q Plot
fontweight="bold"
[Link](Residuals, ax=ax, line="s",

63
 color="darkblue")
Set the label for the x-axis to 'Theoretical Quantiles', font size 14
Sample quantiles
[Link](.9,-.02,
"Elaboración:",
fontsize= 13, fontweight="bold",
ha="right")
[Link](.9,-.08,
Triana, F.\n(2019)
fontsize= 12, ha="right"
[Link]()

As can be seen, most points are located very close to the red reference line, indicating that
the distribution of the data does not differ much from the theoretical (the normal). This result agrees with the one obtained
in the analysis of the empirical distribution and with the results of the three statistical tests used for
verify the fulfillment of the assumption of normality.

64
 This chapter briefly addressed, without examining in detail the theoretical foundations.
given, the fulfillment of some of the assumptions of the classical linear regression model (normal, if considered
we solve the last assumption posed) using specialized Python library tools. The model in
which such assumptions correspond to was illustrated, by means of a particular example, in Chapter 3, in which it
I provided a slight practical approach to the fundamental topic that is generally addressed in higher-level courses.
introductory Econometrics. Now, the difficulty will be slightly increased by proceeding to examine a topic.
somewhat more complex (in a more theoretical than practical sense) to which the next chapter will be dedicated.

65
Chapter 5

Two-Stage Least Squares

One of the assumptions (specifically #2) of the classical linear regression model proposed in the early
The previous chapter establishes that the covariance between the explanatory variables and the error term is zero, that is,
that all (following the notation used previously) is independent of the error term. When
the explanatory variables fulfill this condition, they are called exogenous and when they violate such assumption they
endogenous.
When any of the explanatory variables of the model is (are) endogenous, the estimators obtained
by the Ordinary Least Squares method they lose their property of MELI. To solve such an inconvenience
the method of Instrumental Variables is used, in which variables not included in the equation are employed
original to estimate the included endogenous explanatory variables and address their endogeneity.
When one or more explanatory variables and the error term are correlated, conditions are violated.
of the assumptions of the classical linear regression model, so the OLS estimators lose some
of its attractive characteristics. To solve the issue generated by the endogeneity of some(s)
covariante(s) se requiere de información adicional, la cual se obtiene de una variable observable exógena no
included.
Taking a linear regression model of the form

y1i= β 0+ β 1 xi+ β 2 y2i+ ui

where1it is the dependent variable, x is an exogenous explanatory variable (it does not correlate with the
error term) yy2it is an endogenous explanatory variable (it correlates with the error term), if it is given
a variable that meets the following conditions:

[Link](z,u) = 0

[Link](z,y2 ) = 0

then you are an instrumental variable (or instrument) for y2 .

The first condition established refers to the exogeneity of the instrument and the second to the relevance.
advance of the same; these conditions are useful to guarantee the 'quality' of the Instrumental Variables, the
they are important as they allow addressing endogeneity and overcoming the problems caused by it.
For the interested reader who wishes to delve deeper into this topic, it is suggested to consult Chapter 15, Estimation with Variables.
instrumental variables and two-stage least squares, by Wooldridge (2010).
Up to this point, the reader should know that the correlation between explanatory variables and the error term
it is a problem for obtaining estimators with the best properties by the Least Squares method.
Ordinary drados and that instrumental variables can be used to address the endogeneity of the(s)

66
explanatory variable(s); however, the specific process by which it is done has not yet been indicated.
use of instruments.
One of the methods used to obtain Instrumental Variables estimators is the Least Squares method.
Squares in 2 Stages (MC2E), which receives this name because the process it encompasses is
structured in two well-defined 'stages' in which specific regressions are carried out with objectives
concretes.
Taking the linear regression model mentioned above, and following Wooldridge (2010), the pri-
my stage of the MC2E method consists of performing the regression ˆ 2I= π0ˆ+ π1 xiˆ+ π2 zi , where
ˆ is it ob-
they have the adjusted values and ˆ 2 The second stage is a regression by Ordinary Least Squares of
y1i= β 0+ β 1 xi+ β 2 y2i+ uˆi(note that in this regression, the values of y are used)2instead ofˆthe va-
lores day2 Thus, what it is basically about is estimating the endogenous explanatory variable from
exogenous variables and, subsequently, use the estimated variable as an explanatory variable, along with those of-
more exogenous than the original model (instruments are not included), to estimate the dependent variable.
A very minimal description of Instrumental Variables and the Least Squares method has already been given.
in 2 Stages, so the reader must have a general idea of the topic that is now intended to be addressed in a
practice; task to which, hereinafter, effective start is given with Python.
The dataset that will be used to illustrate the use of instrumental variables and the MC2E method is the one employed
in Example 15.5, Returns of education for the working woman, by Wooldridge (2010). The dataset is
has been obtained from Stata, with the command bcuse, and has been exported as a .csv file for use in Python.
The reader is invited to consult the example, so that they can verify the results obtained.

Environment preparation
The first cell in which the user will write and execute code will have the following content:

In [1]:import numpy as np
import pandas as pd
[Link]
from [Link] importIV2SLS
import [Link] as plt
%matplotlibinline
[Link]("seaborn-white")

The reader should already know that the previous block of code allows loading the necessary tools to
carry out the required actions in the proposed econometric analysis; in case there is no clarity
In this regard, you are asked to consult the section 'Preparing the Environment' of the second chapter, where it is provided
an explanation about such a question.
The code in this cell is exactly the same as that used in the environment preparation of Chapter 2.
except for the following: 1) the Seaborn library is not imported, 2) even more importantly, the function is imported
IV2SLS(...) from the iv module of the linearmodels library; this function is fundamental for the topic to be addressed in
this chapter and its importation is mandatory.
The user must be very careful when executing the previous code cell, as doing so immediately,
Without taking into account the information presented below, the process will not be completely successful;
this, for the reason that is presented next. The library linemodels is not integrated by default into the Distri-
Anaconda distribution (which is the one we are using), so it is necessary to install it on our own.
account to be able to use the tools it provides.
To install a library, simply go to the Command Prompt ('cmd'); once there, you just need to
write conda install and then the name of the library you are interested in. In our

67
In this case, the code to be used hides 'install linearmodels'. Alternatively, '!conda' can be used.
install linearmodels directly in the notebook.

Note: The above code uses conda, the package manager of the Anaconda distribution,
which relies on the Anaconda repository; however, it is possible that the library installation process
through decondano is successful, in which case the package management system can be used
Python standard, pip, which works directly with the Python Package Index (PyPi). The code for the
installation delinearmodels through pip pip install linearmodels in 'cmd' or !pip install
linear models the notebook; running this code should be successfully carried out the ins-
Installation of the library linear models.

The importation of model outlines is fundamental to our work, as it is where the fun-
actions that we will use to make use of instrumental variables and apply the method of Minimum Squares
two in 2 stages; without this library, the estimation process and execution would be exceedingly complex
of the tasks required to carry out an 'acceptable' econometric analysis. The user must verify the
prior installation delinearmodels (if not installed, carry out the installation process)
pending, be in accordance with pip) and, in case of reusing the code block for environment preparation
from the second chapter, you must be careful and ensure that you correctly include the import line of the
function IV2SLS(...) from the module ivdelinearmodels.
Once the code cell for environment setup has been executed without generating any type of error,
one can proceed with complete peace of mind to the next phase, which is fundamentally important in
how much is needed to have the 'raw material' to carry out the analysis: the data.

Data import
Data import is carried out using functions from the pandas library. Depending on the type of
file in which the dataset information is stored, we must resort to a partition function
cellular and a specific set of parameters.
For our specific case, the following code is the one that imports the dataset:

In [2]: data=pd.read_csv('[Link]')

To the reader who does not understand this line of code, it is suggested to refer to Chapter 2, where they will find a
explanation that will surely clarify your doubts.
Once the line of code in this cell has been executed, the dataset should have been imported correctly.
To verify this, we use the method .head(), applied to the object containing the dataset (in this case it
we have assigned the name data):

In [3]:[Link]()

Out[3]: inlf hours kidslt6 kidsge6 age educ wage repwage hushrs husage \
0 1 1610 1 0 32 12 3.3540 2.65 2708 34
1 1 1656 0 2 30 12 1.3889 2.65 2310 30
2 1 1980 1 3 35 12 4.5455 4.04 3072 40
3 1 456 0 3 34 12 1.0965 3.25 1920 53
4 1 1568 1 2 31 14 4.5918 3.60 2000 32

… faminc MTR mother education father education unemployment city experience

nwifeinc

68
 0 … 16310 0.7215 12 7 5.0 0 14 10.910060
1 … 21800 0.6615 7 7 11.0 1 5 19.499981
2 … 21040 0.6915 12 7 5.0 0 15 12.039910
3 … 7300 0.7815 7 7 5.0 0 6 6.799996
4 … 27300 0.6215 12 14 9.5 1 7 20.100060

wage expert
0 1.210154 196
1 0.328512 25
2 1.514138 225
3 0.092123 36
4 1.524272 49

[5 rows x 22 columns]

We can observe that the import process has been successful and the information from the dataset has been
correctly stored in a pandas DataFrame. Now, we can continue calmly with the development-
We are rolling out our analysis, as we have successfully completed the first step. We will proceed to use
instrumental variables using the Two-Stage Least Squares method, however, the user
You must take into account the information provided to you below.
The dataset we just imported contains some missing observations (this is due to the source
original of the information and not to a preference of the author of this work, although we will see how such a situation leads us
will result in something enriching by improving our exploratory and data management capabilities). These values
Missing items can create problems later on, so we will proceed to identify them and address them.
appropriate treatment.
To examine how many rows have a field with a missing value and know its location, use
we remove the methods .isna() and .sum(); the first allows us to identify missing values (through a
boolean criterion) and the second allows for executing a sum that takes into account the type of generated data.
by the first method, it allows knowing the quantity of them.

In [4]:[Link]().sum()

Out[4]:inlf 0
hours 0
kidslt6 0
kidsge6 0
age 0
educ 0
wage 325
repwage 0
hushrs 0
usage 0
huseduc 0
huswage 0
faminc 0
mtr 0
mothereduc 0
fat education 0

69
 I am 0
city 0
exper 0
nwifeinc 0
lwage 325
expersq 0
dtype: int64

It can be observed that the only variables with missing values are 'wage' and 'lwage'.
with 325 data classified as N for each one. This situation may lead to certain problems, so
which, in an effort to minimize the possibility of inconveniences during the estimation phase, will only be used
the records in which there is a specific value (noNa) for the variable 'lwage' (which has the role of
dependent variable in the model). How can we achieve this goal, selecting only records
with non-missing values?
To obtain a DataFrame that only contains records where there are no missing values for
the variable 'lwage' will use selection through brackets '[ ]' and the method .notna() (applied
about the reference variable); this, for the purpose of only considering records that meet the
The specified condition. The obtained result will be assigned to the existing object data, with the aim of counting.
with a single dataset and avoid possible confusions. The code to execute is the following:
In [5]: data=data[data["lwage"].notna()]
Now, we proceed to verify that we have achieved what we set out to do in our DataFrame.
there are no missing values; again, we will use the methods .isna() and .sum():

In [6]:[Link]().sum()
Out[6]:inlf 0
hours 0
kidslt6 0
kidsge6 0
age 0
educ 0
wage 0
repwage 0
hushrs 0
husage 0
useduc 0
huswage 0
faminc 0
mtr 0
mothereduc 0
father education0
we are 0
city 0
exper 0
nwifeinc 0
lwage 0
expersq 0
dtype: int64

70
 We observe that all the variables contained in the DataFrame lack missing values, so
no record is incomplete. This DataFrame, which is named data, is the one we will use to
work with instrumental variables and apply the Two-Stage Least Squares method, which addresses
We will continue from two approaches: 1) manual process and 2) automatic process.

Manual process
The Two-Stage Least Squares method, as its name indicates, is composed of two stages.
of clearly defined estimates: the first consists of returning the endogenous explanatory variable(s)
against the instrument(s) and the exogenous explanatory variable(s); the second stage consists of using
the estimated variable in the first stage instead of the endogenous explanatory variable in the original model. To
reader, if the description seems confusing, it is suggested to refer to the first part of this chapter or consult
Chapter 15, Estimation with Instrumental Variables and Two-Stage Least Squares, from Wooldridge (2010),
where this topic is addressed in detail.
The 'manual' process is named as such because we must carry out the regressions of the first
and the second stage on our own, specifying the set of variables to be used in each one.
of these. Such a procedure is not recommended and most programs and languages that support analysis
statistical (Python is no exception) has functions that perform the execution "automatically",
preventing the user from having to carry out the tasks on their own.
Despite what is stated in the previous paragraph, given the practical approach of this guide and the usefulness it can
to have for the user to carry out the process step by step in the understanding of the topic addressed, we will proceed
to apply the MC2E method 'manually' (the automatic process is addressed in the following section, so
to the reader who is not interested in carrying out both stages on their own, it is suggested to consult that
section directly, without dedicating their time to the study of the present.
The first thing we must do to correctly carry out the first stage regressions and the
The second stage is to appropriately group the available variables, regarding each of them.
plays a specific role. Thus, we will proceed to assign the dependent variable to a specific object, the variable
endogenous explanatory to another particular object and, likewise, the exogenous explanatory variables and the
instruments to other concrete objects. The code to execute is the following:

In [7]: Y1=data["lwage"]Dependent variable

Y2=data["educ"]Endogenous variable
X=data[["exper","expersq"]]Exogenous variables
Z=data[["motheduc","fatheduc"]]Instruments

First stage
During the first stage, the endogenous explanatory variable is regressed against the exogenous explanatory variables.
nas and the instruments. The reader can observe that the exogenous explanatory variables and the instruments have
have been assigned to different objects, therefore, in order to simplify the code, they will now be assigned
together into a single object. To achieve this goal, it is enough to use the function concat(...) of
pandas, using a Python list '[...]' as an argument that includes the DataFrames that contain
the exogenous explanatory variables and the instruments, assigning the value of 1 to the parameter axis (for conca-
to join the DataFrames 'horizontally'). The result generated by applying the function to the list
it will be called (something lengthy, but very informative) ofFirstStageRegressors, like this:

In [8]:frames=[X, Z]
RegressorsFirstStage=[Link](frames, axis= 1)

71
 We verify that the process has been successfully completed by using the [Link]():

In [9]:[Link]()

Out[9]: exper expersq motheduc fatheduc

0 14 196 12 7
1 5 25 7 7
2 15 225 12 7
3 6 36 7 7
4 7 49 12 14

Now, given that our variables are stored conveniently, we build the co- model
corresponding to the first stage, taking the endogenous explanatory variable ('educ') as the dependent variable
in original notation, Y2 in our environment) and as explanatory variables the exogenous explanatory variables
nas ('exper', 'expersq' in original notation) and the instruments ('motheduc', 'fatheduc' in original notation
Note that the exogenous explanatory variables and the instruments are contained in a single object.
(RegressionFirstStage)
The code to execute is the following:

In [10]: FirstStage=[Link](Y2, sm.add_constant(FirstStageRegressors))

ResultadosPrimeraEtapa=[Link]()
print([Link]())

OLS Regression Results

==============================================================================
Dependent Variable: educ R-squared: 0.211
Model: OLS Adj. R-squared: 0.204
Method: Least Squares F-statistic: 28.36
Date: xxx, xx xxx xxxx Prob (F-statistic): 6.87e-21
Time: xx:xx:xx Log-Likelihood: -909.72
No. Observations: 428 AIC: 1829.
Df Residuals: 423 BIC: 1850.
Df Model: 4
Covariance Type: nonrobust
==============================================================================
coef standard error t P>|t| 0.025 0.975
------------------------------------------------------------------------------
const 9.1026 0.427 21.340 0.000 8.264 9.941
exper 0.0452 0.040 1.124 0.262 -0.034 0.124
expert -0.0010 0.001 -0.839 0.402 -0.003 0.001
mothereduc 0.1576 0.036 4.391 0.000 0.087 0.228
fatheduc 0.1895 0.034 5.615 0.000 0.123 0.256
==============================================================================
Omnibus: 10.903 Durbin-Watson: 1.940
Prob(Omnibus): 0.004 Jarque-Bera (JB): 20.371
Skew: -0.013 Prob(JB): 0.0000377
Kurtosis: 4.068 Cond. No. 1550
==============================================================================

72
Warnings:
Standard Errors assume that the covariance matrix of the errors is correctly
specified.
The condition number is large, 1.55e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

To the reader who has closely followed this work, the previous code cell should not seem strange.
it's just an Ordinary Least Squares regression; if the provided code seems confusing to you,
It is suggested to consult the third chapter of this work, where the regression using OLS is the central topic.
We have used instruments such as covariates in a regression, however, we still do not know if di-
which instruments meet the conditions to satisfactorily perform such a role; we have not yet
evaluated the 'quality' of the instrumental variables used, so we will proceed to examine it.
In order for the instrumental variables used to meet the relevance condition (see the first part
from this chapter), the parameters associated with the instruments must be statistically different from zero.
As can be seen in the model's results instance, both for the variable 'motheduc' and for
the variable 'fatheduc' meets this condition, since the associated coefficients are statistically significant
The p-value is less than 0.05 (chosen significance level) and the confidence interval does not...
contains 0.
The conclusion of the previous paragraph is that the instruments meet the condition of relevance, however,
to formally test the validity of this conclusion, an F test for significance can be used
subset of parameters. The F-test is executed in Python through the method .f_test(), applied on
the instance of model results and using as argument the conditions to evaluate (which can be
to express oneself as a 'string'.
In this case, we will assign the hypothesis to an object that will be used as an argument in the method.
.f_test(); this, with the purpose of simplifying the code. Thus:

In [11]:hipotesis="(fatheduc = 0), (motheduc = 0)"

Test_F=FirstStageResults.f_test(hypothesis)
print(Test_F)

<F test: F=array([[55.40030043]]), p=4.268908724630835e-22, df_denom=423, df_num=2>

The result generated by the method .f_test() is an object of type ContrastResults(), which has
with their own attributes. Some of these attributes correspond to the statistic and the p-value, which are
the values of interest for our purposes. These attributes are NumPy ndarray objects, without
embargo, what interests us is not the object as such but its content, which must be expressed as a data
of tipofloat.
To extract the corresponding values, we resort to the use of selection through
brackets "[ ]" and the method .item(), indicating the position corresponding to the element of interest; likewise,
we will use a pandas Series to obtain a somewhat more 'pleasant' presentation of the result. The
the code to be executed is the following:

In [12]:Statistic=Test_F.statistic[0].item(0)
pvalue=Prueba_F.[Link](0)
Nombres=["Estadístico","p-value"]
[Link]([Statistic, round(pvalue, 3)], index=Names)

73
Statistic 55.4003
p-value 0.0000
dtype: float64

Assuming a significance level of 5%, how− value< α then the null hypothesis is rejected and
It concludes that the parameters associated with the instruments are statistically different from zero (working ...
with aα p < 0.05). The reader is invited to compare these results with those presented in Wooldridge (2010),
so that I can verify that they are exactly the same and that the developed process is correct.
Up to this point, the regression of the first stage has already been carried out and it has been verified that the instructions...
employees meet the relevance condition, so it is possible to proceed to the second stage
of the MC2E method, which is addressed in the following section.

Second stage
During the second stage, the dependent variable of the original equation is returned against the ex- variables.
exogenous explanatory variables and the variable estimated in the first stage (which replaces the endogenous explanatory variable
that initially caused the inconveniences) (Note that at this stage no instruments are used).
During the first stage, we estimate the endogenous explanatory variable as a function of the explanatory variables.
exogenous and the instruments; however, we still do not know the values of this variable, which are the
to be used in the second stage. To obtain such values, we use the method .predict() applied
to the results instance of the first stage model (which we have conveniently named
ResultsFirstStage).
The results generated by the application of the method .predict() will be assigned to a new variable.
We will call it predicted_educ from the DataFrame data, so that we can easily use it.
subsequently. The code to be executed is the following:

In [13]:data["predicted_educ"]=[Link]()

We verify that the new variable has indeed been created by using the method .head():

In [14]:[Link]()

Out[14]: inlf hours kidslt6 kidsge6 age educ wage repwage hushrs husage
0 1 1610 1 0 32 12 3.3540 2.65 2708 34
1 1 1656 0 2 30 12 1.3889 2.65 2310 30
2 1 1980 1 3 35 12 4.5455 4.04 3072 40
3 1 456 0 3 34 12 1.0965 3.25 1920 53
4 1 1568 1 2 31 14 4.5918 3.60 2000 32

… meter mother educ father educ unemployment city experience nwifeinc

0 … 0.7215 12 7 5.0 0 14 10.910060
1 … 0.6615 7 7 11.0 1 5 19.499981
2 … 0.6915 12 7 5.0 0 15 12.039910
3 … 0.7815 7 7 5.0 0 6 6.799996
4 … 0.6215 12 14 9.5 1 7 20.100060

wage expert predicted_educ

0 1.210154 196 12.756017
1 0.328512 25 11.733558

74
 2 1.514138 225 12.771979
3 0.092123 36 11.767683
4 1.524272 49 13.914615

[5 rows x 23 columns]

We can observe that the last column corresponds to the variable predicted_educ, which is the one that
we have created previously; the data it contains corresponds to the estimated values during the first
stage.
We now have all the variables to be used in the second stage regression. Now, just like
during the first stage, we will group the corresponding regressors into a single object. Again, create-
we will create a Python list "[...]" whose content will correspond to the regressors of the second stage, that is, the
exogenous explanatory variables of the original model and the endogenous variable estimated in the first stage; this
The list will be used as an argument of the function concat(...) from pandas (remember the assignment of the value of
1 to the parameter axis). Thus:

In [15]:FramesSegundaEtapa=[X, data["predicted_educ"]]
RegresorasSecondStage=[Link](FramesSecondStage, axis= 1)

We verify that the process has been successfully carried out by resorting to the [Link]():

In [16]:[Link]()

Out[16]: exper expersq predicted_educ

0 14 196 12.756017
1 5 25 11.733558
2 15 225 12.771979
3 6 36 11.767683
4 7 49 13.914615

Now, as we have repeatedly pointed out, we build the model using the dependent variable...
the dependent variable of the original equation (y1 , following the notation used at the beginning of the chapter)
and as explanatory variables the exogenous explanatory variables of the original equation (lasx, following the no-
initial estimation) and the variable estimated in the first stage (andˆ 2and not2 !).
Since we have grouped the explanatory variables of the second stage regression into a single object
(it is what we have done in the last code cell), it is enough to use that object as an argument in the
construction of the model (which is carried out with the already known function OLS(...) from Statsmodels).
The code (which should already be familiar) for the construction and estimation of the second stage model
es el siguiente:

In [17]: SecondStage=[Link](Y1, sm.add_constant(RegressorsSecondStage))

ResultadosSegundaEtapa=[Link]()
print([Link]())

OLS Regression Results

==============================================================================
Dep. Variable: lwage R-squared: 0.050
Model: OLS Adj. R-squared: 0.043
Method: Least Squares F-statistic: 7.405

75
Date: xxx, xx xxx xxxx Prob (F-statistic): 0.0000762
Time: xx:xx:xx Log-Likelihood: -457.17
No. Observations: 428 AIC: 922.3
Df Residuals: 424 BIC: 938.6
Df Model: 3
Covariance Type: nonrobust
==================================================================================
coef standard error t P>|t| 0.025 0.975
----------------------------------------------------------------------------------
const 0.0481 0.420 0.115 0.909 -0.777 0.873
exper 0.0442 0.014 3.136 0.002 0.016 0.072
expersq -0.0009 0.000 -2.134 0.033 -0.002 -7.11e-05
predicted_educ 0.0614 0.033 1.863 0.063 -0.003 0.126
==============================================================================
Omnibus: 53.587 Durbin-Watson: 1.959
Prob(Omnibus): 0.000 Jarque-Bera (JB): 168.354
Skew: -0.551 Prob(JB): 2.77e-37
Kurtosis: 5.868 Cond. No. 4410
==============================================================================

Warnings:
Standard Errors assume that the covariance matrix of the errors is correctly
specified.
The condition number is large, 4.41e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

Thus concludes the manual application of the MC2E method in Python. The reader should take into account
that although the parameters estimated "manually" are correct, the standard errors are not (see the
page 522 of Wooldridge (2010)), so it is not advisable to run the first and second regressions
from the stage of self-employment. To avoid this inconvenience, most packages and statistical programs
they offer the possibility to execute the MC2E method automatically, generating correct results. The
The Python programming language is no exception, as its library linear models offers a function designed
specifically for such task, as we will see below.

Automatic process
To avoid having to perform the two-stage regressions of the MC2E method on our own,
the library of Python linear models gives us the possibility to use its functions to perform 'automatic-
I said regressions. The specific function used for this purpose is IV2SLS... (acronym for
Instrumental Variables 2 Stages Least Squares), whose arguments include the dependent variable (dependent), the
exogenous regressors (exog), endogenous regressors (endog) and instruments (instruments).
For estimation, as is familiar, the method .fit() is used, and for the summary of results the
[Link]
The code corresponding to our example is:

In [18]:TwoStageModel=IV2SLS(dependent=Y1, exog=sm.add_constant(X),
endog=Y2, instruments=Z

76
 ResultadosModelo2Etapas=[Link]()
print([Link])

IV-2SLS Estimation Summary

==============================================================================
Dependent Variable: lwage R-squared: 0.1357
Estimator: IV-2SLS Adj. R-squared: 0.1296
No. Observations: 428 F-statistic: 18.611
Date: xxx, xxx xx xxxx P-value (F-stat) 0.0003
Time: xx:xx:xx Distribution: chi2(3)
Cov. Estimator: robust

Parameter Estimates
==============================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
------------------------------------------------------------------------------
const 0.0481 0.4278 0.1124 0.9105 -0.7903 0.8865
exper 0.0442 0.0155 2.8546 0.0043 0.0138 0.0745
experts -0.0009 0.0004 -2.1001 0.0357 -0.0017 -5.997e-05
educ 0.0614 0.0332 1.8503 0.0643 -0.0036 0.1264
==============================================================================

Endogenous: education
Instruments: motheduc, fatheduc
Robust Covariance (Heteroskedastic)
Debiased: False

The reader can observe how the application of the MC2E method is tremendously simple with the use of
the function IV2SLS(...) [Link] that all tasks carried out in the manual process have been
executed perfectly in just 3 lines of code!, without having to create additional variables or use
multiple groups of the same, obtaining, in addition, completely correct results.

Ordinary Least Squares vs. Two-Stage Least Squares

Finally, a common ordinary least squares regression will be performed, without
Keep in mind that the endogenous regressor is endogenous unless considering all regressors as exogenous.
This will be done solely for illustrative purposes, to examine what the outcomes would have been.
Again, we will use the concat(...) function from pandas to group the explanatory variables.
(exogenous and endogenous) in a single DataFrame. Likewise, in order not to unnecessarily extend ourselves, we will take
the creation and estimation of the model is completed in the same code cell:

In [19]: FramesMCO = [X, Y2]

BacktrackersMCO=[Link](FramesMCO, axis= 1)
ModeloMCO=[Link](Y1, sm.add_constant(RegresorasMCO))
ResultadosMCO=[Link](cov_type="HC0")
print([Link]())

77
 OLS Regression Results
==============================================================================
Dependent Variable: lwage R-squared: 0.157
Model: OLS Adj. R-squared: 0.151
Method: Least Squares F-statistic: 27.56
Date: xxx, xx xxx xxxx Prob (F-statistic): 2.68e-16
Time: xx:xx:xx Log-Likelihood: -431.60
No. Observations: 428 AIC: 871.2
Df Residuals: 424 BIC: 887.4
Df Model: 3
Covariance Type: HC0
==============================================================================
coef standard error z P>|z| 0.025 0.975
------------------------------------------------------------------------------
const -0.5220 0.201 -2.601 0.009 -0.915 -0.129
exper 0.0416 0.015 2.734 0.006 0.012 0.071
expersq -0.0008 0.000 -1.940 0.052 -0.002 8.28e-06
educ 0.1075 0.013 8.170 0.000 0.082 0.133
==============================================================================
Omnibus: 77.792 Durbin-Watson: 1.961
Prob(Omnibus): 0.000 Jarque-Bera (JB): 300.917
Skew: -0.753 Prob(JB): 4.54e-66
Kurtosis: 6.822 Cond. No. 2200
==============================================================================

Warnings:
Standard Errors are heteroscedasticity robust (HC0)
The condition number is large, 2.21e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

The reader can observe that the results obtained by the OLS method differ to some extent from the
corresponding to the MC2E method and that the standard errors estimated by MC2E manually are different
to the automatically obtained ones. You can verify this by comparing the summary tables of each of the
elaborated models; however, such an exercise can be tedious, as each of these tables is
situated in a different position within the notebook and comparing information between them can be exhausting
Fortunately, the libreríalinearmodels provides a tool to facilitate this work, as we will see.
right away.

Comparison of models
It is very practical to be able to compare the estimated models by visualizing their corresponding summaries.
instances of results simultaneously; this possibility can materialize thanks to the function
compare(...) that the library provides in its module iv. The first thing we must do to use
the function compare(...) is to have access to it, so we proceed to import it:

In [20]:from [Link] importcompare

78
 Having executed the previous code cell, we will have imported the function compare(...), so
that we will be able to make effective use of it. As an argument we will use a Python dictionary "{...}", in
we will use the names we assign to each model as keys and the instances of
results of the corresponding models.
The user must keep in mind that for the function compare(...) to outline [Link] to not generate a
error, all models to be compared must correspond to instances of results of the type generated by the
functionIV2SLS(...). In this case, since the functionIV2SLS(...) was only used for the last model, not
it is possible to make a comparison with others. How can such an inconvenience be solved? The answer is
very simple actually: all the models we have created can be built using the function
IV2SLS(...), you just need to be careful and include the correct arguments.
The function IV2SLS(...) has specific parameters for the dependent variable, the exogenous variable(s)
the endogenous variable(s) and the instrument(s). Thus, for a common ordinary estimation by the
MCO method, just assign None to the endogeneity parameter and the instruments parameter; for an estimation
por MC2E manual basta con emplear las regresoras de la segunda etapa en el parámetroexogy asignarNone
to the endogeinstruments parameters; for an automatic MC2E estimation, only the
the same arguments as in the model previously developed in the section 'Automatic process'. The code
(the reader is advised to examine it closely) to estimate the corresponding models, using for
The entire function IV2SLS(...), is the following:

In [21]: MCOResults=IV2SLS(dependent=Y1, exog=sm.add_constant(MCORegressors),

endog=None, instruments=None).fit()
ResultadosManual=IV2SLS(dependent=Y1,
exog=sm.add_constant(SecondStageRegressors),
endog=None, instruments=None).fit()
ResultadosAutomatico=IV2SLS(dependent=Y1, exog=sm.add_constant(X),
endog=Y2, instruments=Z).fit()

Now, since all the result instances have been generated by a model of the same class, po-
We can use the function compare(...) to compare the results obtained and contrast the models.
like this:

In [22]:print(compare({"MCO": ResultadosMCO,
"2 Etapas Manual": ResultadosManual,
2 Automatic Stages: AutomaticResults

Model Comparison
======================================================================
2 Automatic Stages 2 Manual Stages MCO
----------------------------------------------------------------------
Dep. Variable lwage wage lwage
Estimator IV-2SLS OLS OLS
No. Observations 428 428 428
Cov. Est. robust robust robust
R-squared 0.1357 0.0498 0.1568
Adj. R-squared 0.1296 0.0431 0.1509
F-statistic 18.611 17.111 82.671
P-value (F-stat) 0.0003 0.0007 0.0000
================== =========== =========== ===========
const 0.0481 0.0481 -0.5220

79
 (0.1124) (0.1071) -2.6010
expert 0.0442 0.0442 0.0416
(2.8546) (2.7045) (2.7344)
expersq -0.0009 -0.0009 -0.0008
-2.1001 -1.9627 (-1.9402)
educ 0.0614 0.1075
(1.8503) (8.1697)
predicted_educ 0.0614
(1.7553)
==================== ============= ============= =============
Instruments mothereduc
fat education
----------------------------------------------------------------------

T-stats reported in parentheses

The output generated by the execution of the code in the previous cell allows us to visualize simultaneously
see the results of the different models built, facilitating their comparison. The reader
you can see that at the bottom of the output it indicates T-stats reported in parentheses, which
indicates that the reported number in parentheses, below each of the coefficients, corresponds to the value
However, on multiple occasions what is reported with respect to a coefficient is not the value.
but the standard error. Is it possible to report the standard error instead of the value? The answer is yes: the function
compare(...) has a parameter precision, default value tstats, which allows specifying the estimator
of precision (within the available ones) to include in the output.
To report standard errors instead of values, simply assign 'std_errors' to the precision parameter.
from the function compare(...), as evidenced below:

In [23]: print(compare({"MCO": ResultadosMCO,

"2 Etapas Manual": ResultadosManual,
"2 Etapas Automático": ResultadosAutomatico},
precision="std_errors"

Model Comparison
=====================================================================
2 Automatic Stages 2 Manual Stages MCO
---------------------------------------------------------------------
Dep. Variable lwage wage wage
Estimator IV-2SLS OLS OLS
No. Observations 428 428 428
Cov. Est. robust robust robust
R-squared 0.1357 0.0498 0.1568
Adj. R-squared 0.1296 0.0431 0.1509
F-statistic 18.611 17.111 82.671
P-value (F-stat) 0.0003 0.0007 0.0000
================== ========== ========== ==========
const 0.0481 0.0481 -0.5220
0.4278 (0.4492) (0.2007)

80
exper 0.0442 0.0442 0.0416
0.0155 (0.0163) (0.0152)
expert -0.0009 -0.0009 -0.0008
0.0004 0.0005 0.0004
educate 0.0614 0.1075
0.0332 (0.0132)
predicted_educ 0.0614
0.0350
==================== ============ ============ ============
Instruments educational mother
fatheduc
---------------------------------------------------------------------

Std. Errors reported in parentheses

Note that the values reported in parentheses have changed and that in the lower section of the output
Standard errors reported in parentheses instead of the T-stats reported in parentheses.
Preliminary case. In this way, what is now reported in correspondence with each coefficient is the standard error.
and not its value.
At this point, the topic of instrumental variables and the Two-Stage Least Squares method concludes.
In Python. The reader will notice that, so far, only dependent variables have been worked with.
you continue and the regressions carried out in the different examples are linear; this will change slightly
in the next chapter, where there will be a very brief exposure of another type of models.

81
Chapter 6

Discrete dependent variable - Binary case

In the previous chapters, linear regression models have been worked on using a continuous variable.
as a dependent variable; now, we will step away a bit from models with such characteristics, exploring
nonlinear regressions with discrete dependent variables.
In the first models we will address, the dependent variable is a categorical variable, with only two
categories, coded with the values 0 and 1. This variable only takes these two values, which do not have a
signified by oneself (in contrast to the case of continuous variables), as far as its decoding functions:
it takes the value of 1 when the observation belongs to a specific group and 0 when it does not belong to such a group.
In somewhat simpler words: the variable takes the value of 1 when a condition is met (membership
to a specific classification) and 0 when this condition is not met.
It is possible that the reader does not have complete clarity about the characteristics of the variables they use.
we hope that the previous description was not entirely satisfactory; in that case, it may be that some
examples be helpful:

If we have a group made up of people of both sexes, there will be men and ha-
brá women; if we use a 'female' variable to perform a classification by sex, this variable
it will take the value of 1 when the information corresponds to a woman and 0 when it corresponds to a
man. 'female' is a binary variable, as it takes a value when a condition is met (1 when
the condition of female sex is met) and takes another value when the condition is not met (0
when the condition of being female is not met.

If we have a group made up of men, it is possible that some of these men are married.
while others do not, thus, we can use a binary variable 'married' to classify them in the group
corresponding. The variable 'married' takes the value of 1 when it refers to a married man and takes
the value of 0 when it corresponds to an unmarried man.

En un grupo de estudiantes universitarios algunos de estos habrán perdido asignaturas mientras que
others do not. In this way, we can use a binary variable that takes the value of 1 when it corresponds
ponda a un estudiante que ha perdido asignaturas y de 0 cuando corresponda a un estudiante que no ha
lost subjects.

The reader can observe that the binary variable takes the value of 1 when the specific condition is met.
and 0 otherwise; such assignment of values makes no numerical sense and does not affect the regression. Calm down-
It is possible to assign the value of 0 when the condition is met and 1 otherwise; however, for
convention, and greater ease of interpretation, 1 is used for the fulfillment of the condition and 0 for
the non-compliance with it.

82
 Up to this point, the reader has a general idea of what a binary variable is; however, it has not yet been explained.
He has informed how that variable will be used in the regression. It is very important to note that in the models
What we will address next is what is actually estimated is the probability that the dependent variable
take the value of 1 (to the reader interested in knowing why, it is recommended to consult Chapter 15, Models
of qualitative response regression, by Gujarati and Porter (2010); this probability is given by the parameters
(β) and the explanatory variables (x).
In somewhat more technical terms:

Prob[yi= 1| xi ] = P= Fi ( x Tβ ) i

in dondexTi corresponds to the transposed vector of explanatory variables andβ to the parameter vector. Thus,
the coefficients obtained during estimation are associated with the impact that a certain variable has
explanatory about the probability that the dependent variable takes the value of 1; in other words, the coefficients
clients are related (not necessarily quantified directly, due to the functional form adopted,
as we will see later) to the variation in the probability that the dependent variable takes the value of
1.
Some of the models used to work with binary response variables are: 1) Linear Model
of Probability, 2) Logit Model and 3) Probit Model. It is possible to easily construct and estimate each one of them.
These models in Python with the help of specialized libraries. Such tasks are what we will carry out.
continuation, but not before making a very brief explanation of the ideas behind each of the models
treaties.
The dataset we will use to develop the various models will be the one used in Example 15.7.
To smoke or not to smoke, by Gujarati and Porter (2010). This set of information has been obtained from the eco website.
[Link] of SHAZAM Analytics Ltd.

Linear probability model

The Linear Probability Model (or LPM) is, using a non-technical definition, basically a mo-
regression work in which the dependent variable is a binary variable and the estimation is carried out by the
Ordinary Least Squares Method.
Following the notation we have used throughout the work, in the MLP we have:

Prob[yi= 1] = β 0+ β 1 x1i+ ... + β k xwhat+ ui

where the estimators of theβ are obtained by the Ordinary Least Squares method.
Since it is a common ordinary least squares regression, estimating a MLP in Python does not
It should seem strange to the reader, as it is enough to use the same functions employed for the re- model.
multiple linear regression addressed in Chapter 3. In fact, the binary response variable models that are
They address in this chapter should not pose a greater inconvenience to the user regarding their construction.
it is relatively similar to the linear regression model and does not require resorting to different libraries than the
that we have used in previous chapters.

Preparation of the environment

The first cell where the user will write and execute code will have the following content:

In [1]:import numpy as np
import pandas as pd

83
 import [Link] as sm
[Link]
import [Link] as plt

The previous block of code allows loading the necessary tools to carry out the actions.
required in the proposed econometric analysis; in case the reader is not clear about this, they are
please refer to the section 'Preparing the environment' of Chapter 2, where an explanation is provided about
that issue.

Data import
The import of the data is carried out with functions from the pandas library. Depending on the type of
file in which the dataset information is stored, we will need to resort to a partitioning function-
cellular and a specific set of parameters.
The dataset we will import comes from a website, has no header and is space-separated.
therefore we must carefully specify each of the required arguments. Specifically, the
The following code is what should be executed to import the dataset:
In [2]:data=pd.read_csv("[Link]
header=None,
delim_whitespace=True,
names=["Fumador","Edad","Escolaridad","Ingreso",
Pcigs79

To the reader who is not familiar with the data import process, it is suggested to refer to Chapter
2, where you will find an explanation that will surely clarify your doubts. To the reader who is familiar with
With this process, it is noteworthy that, unlike other occasions, additional arguments are used, requiring
Given the characteristics of the original dataset: the header parameter allows specifying whether or not to include a
header; the parameter delim_whitespace allows indicating if the delimiter is whitespace (it is a pa-
boolean metro); the parameter names allows assigning names to the columns of the DataFrame through
from a Python list '[...]'.
Once the line of code in this cell has been executed, the dataset should have been imported correctly.
Otherwise, if the dataset cannot be imported directly from the website, it is possible to upload it from
the computer (as we have done in all previous chapters).
The code to be used in case the import process from the website fails is as follows:
In [3]: data=pd.read_csv("[Link]",
index_col= 0)

Having executed any of the previous code cells, the information of the dataset to be used...
it should have been imported correctly. To verify this, we resort to the [Link](), applied to the object
what the dataset contains (in this case we have assigned it the name data):

In [3]:[Link]()
Out[3]: Fumador Edad Escolaridad Ingreso Pcigs79
0 0 21 12.0 8500 60.6
1 1 28 15.0 12500 60.6
2 0 67 10.0 12500 60.6
3 0 20 12.0 12500 60.6
4 1 32 12.0 20000 60.6

84
 We can observe that the import process has been successful and the dataset information has been soul-
it was correctly centered in a pandas DataFrame. Now, we can continue calmly with the development
From our analysis, we have successfully completed the first step.

Data preparation
Before assigning the variables to the corresponding objects, the descriptive statistics will be examined, to
what is enough is to use the method .describe():

In [4]: [Link]()
Out[4]: Smoker School Age Income Pcigs79
count 1196.000000 1196.000000 1196.000000 1196.000000 1196.000000
mean 0.380435 41.806856 12.221154 19304.765886 60.984950
std 0.485697 17.056941 3.275847 9083.511331 4.848666
min 0.000000 17.000000 0.000000 500.000000 46.300000
25% 0.000000 27.000000 10.000000 12500.000000 59.225000
50% 0.000000 39.000000 12.000000 20000.000000 62.100000
75% 1.000000 56.000000 13.500000 30000.000000 63.800000
max 1.000000 88.000000 18.000000 30000.000000 69.800000
Focus your attention on the column corresponding to the dependent variable ('Smoker'). This variable is
binary and only takes the value of 0 or the value of 1; although its average is 0.38, there is no observation
for which the value of 'Smoker' is 0.38, in fact, there is no observation for which the value of this
a binary variable can be different from 0 or 1. Despite this, the average of a binary variable provides valuable information:
indicates a certain idea about the way the categories are distributed.
Since there are only zeros and ones, the average of this type of variable is equal to the sum of the ones divided by
by the total number of observations (n); therefore, the average of a binary variable corresponds to the
percentage of observations that meet the specific condition. For this specific example, that the
an average of 'Smoker' of 0.38 means that 38% of the examined individuals are smokers and, therefore,
The remaining 62% are not. To have greater clarity on this class distribution, a brief may be helpful.
appropriate visual examination result.
Since the dependent variable is a discrete variable with only two categories, all the data of that variable
will be grouped into only two points (0 and 1) and, according to descriptive statistics, the amount of
data corresponding to 1 will be slightly more than half of that corresponding to 0, so the bar for 1 will be
somewhat higher than 1/2 of the bar corresponding to 0. Let's examine this through a histogram, using
the function subplots(...) from [Link] and the method .hist():

In [5]:fig, ax=[Link]()
Figure title: 'Distribution of Smokers', size=17
fontweight="bold"
[Link](data["Smoker"], color="darkblue")
[Link](.9,-.02,
Preparation:
fontsize= 12, fontweight="bold"
right
[Link](.9,-.1,
Triana, F.\n(2019)
fontsize= 11, ha="right"
[Link]()

85
 The reader who has followed this work to this point will have noticed that in the previous chapters we have
storing various variables in specific objects. We have done this to group these variables of
according to your particular role in the regression, trying to simplify the code to be used later. Now,
For the different models that we will address in this chapter, the same variables will be used, grouped from the
same way for all. We will proceed to perform the grouping through the brackets "[ ]" and the lists
Python "[...]", as we have done in previous chapters:

In [6]: Y=data["Smoker"]
X=data[["Edad","Escolaridad","Ingreso","Pcigs79"]]

The first model to be examined is the Linear Probability Model, which is basically a model of
linear regression with a binary dependent variable and estimation by the Ordinary Least Squares method.
Linear regression and the OLS method are discussed in Chapter 3 of this work, so anyone who is not clear about it
In this regard, it is suggested to consult that chapter.
For the construction of the model and estimation of the coefficients, we use (as we have done in oca-
previous sessions) the OLS(...) function from Statsmodels and the .fit() method; for the visualization of the results,
the function print(...) and the method .summary(). The code to be executed is as follows:

In [7]: MLP=[Link](Y, sm.add_constant(X))

ResultadosMLP=[Link]()
print([Link]())

86
 OLS Regression Results
==============================================================================
Dependent variable: Smoker R-squared: 0.039
Model: OLS Adj. R-squared: 0.036
Method: Least Squares F-statistic: 12.01
Date: xxx, xx xxx xxxx Prob (F-statistic): 1.43e-09
Time: xx:xx:xx Log-Likelihood: -809.19
No. Observations: 1196 AIC: 1628.
Df Residuals: 1191 BIC: 1654.
Df Model: 4
Covariance Type: nonrobust
===============================================================================
coef standard error t P>|t| 0.025 0.975
-------------------------------------------------------------------------------
const 1.1231 0.188 5.963 0.000 0.754 1.493
Age -0.0047 0.001 -5.701 0.000 -0.006 -0.003
Education -0.0206 0.005 -4.465 0.000 -0.030 -0.012
Income 1.026e-06 1.63e-06 0.629 0.530 -2.18e-06 4.23e-06
Pcigs79 -0.0051 0.003 -1.799 0.072 -0.011 0.000
==============================================================================
Omnibus: 37.223 Durbin-Watson: 1.944
Prob(Omnibus): 0.000 Jarque-Bera (JB): 173.523
Skew: 0.449 Prob(JB): 2.09e-38
Kurtosis: 1.364 Cond. No. 291000
==============================================================================

Warnings:
Standard Errors assume that the covariance matrix of the errors is correctly
specified.
The condition number is large, 2.91e+05. This might indicate that there are
strong multicollinearity or other numerical problems.

One of the most important aspects to consider in binary dependent variable models is the
marginal effect of each of the explanatory variables; such effect corresponds to the variation in the probability
dad (that the dependent variable takes the value of 1) regarding the variation in an explanatory variable. In
xβiT)
somewhat more technical terms, the marginal effect of a variable xjcorresponds to the value of∂F(∂x , where
ji
x jit is a variable that is part of vectorx.
Thus, for the Linear Probability Model, the marginal effects are nothing more than the coefficients, because:

F ( xTi β) = β 0+ β 1 x1i+ ... + β k xwho

∂F ( xTi β)
= βj
∂x hi
Therefore, to know the marginal effects of the Linear Probability Model in Python, it is enough to
to know the coefficients associated with the variables. These coefficients correspond to one of the attributes (that
can be consulted with the function dir(...) of the results instance, and it is identified as .params.

87
We can assign the application of this attribute to a specific object and visualize its content using the
function print(...). Like this:

In [8]: MLPMargEff = [Link]

Marginal Effects MLP\n", MLPMargEff)

Marginal Effects MLP

const 1.123089
Age -0.004726
Education -0.020613
Income 0.000001
Pcigs79 -0.005132
dtype: float64

Thus, we observe, for example, that an increase of one year in age is associated, on average, with a
a decrease of 0.47 percentage points in the probability of being a smoker, and that an increase of 1 year in the
education is associated, on average, with a decrease of 2 percentage points in the probability of being
smoker.
In reality, for our purposes, there are no more aspects to discuss in detail in the Linear Model.
of Probability; it is enough to know its construction and estimation process and the obtaining of the effects mar-
Now, as mentioned at the beginning of the chapter, we will move away from linear regression models and
we will explore other possibilities: the first, the Logit Model.

Logit Model
The Linear Probability Model, although very simple, is not widely used in practice;
the reason for such a situation lies mainly in the fact that this model can lead to the
obtaining probabilities lacking meaning (for example, greater than 1 or that violate the axiom that the
probability cannot be negative.
Although it was not mentioned in the previous section, the MLP can be very problematic, as it can
generate probabilities less than 0 or greater than 1, which detracts from attractiveness and convenience. To overcome the
deficiencies of the Linear Probability Model resort to nonlinear regression models, of which the
the most popular are Logity and Probit.
In the Logit Model, the probability that the dependent variable takes the value of 1 is expressed as a
function of the following form:

1
Pi= FxTβ =i
1+ e−( β +0β x+1...+1iβ x) k who

Alternatively, taking Zi= β 0 + β 1 x1i+ ... + β k xwho , the function can be expressed as follows
way:

1
Pi= F ( Zi ) =
1+ e−Z i

Another way to write the previous expression is PI= e Zi 1 Thus, it is established that:
1+e Zi , therefore1− P= i 1+ e Z i

Pi
= eZ i
1− P i

88
 Applying the logarithm to the previous expression yields:

PI
L= ln = Zi= β 0+ β 1 x1i+ ... + β k xwho
1− P i
The denominality corresponds to the logarithm of the odds ratio. Thus, "while the MLP
it is assumed that
i Pest is linearly related to xi the logit model assumes that the logarithm of the ratio of
probabilities are linearly related to xi (Gujarati and Porter, 2010, p. 555).
The estimation of the parameters (βIn the Logit Model, it is carried out through the Maximum method.
Likelihood (Maximum Likelihood) or ML. To the reader interested in learning more about the Logit Model
and its estimation process, it is suggested to consult Chapter 15, Models of regression of qualitative response,
Gujarati and Porter (2010).
One of the points to consider in the Logit Model is that the interpretation of the coefficients does not result from
greater utility, as they quantify the impact that a change in an explanatory variable has on
logarithm of the odds ratio, something that, at first glance, has little illustrative value. Due to this, it
it relies on marginal effects to quantify the impact of the explanatory variables.
As is evident, the marginal effect of an explanatory variable in a Logit Model, unlike the
In the case of the Linear Probability Model, it does not correspond to the coefficient associated with that variable, but rather it is
given by:

∂FxTi β e Zi
= βj
∂x yes (1+ e ZI)2
Thus, it is observed that the marginal effect of an explanatory variable depends on the set of explanatory variables.
captives and is not given directly by the coefficient associated with it. Likewise, it is important to note that in the
the previous expression has found the marginal effect of variable xjfor observation, but it must be taken into account
consideration that the dataset contains information from multiple observations and the values of the variables are
differences between these observations, so there is no single value for the marginal effect of the variable
x jotherwise (number of observations) values for this.
Since it is not very practical to know the magnitude of the marginal effect of a variable for each one of
the observations, it is necessary to have a 'general' measure of the marginal effect of a certain variable,
But how is this 'general' measure obtained? One possible answer to this question is: there are several
"Paths." The most popular ones are presented below.
The first path (known as marginal effects at the mean) to obtain a 'general' measure of
The marginal effect of a variable is to use the average observation to calculate it; that is to say, to use the average.
of each variable as the value of the corresponding variable within Z.
Defining the vector that contains the average values of the explanatory variables as x̄, and having in
it counts that, consequently, x̄β = Z,T ¯
the marginal effect of the variable xj is given by:
¯
∂P eZ
∂x j 2
Z¯
= β j 1+ e
Thus, the 'general' marginal effect of a certain variable corresponds to the marginal effect of that variable.
calculated with the values corresponding to the average observation.
The other way (known as average marginal effects) to obtain a ‘general’ measure of the effect
The marginal of a variable is obtaining the average of the individual marginal effects corresponding to it.
variable; that is to say, to calculate the marginal effect of the variable xjfor each of the observations (in total it
hold observations), and find the average of those marginal effects. Specifically, the marginal effect
from the variablexj is given by:

89
 ∑ in=1
∂FxTi β
∑ in=1 e Zi
2
∂P ∂x yes 1+e Zi
= = βj
∂x j n n
Thus, the 'general' marginal effect of the determined variable corresponds to the average of the marginal effects.
individuals for said variable.
Up to this point, a very brief description has been given of the general idea that the Model is based on.
Logit; now, it's time to put into practice the knowledge acquired by carrying out the construction.
and estimation of a logit model in Python.

Logit Model in Python

The dataset that will be used for the construction and estimation of the Logit Model has been imported prior to
mind, in the section of this chapter that addresses the Linear Probability Model; likewise, the preparation
The environment and data have been successfully carried out before. Given that the construction and esti-
Logit model estimation is achieved with Statsmodels functions and does not require instruments from.
Apart from the libraries that have already been imported, all the tools are available to work adequately.
mind. Users who have not completed the exercise corresponding to the MLP section are invited to execute the
code cells #1, #2, and #6, which are essential for the work that is now intended to be developed.
To build a logit model in Python, the Logit(...) function from Statsmodels is used; its arguments
includes the dependent variable, the explanatory variables, and the missing parameter (with default value None)
which indicates the action to take in case there are missing values.
For the estimation, as in the models we have discussed in previous chapters, the method is used
.fit(); for the visualization of the summary of results the method .summary() (applied on the instance
of results) and the function print(...) (taking as an argument the application of the method .summary()).
The code corresponding to our case is:

In [9]:LogitModel=[Link](Y, sm.add_constant(X))
ResultadosLogit=[Link]()
print([Link]())

Optimization terminated successfully.

Current function value: 0.644516
Iterations 5
Logit Regression Results
==============================================================================
Dependent Variable: Smoker No. Observations: 1196
Model: Logit Df Residuals: 1191
Method: MLE Df Model: 4
Date: xxx, xx xxx xxxx Pseudo R-squ.: 0.02975
Time: xx:xx:xx Log-Likelihood: -770.84
converged True LL-Null: -794.47
LLR p-value: 1.341e-09
===============================================================================
coef standard error z P>|z| 0.025 0.975]
-------------------------------------------------------------------------------
const 2.7451 0.829 3.311 0.001 1.120 4.370
Age -0.0209 0.004 -5.577 0.000 -0.028 -0.014

90
Schooling -0.0910 0.021 -4.402 0.000 -0.131 -0.050
Income 4.72e-06 7.17e-06 0.658 0.510 -9.33e-06 0.0000188
Pcigs79 -0.0223 0.012 -1.789 0.074 -0.047 0.002
===============================================================================

The reader must remember that the direct interpretation of the coefficients in the Logit Model is not valid.
particularly useful at first glance, as these coefficients quantify the impact of a certain variable
about the logarithm of the odds ratio and not directly about the probability of the dependent variable
I took the value of 1.
To understand the impact of a specific variable in a logit model, as mentioned at the beginning
At the beginning of this chapter, we refer to marginal effects. These effects are obtained through the method
.get_margeff(), applied to the model's results instance. The reader should keep in mind that ...
To obtain a 'general' measure of the marginal effect of a specific variable, the existence of two was pointed out.
paths (which are not the only ones): 1) obtaining the marginal effect using the average observation and 2) obtaining the
average of the individual marginal effects.
To obtain the marginal effects using the first 'path', it is enough to assign the value 'mean' to
parameter of the method .get_margeff(); to obtain them by the second 'path', it is enough to assign it the
value of 'overall'. The parameter of the method .get_margeff() corresponds to 'overall' by default, so
that the average marginal effects are those obtained in case the value of at is not specified.
The average marginal effects (second path) are:
In [10]: LogitMargEff = ResultadosLogit.get_margeff()
print([Link]())
Logit Marginal Effects
=====================================
Dep. Variable: Smoker
Method: dydx
At: overall
===============================================================================
dy/dx standard error z P>|z| 0.025 0.975
-------------------------------------------------------------------------------
Age -0.0047 0.001 -5.864 0.000 -0.006 -0.003
Education -0.0206 0.005 -4.539 0.000 -0.030 -0.012
Income 1.069e-06 1.62e-06 0.659 0.510 -2.11e-06 4.25e-06
Pcigs79 -0.0051 0.003 -1.798 0.072 -0.011 0.000
===============================================================================

This table, in addition to the specific value of the marginal effects, also presents the confidence intervals for
corresponding bonds and indicate which method is used to calculate them (At:). Thus, we observe, for example,
that an increase of one year in age is associated, on average, with a decrease of 0.47 points
percentages in the probability of being a smoker and that an increase of 1 year in education is associated, in
on average, a decrease of 2 percentage points in the probability of being a smoker.
For the marginal effects at the mean (first path), it is only necessary to assign 'mean' to the parameter.
The code is as follows:
In [11]: LogitMargEff = ResultadosLogit.get_margeff(at="mean")
print([Link]())

91
 Logit Marginal Effects
=====================================
Dependent Variable: Smoker
Method: dydx
At: mean
===============================================================================
dy/dx standard error z P>|z| 0.025 0.975
-------------------------------------------------------------------------------
Age -0.0049 0.001 -5.598 0.000 -0.007 -0.003
Schooling -0.0213 0.005 -4.411 0.000 -0.031 -0.012
Income 1.107e-06 1.68e-06 0.658 0.510 -2.19e-06 4.4e-06
Pcigs79 -0.0052 0.003 -1.790 0.073 -0.011 0.000
===============================================================================

Thus, we observe, for example, that an increase of one year in age is associated, on average, with a
a decrease of 0.49 percentage points in the probability of being a smoker and that an increase of 1 year in the
education is associated, on average, with a decrease of 2.13 percentage points in the likelihood of
being a smoker.

Probit Model
Similar to the Logit Model, a popular alternative to the MLP is the Probit Model. In this model, the pro-
The probability that the dependent variable takes the value of 1 is given by the following:
xTi
Pi= FxTβ =i Φ x Tβ = I Z φ(z)dz
−∞
So, FxTi β corresponds to the Cumulative Distribution Function of the Standard Normal Distribution. To
just like with the Logit Model, the probabilities obtained are not less than 0 or more than 1, and the method of
the parameter estimation is Maximum Likelihood (ML). To the reader interested in learning more
detail the Probit Model and its estimation process, it is recommended to consult Chapter 15, Models of
qualitative response regression, by Gujarati and Porter (2010).
As in the case of the Logit Model, the direct interpretation of the coefficients associated with the variables
explanatory ones are not very useful and the corresponding marginal effects differ from these. For the
Probit model, the marginal effect of the variable [Link] given by:

∂FxTi β
= φ xTi β β j
x yes
whereφ(.) It is the Probability Density Function of the Standard Normal Distribution. How can it
It can be observed that the marginal effects in the Probit Model, just like in the Logit Model, are not given directly.
They do not only depend on the coefficient associated with the variable, but also on the values of the explanatory variables.
Thus, in a dataset with observations, values are obtained for the marginal effect of the same variable;
To obtain a 'general' measure, the two previously indicated paths can be followed in the case of the Model.
Logit.
In the first approach, the average values of the variables are used in order to find a single
marginal effect for a certain variable. Defining again the vector that contains the average values
God of the explanatory variables comox¯, the marginal effect of the variable xjis given by:

92
 ∂P
= φ xbarITβ
∂x j
Thus, the 'general' marginal effect of a particular variable corresponds to the marginal effect of that variable.
calculated with the values corresponding to the average observation.
In the second path, the marginal effects corresponding to the same variable are calculated.
average is obtained. Specifically:

∂P ∑ in=1φ x Tβi
= βj
∂x j n
Thus, the 'general' marginal effect determined by the variable corresponds to the average of the marginal effects.
individuals for that variable.
Up to this point, a very brief description has been given of the general idea on which the Model is based.
Probit; now, it is time to put the acquired knowledge into practice by carrying out the construction
and estimation of a probit model in Python.

Probit Model in Python

The dataset that will be used for the construction and estimation of the Probit Model is the same as that used for
the MLP and the Logit Model. Since the construction and estimation of probit models is achieved with functions
deStatsmodels do not require instruments from libraries other than those already imported, there is
with all the tools to be able to work properly. To the user who has not followed the exercise.
corresponding to the MLP section or to the Logit Model section, you are invited to run the code cells
#1, #2 and #6, which are essential for the work that is now intended to be developed.
To build a probit model in Python, the Probit(...) function from Statsmodels is used; its arguments...
Mentos include the dependent variable, the explanatory variables, and the missing parameter (with default value)
’none’) which allows indicating the action to take in case there are missing values.
For the estimation, the method is used, just like in the models we have discussed in previous chapters.
.fit(); for the visualization of results, the method .summary() (applied on the instance of results-
two) and the function print(...) (taking as argument the application of the method .summary()). The code
corresponding to our case is:

In [12]: ProbitModel=[Link](Y, sm.add_constant(X))

ResultadosProbit=[Link]()
print([Link]())

Optimization terminated successfully.

Current function value: 0.644304
Iterations 5
Probit Regression Results
==============================================================================
Dependent Variable: Smoker No. Observations: 1196
Model: Probit Df Residuals: 1191
Method: MLE Df Model: 4
Date: xxx, xx xxx xxxx Pseudo R-squ.: 0.03007
Time: xx:xx:xx Log-Likelihood: -770.59
converged True LL-Null: -794.47
LLR p-value: 1.052e-09

93
===============================================================================
coef std err P>|z| z [0.025 0.975]
-------------------------------------------------------------------------------
const 1.7019 0.511 3.333 0.001 0.701 2.703
Age -0.0130 0.002 -5.655 0.000 -0.017 -0.008
Schooling -0.0562 0.013 -4.450 0.000 -0.081 -0.031
Income 2.72e-06 4.4e-06 0.619 0.536 -5.9e-06 0.0000113
Pcigs79 -0.0138 0.008 -1.792 0.073 -0.029 0.001
===============================================================================

The reader should remember that the direct interpretation of the coefficients in the Probit Model is not valid.
particularly useful at first glance and to understand the impact of a specific variable in a model
probit, as mentioned at the beginning of this chapter, marginal effects are used.
Marginal effects are obtained through the method .get_margeff(), applied to the instance of re-
results of the model. The reader must keep in mind that to obtain a 'general' measure of the marginal effect
From a specific variable, the existence of two paths (which are not the only ones) was pointed out: 1) to obtain the effect
marginal using the average observation and 2) obtain the average of the individual marginal effects.
To obtain the marginal effects using the first 'path', it is enough to assign the value 'mean' to the parameter.
Call the method .get_margeff(); to obtain them through the second "path", it is enough to assign it the value.
of 'overall'. The parameter of the method .get_margeff() corresponds to 'overall' by default, so the
Average marginal effects are those obtained in case the value of at is not specified.
The average marginal effects (second path) are:

In [13]:ProbitMargEff=ResultadosProbit.get_margeff()
print([Link]())

Probit Marginal Effects

=====================================
Dependent Variable: Smoker
Method: dydx
At: overall
===============================================================================
dy/dx standard error z P>|z| 0.025 0.975
-------------------------------------------------------------------------------
Age -0.0048 0.001 -5.899 0.000 -0.006 -0.003
Schooling -0.0207 0.005 -4.565 0.000 -0.030 -0.012
Income 1.002e-06 1.62e-06 0.619 0.536 -2.17e-06 4.18e-06
Pcigs79 -0.0051 0.003 -1.800 0.072 -0.011 0.000
===============================================================================

The information presented in this table is the same as that obtained in the case of the Logit model. Thus, ob-
we observe, for example, that an increase of one year in age is associated, on average, with a decrease
by 0.48 percentage points in the probability of being a smoker and that an increase of 1 year in education
is associated, on average, with a decrease of 2 percentage points in the probability of being a smoker.
For the marginal effects at the mean (first path), it is only necessary to modify the value of the parameter.
the code is as follows:

94
In [14]:ProbitMargEff=ResultadosProbit.get_margeff(at="mean")
print([Link]())

Probit Marginal Effects

=====================================
Dependent Variable: Smoker
Method: dydx
At: mean
===============================================================================
dy/dx standard error z P>|z| 0.025 0.975
-------------------------------------------------------------------------------
Age -0.0049 0.001 -5.666 0.000 -0.007 -0.003
Schooling -0.0213 0.005 -4.455 0.000 -0.031 -0.012
Income 1.032e-06 0.00000167 0.619 0.536 -2.24e-06 4.3e-06
Pcigs79 -0.0052 0.003 -1.792 0.073 -0.011 0.000
===============================================================================

Using this method, we observe, for example, that an increase of one year in age is associated with,
on average, a decrease of 0.49 percentage points in the probability of being a smoker and that a
A one-year increase in schooling is associated, on average, with a decrease of 2.13 percentage points.
in the probability of being a smoker.

Model comparison
It is useful to be able to compare the results of the different models studied when having their availability.
information in a single space; that is, having a summary table in which it is presented simultaneously-
Again, the basic information of the results instance corresponding to each model.
It is possible to compare results of models created with Statsmodels functions using the function
summary_col(...) from the [Link].summary2 module. This function allows summarizing multiple
instances of results presenting them side by side; some of its arguments are the list of ins-
instances of results to compare (results), the list of names to assign to the summary of each instance (model_na-
and the parameter stars (of boolean type), which indicates whether or not to visualize the significance through
asterisks (*)
The first step to take in order to compare the instances of results from the models is to import the fun-
action that allows carrying out such a task. The import line of the function summary_col(...) from Statsmo-
[Link].summary2es the following:

In [15]:from [Link].summary2 importsummary_col

Now, having the required function available, it is enough to specify the corresponding arguments.
and run the code to get the summary table (which we will call MyTable):

In [16]: MyTable = summary_col(results=[MLPResults, LogitResults,

Probit Results
stars=True,
model_names=["Modelo MLP","Modelo Logit",
Probit Model
print(MyTable)

95
=================================================
Modelo MLP Modelo Logit Modelo Probit
-------------------------------------------------
const 1.1231*** 2.7451*** 1.7019***
(0.1884) (0.8292) 0.5106
Age -0.0047*** -0.0209*** -0.0130***
(0.0008) (0.0037) 0.0023
Education -0.0206*** -0.0910*** -0.0562***
(0.0046) (0.0207) 0.0126
Income 0.0000 0.0000 0.0000
(0.0000) 0.0000 (0.0000)
Pcigs79 -0.0051* -0.0223* -0.0138*
(0.0029) 0.0125 (0.0077)
=================================================
Standard errors in parentheses.
* p<.1, ** p<.05, ***p<.01

The reader should remember that the direct interpretation of the coefficients of the Logit and Probit models does not
it is particularly useful and it is not possible to compare its magnitude with that of the corresponding to the Model
Linear Probability, as it quantifies different relationships. Due to this, effects have been resorted to.
marginales to obtain a comparable measure, between models, of the effect of an explanatory variable
about the probability that the dependent variable takes the value of 1. The summary table we have just
present the coefficients of each model and the marginal effects, does this mean that it is not possible
compare the information presented?
Although it is not possible to directly compare the coefficients of the MLP, Logit, and Probit models, and for
Thus, we use the marginal effects, the reader should keep in mind that there is an approximate relationship.
among the magnitudes of those coefficients.
Following Katchova (2013), the coefficients of the MLP, Logit, and Probit models are characterized by the
following relationships:

β Logit' 4β MLP

β Probit' 2.5β MLP

β Logit' 1.6β Probit

Thus, it is enough to apply the corresponding transformations to the coefficients of the Logit and Probit models.
to obtain directly comparable measures with the coefficients of the Linear Probability Model; without
embargo, duties to rise from said transformations are based on approximate relationships. Therefore, if you want
to compare the models, it is most appropriate to base such comparison on the marginal effects.
Statsmodels does not offer a function or method (as of the publication of this work and to the knowledge of
author) that allows comparing the marginal effects of several models; however, this situation should not dis-
encourage the user, as it is still possible to obtain a summary table that presents the marginal effects.
Now, something very important to keep in mind is that the work associated with achieving this objective
it can be somewhat tedious (or even complicated for the novice) and must be developed carefully.

96
 The goal we aim for is to obtain a summary table that presents the marginal effects.
corresponding to the Linear Probability Model, to the Logit Model, and to the Probit Model, so the data
the essentials we need are the values of such marginal effects. The first thing we will do is create
DataFrames to store the corresponding data.
The information on the marginal effects of the Logit Model and the Probit Model has been stored with
specific objects' anteriority; now, we will extract only the magnitudes of these marginal effects
(with the [Link], applied to the object that contains the information on the marginal effects regarding
you) and we will use them as values for the parameter data in the function DataFrame(...) of pandas. Likewise,
For the parameter index, we will use a Python list "[...]" that contains the names of the co-variables.
respondents and we will assign to the parameter column the name of the respective model. The procedure is the
the same for the Logit and Probit models, as evidenced below:

In [17]:Coeficientes=["Edad","Escolaridad","Ingreso","Pcigs79"]
EfectosLogit=[Link]([Link], index=Coeficientes,
columns=["Logit"])
EffectProbit=[Link]([Link], index=Coefficients,
columns=["Probit"])

We already have the data on the marginal effects of the Logit and Probit models, but what about the
marginal effects of the MLP? The reader should remember that the marginal effects of the Linear Probability Model
likelihood corresponds directly to the coefficients, so unlike the Logit and Probit models, their
The result of applying this attribute differs.
that is obtained with the application of the method .get_margeff() of the Logit and Probit models, so that
the procedure to extract the magnitudes of the marginal effects (coefficients) will be somewhat different from the
that we have developed previously.
The first thing the user should notice is that the [Link] generates a result that includes
the coefficient corresponding to the intercept term, however, for the Logit and Probit models, as re-
It is logical, no marginal effect corresponding to this is reported. Therefore, in order to compare
the marginal effects of the same variables, the intercept term will be ignoredβ 0 ).
To carry out the 'extraction' of the magnitudes of all the coefficients except for the term of the inter-
the brackets "[ ]" are used, applied to the object that contains the parameters of the MLP (previously
we store it as a Pandas series). Within the brackets, we will use the sequence 1:5; this,
with the purpose of ignoring the term of the intercept. Additionally, the name of the model will be assigned to
parameter columns.
Note: To the reader who is unaware, it is informed that a Pandas Series indexes (by default) its ob-
observations with integers starting from 0 and that the selection by brackets "[]" with integers is inclusive
in the first value and exclusive in the last, so the use of a sequence 1:5 within cor-
[ ] selection means 'ignore the first observation (identified by 0) and include
all the others up to the one indexed with the number 4 (does not include the one with index 5).

The code to execute for the MLP case is as follows:

In [18]: EfectosMLP=[Link](MLPMargEff[1:5], columns=["MLP"])

The DataFrames containing the magnitudes of the marginal effects of the models are already available.
Logit, Probit, and MLP. Now, it is time to group these DataFrames into one; to achieve this, it is used
the concat(...) function of pandas; in its arguments, a Python list '[...]' is included that contains the
DataFrame to concatenate and the axis parameter with a value of 1, so that the concatenation is achieved by placing
a DataFrame next to another:

97
In [19]:[Link]([EfectosMLP, EfectosLogit, EfectosProbit], axis= 1)

Out[19]: MLP Logit Probit

Age -0.004726 -0.004890 -0.004920
Schooling -0.020613 -0.021334 -0.021340
Income 0.000001 0.000001 0.000001
Pcigs79 ["-0.005132","-0.005234","-0.005235"]

Now there is a single summary table that presents the marginal effects of the different
two models: magnitudes that are indeed directly comparable. We can observe that the estimated effects
do not differ significantly: for example, the marginal effect of the variable Schooling of -2.06 p.p in the
MLP, -2.13 p.p in the Logit model and -2.13 p.p in the Probit model; for the Age variable the marginal effect
ranging from -0.4726 p.p (MLP) to -0.492 p.p (Probit).
At this point, this chapter closes, in which we have dealt with models with a binary dependent variable;
Now, we will extend the scope of discrete dependent variable models by including regressors.
which take only two values, that is to say, binary variables.

98
Chapter 7

Discrete dependent variable - Non-binary case

In the previous chapter, we worked with models where the dependent variable was discrete and, furthermore, only
it could take two values; in this chapter, the topic of the discrete dependent variable will be continued, but
the restriction that it can only take two values will be removed, that is, we will address the topic of variables
discrete non-binary.
The dependent variables we will work with will be categorical, but they will not be limited to only
two categories, so they will not take only two values (0 and 1) but as many values as alternatives
have that variable.
The brief exposition of the theoretical foundations of the model to be developed in this chapter is based on
in its entirety, in Katchova (2013, 2), therefore, if necessary, the reader is invited to consult it.
source for greater clarity on this matter.
The dependent variables to be used are non-binary discrete, so they take values from a set.
finished composed of more than two options. It is possible that some examples may have certain value to illustrate the
nature of the mentioned variables:

When a person wants to consume a beer, they have a finite and well-defined set at their disposal.
nest of beer brands that you can acquire. These brands can be coded using numbers,
assigning a different number to each one, so the resulting variable corresponds to a dis- variable
non-binary (assuming there are more than two brands of beer available), as there are
multiple, but limited, possibilities.

When someone attends a movie theater, they have different food 'combos' available to them. There is a
well-defined set of 'combos' and their quantity is limited; in this way, a number can be assigned
different from each combo, with the purpose of identifying it. Thus, the resulting variable from the coding of
The 'combos' is a non-binary discrete variable (assuming there are more than two combos).

A discrete variable with categories (or alternatives) undergoes a coding process using
numbers; these numbers have no meaning by themselves, as their magnitude is not interpretable and
its function is exclusively of coders. Likewise, it should be noted that the order in which they are assigned
numbers are free, as they lack their own meaning and do not have ordinal nature, so they do not
they establish a defined order for the alternatives of the variables.
The variable has multiple categories, but it is important to note that a specific individual (entity)
only a single alternative corresponds; that is, each individual or entity, of the multiple alternatives to dis-
position, you can choose only one. Returning to the previous examples, for the case of beers, the
The individual can only choose one brand of beer that they will consume at this moment and, in the case of the
food combos, the person should only choose a single combo.

99
 For this type of models, the information is presented in two specific formats: wide and long. In the format
wide, the information for each individual (entity) is presented in a single row, therefore= j( j=
1, 2, 3, ..., m )In the long format, the information for each individual (entity) is presented in rows, each one
of which corresponds to one of the alternatives of the dependent variable, which takes the value of 1 one
only once and from 0m−
Once for each individual (entity).
Using a somewhat more illustrative description, the same information is presented below,
in wide format in long format. Reverting to the example of beers, considering only 5
brands (coded as follows: Aguila: 1, Poker: 2, [Link]: 3, Corona: 4, Heineken:
5), the following information is presented in wide format:

This same information, in long format, corresponds to the following table:

Thus, the reader can appreciate that for the wide format= j, while for the formatlong:

100
 1, I am= j
y j= (
0, I= j
Something important to mention is that in models of non-binary discrete dependent variable,
Just like in models of discrete binary dependent variables, what matters is not the value of the variable.
dependent as such but rather the probability that it takes a specific value. Thus, in the Logit Model,
since the dependent variable was binary, there were only two possibilities, so the probability was estimated
the probability that the dependent variable will take the value of 1, and the probability that it will take the value of 0
it simply corresponded to one less the probability of it taking the value of 1. However, in the model
Multinomial Logit (which is what we will work with from now on) does not only have two alternatives, so,
What probability is estimated?
The reader should know that in the Multinomial Logit model what is sought is to estimate the probability that
choose the alternative within the alternatives that the dependent variable has. For variable models
discrete non-binary dependent, the probability that the individual chooses the alternative is given by:

Pyes= Prob(yi= j) = Fj( xi γ)

Depending on the functional form that F adoptsj a different model will be obtained: Multinomial Logit,
Probit Multinomial, Conditional Logit, Ordered Logit, etc. This work will focus solely on the Model
Multinomial Logit; this is due to the availability of effective and easy-to-use tools in the libraries.
Specialized Python.
Before proceeding to the practical work in Python, it is important for the reader to take into account the types
of explanatory variables used in non-binary discrete dependent variable models, in terms of
its classification is relevant; the explanatory variables can be alternative-invariant regressors (also
case-specific regressors or alternative-variant regressors (also called alternative-specific regressors).
The explanatory variables of the type alternative-invariant regressors have no direct correspondence with
alternative finished; the values of these variables vary among individuals (entities) but do not vary among
alternatives for the same individual (entity). Picking up the example of the beer that we have discussed with an-
terioridad, la variableIngresoes una variablealternative-invarianten cuanto tiene el mismo valor para el mismo
individual and does not change depending on the chosen alternative; for Individual 1, the income remains COP
2,300,000 regardless of which brand of beer is chosen; for Individual 2, the income remains COP
3,000,000 no matter which beer brand you choose, etc.
The explanatory variables of the alternative-variant type regressors have a direct correspondence with determi-
no alternative, so its value depends on it; the values of these variables vary between alternatives and
they can also vary between individuals. Once again, referring to the example of beer consumption, the pre-
cio es una variablealternative-variant, en cuanto tiene distintos valores para distintas alternativas; asimismo, se
shows some variation among individuals. Being somewhat more explicit, the reader can observe that the pre-
the 'Eagle' die is not the same as the 'Crown' or the 'Poker'; this is referred to with alternative-variant,
Well, although it is about the same characteristic (price), its value depends on the corresponding alternative.
An additional aspect to consider is that the explanatory variables type alternative-invariant regressors
they speak in the Multinomial Logit Model and the explanatory variables type alternative-variant regressors are used in
the Conditional Logit Model and Mixed Logit. In the practical work that we will carry out with Python, as has been
As previously indicated, we will only address the Multinomial Logit Model, so only
explanatory variables type alternative-invariant during the construction and estimation process.
Returning to the description of the Multinomial Logit Model, the adopted functional form is as follows:
T
ewI γj
Pyes= Prob(yI= j) =
∑m wi γh T
h =1 e

101
 whereyesit is the probability that the individual (entity) chooses the alternative j, wi it is the vector of
alternative-invariant variables for the individual (entity) yγ jit is the set of corresponding coefficients
to the alternatives. To achieve the probability estimation, one of the sets of coefficients is normalized.
to 0, so that the interpretation of others is done with reference to the base category (to which it corresponds
the set of coefficients 0).
One point to highlight is that the magnitude of the coefficients of the Multinomial Logit Model, just like that of
the Logit and Probit should not be interpreted directly; these coefficients only indicate whether it increases (when
the coefficient is positive) or decreases (when the coefficient is negative) the probability that is chosen the
alternatives, in comparison to the base category, but do not quantify the variation.
To have a measure that quantifies the impact of a certain variable on the probability of
to choose the alternative, resort to the marginal effects, which, for the Multinomial Logit Model, are
given by the following:
T
ewi γj
Phi= T T T
ewi γalt1 + e wi γalt2 + ... + e wi γaltm
where it refers to the alternative j (alt1 refers to alternative 1, alt2 refers to alternative
2 and so on). Now, defining wTi γalthhowhi, it is given:

ezyes
Pyes=
ez1i+ ez2i+ ... + ezme
Therefore, ashi= wTi γalth= γh0+ γh1 w1i+ ... + γhr writhe first subscript (in the coefficients
It refers to the alternative to which the coefficients correspond (malternatives), the second subin-
the coefficients (first subscript for the variables) refers to the alternative-invariant variable
(alternative-invariant variables) and the second subscript of the variables refers to the individual (entity),
then it must be:

∂Pyes γ jefz yes( h=∑1emzhi)−ez yes ∑m

h=1γh fe
z hi
= 2
∂w f (∑m z
h=1e hi )

∑)mez zhi
∂Pyes γ jef( zhji=1and yes
∑m
h=1γh fe
z hi
= 2 − 2
∂w f (∑m z
h=1e hi ) (∑m z
h=1e hi )

∂Pyes ezyes ∑m
h=1γh fe
z hi
= m z γj f −
∂w f ∑ h=1e hi ∑m z
h=1e hi

∂Phello m

∂w f = Pyesγ j f − ∑ Phi γh f !
h =1

It should be noted that marginal effects do not necessarily have the same sign as the coefficients; for
Therefore, the coefficients indicate (do not measure) whether the probability of choosing the alternative increases or decreases.
comparison with the base category, and the marginal effects quantify the variation in the probability of choosing
a specific alternative in response to the variation in an explanatory variable (alternative-invariant).
As with the Logit and Probit models studied in the previous chapter, in the Multinomial Logit Model-
mial, for the same variable different marginal effects are obtained (there are individuals (entities)). For the
Thus, to obtain a general measure of the marginal effect of a particular variable, one can use the
two "paths" presented in Chapter 6 of this work (to the reader who is unclear about this, it is suggested
consult that chapter, as the logic applied to the current one is the same).
It is important to note that the marginal effects of the same variable must sum to zero among the alternates.
tivas, that is to say:

102
 m ∂Pj
∑∂w f =0
j =1

∂Pj
where it is the marginal effect of the variablef about the probability of choosing the alternative and there is
∂w f
ten alternatives. It is also important to bring up that regardless of the alternative chosen as
base category, although the coefficients may change, the marginal effects will be the same.
Up to this point, a very brief exposition of the theoretical foundations on which it is based has been given.
structure the model that we now intend to address practically in Python. It is possible that the structure
the mathematics on which the Multinomial Logit Model is built is not completely understandable for the
reader; however, the practical work may help you gain greater clarity on certain matters.
aspects.

Multinomial Logit Model in Python

The dataset that will be used for the practical work in Python is the one used in Katchova (2013, 3) (the data-
The sets used by Professor Katchova come, in their original version, from Herriges and King (1999) and Cameron.
y Trivedi (2005)). The dataset has been obtained from the Econometrics Academy website [Link]
What is intended to study is how income, which is a variable, affects
of alternative-invariant type, to the choice of the fishing type ('mode') of the individuals.
The functions for working with Multinomial Logit models belong to the Statsmodels library, which
not integrated into the Anaconda distribution; likewise, the other libraries that provide support for analysis
econometric that will be advanced (pandasyNumPy) are also included in this distribution, therefore
No additional library installation is required from the user.

Preparation of the environment

The first cell in which the user will write and execute code will have the following content:

In [1]: import numpy as np

import pandas as pd
[Link]

The previous block of code allows loading the necessary tools to carry out the actions.
required in the proposed econometric analysis; in case the reader is not clear about this, they are
please refer to the section 'Preparing the Environment' of Chapter 2, where an explanation is provided about
such matter.

Data importation
The import of data is carried out using functions from the pandas library. Depending on the type of
file in which the information of the dataset is stored, we will need to resort to a partition function-
cellular and a set of specific parameters.
For our specific case, the following code is the one that imports the dataset:

In [2]: data = pd.read_csv("multinomial_fishing1.csv")

103
 Once the line of code in this cell has been executed, the dataset should have been imported correctly.
To check it, we resort to the .head() method, applied to the object that contains the dataset (in this case he
we have assigned the name data):

In [3]:[Link]()

Out[3]: mode price catchrate [Link] [Link] [Link] [Link] \

0 charter 182.930 0.5391 0 0 0 1
1 charter 34.534 0.4671 0 0 0 1
2 private 24.334 0.2413 0 0 1 0
3 pier 15.134 0.0789 0 1 0 0
4 private 41.514 0.1082 0 0 1 0

[Link] [Link] [Link] [Link] [Link] \

0 157.930 157.930 157.930 182.930 0.0678
1 15.114 15.114 10.534 34.534 0.1049
2 161.874 161.874 24.334 59.334 0.5333
3 15.134 15.134 55.930 84.930 0.0678
4 106.930 106.930 41.514 71.014 0.0678

[Link] [Link] [Link] income mode1

0 0.0503 0.2601 0.5391 7.083332 4
1 0.0451 0.1574 0.4671 1.250000 4
2 0.4522 0.2413 1.0266 3.750000 3
3 0.0789 0.1643 0.5391 2.083333 2
4 0.0503 0.1082 0.3240 4.583332 3

We can observe that the import process has been successful and the information from the dataset has been soul-
it was correctly done in a pandas DataFrame. Now we can calmly continue with the development
de nuestro análisis, en cuanto hemos concluido satisfactoriamente el primer paso. Procederemos a construir y
to estimate a multinomial logit model, which takes as the dependent variable a non-binary discrete variable;
Therefore, it is useful to know the alternatives that this variable contains.
To know the alternatives of the dependent variable, it is enough to select it (using the brackets
'[ ]') and apply the .unique() method, which allows you to know the unique values that an object has. The code
to execute is:

In [4]: data["mode"].unique()

["charter","private","pier","beach"]

The result we obtain is a NumPy ndarray, in which we are informed of the unique values.
(alternatives) that the dependent variable takes are 'charter', 'private', 'pier' and 'beach'; therefore, we count
with four alternatives and, consequently, we will obtain three sets of coefficients (the reader should already know
the why.
For our case, it is very simple to know the number of alternatives available, as it is enough to
list them; however, if the variable had many alternatives, it would be tedious to count them one by one.
pandas offers an effective solution to this issue: instead of using the .unique() method, the
mé[Link](), which counts the number of unique values. To verify that indeed there are.
4 alternatives, the following code is executed:

104
In [5]: data['mode'].nunique()

Out[5]:4

Data preparation
The reader who has followed this work up to this point will have noticed that in the previous chapters we have
stored various variables in specific objects. We have done this to group these variables of
according to its particular role in the regression and trying to simplify the code to be used afterwards. Now,
In this exercise, although only two variables will be used, the variables will be assigned to specific objects.
to simplify the syntax a bit and facilitate the understanding of the code. We will proceed to make the assignment
extracting the variables using the brackets "[ ]", just as we have done in previous chapters:

In [6]: Y=data["mode"]
X=data["income"]

Construction and estimation

For the construction and estimation of the model, the function MNLogit(...) from Statsmodels is used.
[Link](); for visualizing the results, the function print(...) and the method summary(). The
the code to execute is the following:

In [7]:ModelLogitMN=[Link](Y, sm.add_constant(X))
ResultadosLogitMN=[Link]()
print([Link]())

Optimization terminated successfully.

Current function value: 1.249704
Iterations 5
MNLogit Regression Results
==============================================================================
Dependent Variable: mode No. Observations: 1182
Model: MNLogit Df Residuals: 1176
Method: MLE Df Model: 3
Date: xxx, xx xxx xxxx Pseudo R-squared: 0.01374
Time: xx:xx:xx Log-Likelihood: -1477.2
converged True LL-Null: -1497.7
LLR p-value: 6.093e-09
================================================================================
mode=charter coef standard error z P>|z| 0.025 0.975
--------------------------------------------------------------------------------
const 1.3413 0.195 6.896 0.000 0.960 1.723
income -0.0316 0.042 -0.756 0.450 -0.114 0.050
--------------------------------------------------------------------------------
mode=pier coef standard error z P>|z| 0.025 0.975
------------------------------------------------------------------------------
const 0.8142 0.229 3.561 0.000 0.366 1.262
income -0.1434 0.053 -2.691 0.007 -0.248 -0.039
------------------------------------------------------------------------------

105
mode=private coef standard error z P>|z| 0.025 0.975
--------------------------------------------------------------------------------
const 0.7389 0.197 3.756 0.000 0.353 1.125
income 0.0919 0.041 2.260 0.024 0.012 0.172
================================================================================

Note: The information with which the function MNLogit(...) of Statsmodels works must be contained
in a wide format dataset. In case of having a long format, the dataset must be reorganized from
so that, before using your information in the function MNLogit(...) of Statsmodels, its format is
wide.

The reader should remember that it was mentioned that, for the estimation, the set of coefficients of an alternative
it is normalized to zero and the interpretation of the obtained coefficients is made with reference to the category
base, however, what is the base category?
For this specific example, the base category is the alternative 'beach', but why? The answer is
it is in the logic that applies the function MNLogit(...) from Statsmodels: when passed as a dependent variable
given a variable not encoded by numbers (the reader should note that the variable mode is not encoded),
the function automatically executes the encoding alphabetically and takes the alternate as the base category
that takes the first position in that order. Thus, as the alternatives of the variable mode are 'charter',
'private', 'beach' and 'pier', the base category is 'beach', as it takes the first position if organized
alphabetically such alternatives.
Now, the interpretation of the coefficients must be made with respect to the alternative 'beach'. Thus, we have
queen comparison with fishing on the beach, a higher income is associated with a lower probability
fishing on a rented boat ('charter') or at the pier ('pier') and a higher probability of fishing on a private boat
private
We have said that a higher income is associated, in comparison to fishing on the beach, with a probability
more or less of certain alternatives, but how much is this greater or lesser probability? Remember
that the coefficients do not quantify the relationship, they only indicate its direction; to obtain a measure of the im-
the pact of a certain variable on the probability of choosing a specific alternative employs the effects
marginals.
To obtain the marginal effects of a multinomial logit model, it is sufficient to apply the method
.get_margeff() to the results instance of that model; to visualize these marginal effects, so
only the use of the method .summary() and the function print(...) is required, as follows:

In [8]: LogitMNMargEff = ResultadosLogitMN.get_margeff()

print([Link]())

MNLogit Marginal Effects

=====================================
Dependent Variable: mode
Method: dydx
At: overall
================================================================================
mode=beach dy/dx standard error z P>|z| 0.025 0.975
--------------------------------------------------------------------------------
income 0.0002 0.004 0.044 0.965 -0.007 0.008
--------------------------------------------------------------------------------
mode=charter dy/dx standard error P>|z| [0.025z 0.975]

106
--------------------------------------------------------------------------------
income -0.0112 0.006 -1.876 0.061 -0.023 0.000
--------------------------------------------------------------------------------
mode=pier dy/dx standard error P>|z| [0.025z 0.975]
------------------------------------------------------------------------------
income -0.0208 0.005 -4.040 0.000 -0.031 -0.011
------------------------------------------------------------------------------
mode=private dy/dx standard error P>|z| [0.025z 0.975]
--------------------------------------------------------------------------------
income 0.0318 0.005 6.039 0.000 0.021 0.042
================================================================================

In this way, an increase of one unit in income is associated with fishing in the
beach 0.02 percentage points more likely; than fishing on a rented boat (charter) being 1.12
2.08 percentage points less likely; that fishing at the pier is 2.08 percentage points less.
likely; and that fishing on a private boat is 3.18 percentage points more likely. (The reader
you can verify that the sum of these marginal effects is zero
The reader should keep in mind that there are several ways to obtain a general measure of the effect.
the marginal effect of a certain variable; the marginal effects that are calculated by default by the method
.get_margeff() are the average marginal effects. If what you want to know are the marginal effects
in the average, the value of 'mean' should be assigned to the parameter of the method .get_margeff(), as follows:

In [9]: LogitMNMargEff = ResultadosLogitMN.get_margeff(at='mean')

print([Link]())

MNLogit Marginal Effects

=====================================
Dependent Variable: mode
Method: dydx
At: mean
================================================================================
mode=beach dy/dx standard error z P>|z| 0.025 0.975
--------------------------------------------------------------------------------
income 7.496e-05 0.004 0.019 0.985 -0.008 0.008
--------------------------------------------------------------------------------
mode=charter dy/dx standard error P>|z| [0.025z 0.975]
--------------------------------------------------------------------------------
income -0.0120 0.006 -1.977 0.048 -0.024 -0.000
--------------------------------------------------------------------------------
mode=pier dy/dx standard error P>|z| [0.025z 0.975]
------------------------------------------------------------------------------
income -0.0207 0.005 -4.239 0.000 -0.030 -0.011
------------------------------------------------------------------------------
mode=private dy/dx standard error P>|z| [0.025z 0.975]
--------------------------------------------------------------------------------
income 0.0326 0.006 5.727 0.000 0.021 0.044
================================================================================

107
 Now, it is known that an increase of one unit in income is associated with fishing on the beach.
('beach') sea 0.008 p.p more likely; that fishing on a rented boat ('charter') is 1.2 p.p less likely;
that fishing from the pier is 2.07 p.p less likely; and that fishing from a private boat is
3.26 p.p more likely. (Again, the reader can verify that the sum of these marginal effects is
zero).
If the reader compares the obtained results with those of Katchova (2013, 3), they will notice that the ones she obtains
they are different; this is because in the exercise developed by Professor Katchova it is not taken as
base category 'beach' if not 'charter', is it possible for us to do the same? The answer is yes.
The function MNLogit(...) from Statsmodels takes as the base category the alternative that occupies the first place.
position in alphabetical or numerical order and it is not possible to modify such logic in its operation. However,
if the alternatives are reclassified so that the one desired as the base category occupies the first position
position, then the function MNLogit(...) of Statsmodels will assume that alternative as the category
base.
For example, by using pandas tools, a categorical variable can be decoded using
the method .map(...) and using a Python dictionary "{...}" as an argument; within that dictionary
the keys will correspond to the alternatives and the values to the codes assigned to each one.
For example, below a variable 'mode2' is created that corresponds to the encoded variable 'mode'.
so that the 'charter' alternative is in the first position in numerical order and, therefore, is
taken as a base category by the function MNLogit(...) of Statsmodels:

In [10]:data["mode2"]=data["mode"].map({"charter":"1_charter",
"beach":"2_beach",
"pier":"3_pier",
"private":"4_private"})

Using the same code structure that we have used previously, the following is obtained
model:

In [11]: Y=data["mode2"]
LogitModelMN=[Link](Y, sm.add_constant(X))
LogitMNResults = [Link]()
print([Link]())

Optimization terminated successfully.

Current function value: 1.249704
Iterations 5
MNLogit Regression Results
==============================================================================
Dep. Variable: mode2 No. Observations: 1182
Model: MNLogit Df Residuals: 1176
Method: MLE Df Model: 3
Date: xxx, xx xxx xxxx Pseudo R-squared: 0.01374
Time: xx:xx:xx Log-Likelihood: -1477.2
converged True LL-Null: -1497.7
LLR p-value: 6.093e-09
===================================================================================
mode2=2_beach coef standard error z P>|z| 0.025 0.975
-----------------------------------------------------------------------------------
const -1.3413 0.195 -6.896 0.000 -1.723 -0.960

108
income 0.0316 0.042 0.756 0.450 -0.050 0.114
-----------------------------------------------------------------------------------
mode2=3_pier coef standard error P>|z| [0.025z 0.975]
--------------------------------------------------------------------------------
const -0.5271 0.178 -2.965 0.003 -0.876 -0.179
income -0.1118 0.044 -2.541 0.011 -0.198 -0.026
--------------------------------------------------------------------------------
mode2=4_private coef standard error z P>|z| 0.025 0.975
-----------------------------------------------------------------------------------
const -0.6024 0.136 -4.426 0.000 -0.869 -0.336
income 0.1235 0.028 4.426 0.000 0.069 0.178
===================================================================================

As can be appreciated, the magnitudes of the coefficients have changed and, moreover, their interpretation now
it should not be done with reference to 'beach' but to 'charter'. Thus, it has to be, compared to boat fishing
chartered, a higher income is associated with a greater likelihood of fishing on the beach
and in a private boat ('private') and a lower probability of fishing at the pier ('pier'). What about the effects
marginal?
The reader should remember that, regardless of what the base category is, the marginal effect of a variable
specificity, unlike the coefficient, will always be the same. To verify this, it is enough to apply the method
.get_margeff() to the model's results instance:

In [12]: LogitMNMargEff = ResultadosLogitMN.get_margeff()

print([Link]())

MNLogit Marginal Effects

=====================================
Dependent Variable: mode2
Method: dydx
At: overall
===================================================================================
mode2=1_charter dy/dx standard error z P>|z| [0.025 0.975
-----------------------------------------------------------------------------------
income -0.0112 0.006 -1.876 0.061 -0.023 0.000
-----------------------------------------------------------------------------------
mode2=2_beach dy/dx standard error P>|z| [0.025z 0.975]
---------------------------------------------------------------------------------
income 0.0002 0.004 0.044 0.965 -0.007 0.008
---------------------------------------------------------------------------------
mode2=3_pier dy/dx standard error P>|z| [0.025z 0.975]
--------------------------------------------------------------------------------
income -0.0208 0.005 -4.040 0.000 -0.031 -0.011
--------------------------------------------------------------------------------
mode2=4_private dy/dx standard error P>|z| [0.025z 0.975]
-----------------------------------------------------------------------------------
income 0.0318 0.005 6.039 0.000 0.021 0.042
===================================================================================

109
 It can be observed that the marginal effects obtained from the model with the base category 'charter' are exact-
exactly the same as those obtained from the model with the base category 'beach' (you can examine them one by one if you wish
It wishes). Remembering that the method .get_margeff() calculates average marginal effects by default, if it
what is desired to know are the marginal effects on the mean, it is only necessary to assign the value of 'mean' to
parameterat

In [13]: LogitMNMargEff = ResultadosLogitMN.get_margeff(at='mean')

print([Link]())

MNLogit Marginal Effects

=====================================
Dependent Variable: mode2
Method: dydx
At: mean
===================================================================================
mode2=1_charter dy/dx standard error z P>|z| 0.025 0.975
-----------------------------------------------------------------------------------
income -0.0120 0.006 -1.977 0.048 -0.024 -0.000
-----------------------------------------------------------------------------------
mode2=2_beach dy/dx standard error P>|z| [0.025z 0.975]
---------------------------------------------------------------------------------
income 7.496e-05 0.004 0.019 0.985 -0.008 0.008
---------------------------------------------------------------------------------
mode2=3_pier dy/dx standard error P>|z| [0.025z 0.975]
--------------------------------------------------------------------------------
income -0.0207 0.005 -4.239 0.000 -0.030 -0.011
--------------------------------------------------------------------------------
mode2=4_private dy/dx standard error P>|z| [0.025z 0.975]
-----------------------------------------------------------------------------------
income 0.0326 0.006 5.727 0.000 0.021 0.044
===================================================================================

Again, the marginal effects at the mean are exactly the same for the model with the base category.
'beach' and for the model with the base category 'charter', equality that holds for any other alternative.
of the dependent variable (using, obviously, the same data).
At this point, this chapter concludes, with which this brief work comes to an end, which has sought to offer
illustrative examples of basic econometric analysis using the Python programming language
that some of the most relevant tools for this exercise have been presented, trying to provide to the
brief but useful descriptions of the operation of each of these.

110
References

Chatterjee, S., y Simonoff, J. (2013).Handbook of Regression Analysis. Estados Unidos: John Wiley & Sons,
Inc.

Dormann, C.F., Elith, J., Bacher, S., Buchmann, C., Gudrun, C., Carré, G., García, J.R., Gruber, B., Lafourcade,
B., Leitão, P.J., Münkemüller, T., McClean, C., Osborne, P.E., Reineking, B., Schröder, B., Skidmore, A.K.
Zurell, D., and Lautenbach, S. (2013). Collinearity: a review of methods to deal with it and a simulation
study evaluating their performance. Ecography, 36, 27-46. doi: 10.1111/j.1600-0587.2012.07348.x

Fahrmeier, L., Kneib, T., y Lang, Stefan. (2007).Regression. Modelle, Methoden und Anwendungen. Springer.

Greene, W.H. (2003).Econometric Analysis. New Jersey, Estados Unidos: Pearson Education, Inc.

Gujarati, D.N. and Porter, D.C. (2010). Econometrics. Mexico: McGraw-Hill/Interamericana Editors, S.A.
C.V.

James, G., Witten, D., Hastie, T., and Tibshirani, R. (2013). An Introduction to Statistical Learning with Applications
tions in R. doi: 10.1007/978-1-4614-7138-7

Katchova, A. (2013).Econometrics – Probit and Logit Models. Recuperado de: [Link]

le/d/0BwogTI8d6EEidjg2OGl4b3dxVmc/edit

Katchova, A. (2013,2). Econometrics – Multinomial Probit and Logit Models. Retrieved from:
[Link]

Katchova, A. (2013, 3). Econometrics – Multinomial Probit and Logit Models, Conditional Lo-
git Model, Mixed Logit Model Examples. Retrieved from: [Link]
TI8d6EEiLVZ0N1h3LTBLYlk/edit

Kleiber, C., and Zeileis, A. (2008). Applied Econometrics with R. United States: Springer Science+Business
Media, LLC.

Murray, L., Nguyen, H., Lee, Y-F., Remmenga, M., and Smith, D.W. (2012). Variance Inflation Factors in Re-
Regression Models with dummy variables. Conference on Applied Statistics in Agriculture. doi: 10.4148/2475-
7772.1034

Pandas. (n.d). Python Data Analysis Library. Retrieved from: [Link]

Perktold, J., Seabold, Skipper., and Taylor, J. (2018). Welcome to Statsmodels’s Documentation. Retrieved
from: [Link]

Project Jupyter. (2019). The Jupyter Notebook. Retrieved from: [Link]

111
Python Software Foundation. (n.d.). What is Python? Executive Summary. Retrieved from: [Link]
[Link]/doc/essays/blurb/

Razali, N., and Wah, Y. (2011). Power comparisons of Shapiro-Wilk, Kolmogorov-Smirnov, Lilliefors
and Anderson-Darling tests. Journal of Statistical Modeling and Analytics, 2(1), 21-33. Retrieved from:
[Link]

The Matplotlib development team. (2018). Matplotlib. Retrieved from: [Link]

The SciPy community. (2019). What is NumPy? Retrieved from: [Link]

docs/user/[Link]

Vaart, A.W. van der. (1998).Asymptotic Statistics. Reino Unido: Cambridge University Press.

Verbeek, M. (2004).A guide to modern econometrics. Inglaterra: John Wiley & Sons Ltd.

Waskom, M. (2018). An introduction to seaborn. Retrieved from: [Link]

[Link]

Wooldridge, J.M. (2010).Introducción a la econometría. Un enfoque moderno. México: Cengage Learning

Editors, S.A. de C.V.

Würtz, D., and Katzgraber, H.G. (2009). Precise finite-sample quantiles of the Jarque-Bera adjusted Lagran-
ge multiplier test. ETH Econophysics Working and White Papers Series. Retrieved from: [Link]
[Link]/sites/default/files/[Link]

112