14.
170: Programming for
Economists
1/12/2009-1/16/2009
Melissa Dell
Matt Notowidigdo
Paul Schrimpf
What is 14.170?
• 6.170 is a very popular undergraduate course at MIT titled “Introduction
to Software Engineering.” Goals of the course are threefold:
1. Develop good programming habits
2. Learn how to implement basic algorithms
3. Learn various specific features and details of a popular programming language
(currently Java, but has been Python, Scheme, C, C++ in the past)
• We created a one-week course 14.170 with similar goals:
1. Develop good programming habits
2. Learn how to implement basic algorithms
3. Learn various specific features and details of several popular programming
language (Stata, Perl, Matlab, C)
• Course information
– COURSE WEBPAGE: [Link]
– E-mails ([at] mit dot edu):
– mdell
– noto
– paul_s
COURSE OVERVIEW
Today (MATT): Basic Stata, Intermediate Stata,
MLE and NLLS in Stata
Tuesday (MATT): Mata, GMM, Large data sets,
numerical precision issues
Wednesday (PAUL): Basic, Intermediate and
Advanced Matlab
Thursday (MATT, MELISSA): Perl, GIS
Friday (PAUL): Intro to C, More C (C + Matlab, C +
Stata)
(see syllabus for more details)
Lecture 1, Basic Stata
Basic Stata overview slide
• Basic data management
– Reading, writing data sets
– Generating, re-coding, parsing variables (+ regular expressions, if time is
permitting)
– Built-in functions
– Sorting, merging, reshaping, collapsing
• Programming language details (control structures, variables,
procedures)
– forvalues, foreach, while, if, in
– Global, local, and temporary variables
– Missing variables (worst programming language design decision in all of
Stata)
• Programming “best practices”
– Comments!
– Assertions
– Summaries, tabulations (and LOOK at them!)
• Commonly-used built-in features
– Regression and post-estimation commands
– Outputting results
Data Management
• The key manual is “Stata Data Management”
• You should know almost every command in the book
very well before you prepare a data set for a project
• Avoid re-inventing the wheel
• We will go over the most commonly needed commands
(but we will not go over all of them)
• Type “help command” to find out more in Stata, e.g.
“help infile”
• Standard RA “prepare data set” project
1. Read in data
2. Effectively summarize/tabulate data, present graphs
3. Prepare data set for analysis (generate, reshape, parse,
encode, recode)
4. Preliminary regressions and output results
Getting started
• There are several ways to use Stata …
DO file editor
Stata user Text editor (e.g.
operating system interface emacs, TextPad)
Windows (A) (B)
UNIX / Linux (C) (D)
• I recommend starting with (A)
• I use (D) because I find the emacs text editor to be
very effective (and conducive to good
programming practice)
Getting started, (A)
Getting started, (A), con’t
Press
“Ctrl-8”
to open editor!
Reading in data
If data is already in Stata file format (thanks NBER!), we are all set …
clear
set memory 500m
use “/proj/matt/[Link]”
If data is not in Stata format, then can use insheet for tab-delimited files or infile or infix for
fixed-width files (with or without a data dictionary). Another good option is to use
Stat/Transfer
clear
set memory 500m
insheet using “/proj/matt/cricket/[Link]”, tab
clear
set memory 500m
infix ///
int year 1-4 ///
byte statefip 14-15 ///
byte sex 30 ///
byte hrswork 53-54 ///
long incwage 62-67 ///
using [Link]
insheet data
player year round pick position height weight season salary estimated
Andrew Bogut 2005 1st 1st C 84 245 2006 4340520 0
Marvin Williams 2005 1st 2nd F 81 230 2006 3883560 0
Deron Williams 2005 1st 3rd G 75 210 2006 3487400 0
Chris Paul 2005 1st 4th G 75 175 2006 3144240 0
Raymond Felton 2005 1st 5th G 73 198 2006 2847360 0
Martell Webster 2005 1st 6th G-F 81 210 2006 2586120 0
Charlie Villanueva 2005 1st 7th F 83 240 2006 2360880 0
Channing Frye 2005 1st 8th F-C 83 248 2006 2162880 0
Ike Diogu 2005 1st 9th F 80 250 2006 1988160 0
Andrew Bynum 2005 1st 10th C 84 285 2006 1888680 0
Yaroslav Korolev 2005 1st 12th F 81 203 2006 1704480 0
Sean May 2005 1st 13th F 81 266 2006 1619280 0
Rashad McCants 2005 1st 14th G 76 207 2006 1538400 0
Antoine Wright 2005 1st 15th G-F 79 210 2006 1461360 0
Joey Graham 2005 1st 16th G-F 79 225 2006 1388400 0
Danny Granger 2005 1st 17th F 80 225 2006 1318920 0
Gerald Green 2005 1st 18th F 80 200 2006 1253040 0
Hakim Warrick 2005 1st 19th F 81 219 2006 1196520 0
Julius Hodge 2005 1st 20th G 79 210 2006 1148760 0
infix data
1965025135811090025135801016611001 11003341 000002488000000
1965025135811090025135802015821001 11003102 400002180002000
1965025135811090026589103011621006 13105222 000000030000000
1965025135811090032384303013411006 15007102 400007259005250
1965025135811090025135801016511001 01001341 000005000005000
1965025135811090025135802015521001 10003102 400004200004200
1965024645911090024645901015611001 11003102 540004500004500
1965024645911090024645902015321001 12004311 000000000000000
1965024645911090022282003011811006 14106331 000000000000000
1965025633611090025633601016811001 06002212 000005067004827
1965025633611090025633602016021001 10103311 000000000000000
1965022075111090022075101014712001 10003212 000002100002100
1965022075111090022075102014322001 13005102 232002000002000
Stat/Transfer
Describing and summarizing data
describe
summarize
list in 1/100
list if exptot > 1000000 | paytot > 1000000
summarize exptot paytot, detail
tabulate ctscnhos, missing
tabulate cclabhos, missing
Stata data types
id str7 %9s A.H.A identification
number
reg byte %8.0g region code
stcd byte %8.0g state code
hospno str4 %9s hospital number
ohsurg82 byte %8.0g open heart surgery
nerosurg byte %8.0g neurosurgery
bdtot long %12.0g beds set up
admtot double %10.0g total admissions
ipdtot double %10.0g total inpatient days
Stata data types, con’t
• Good programming practices:
– Choose the right data type for your data (“admissions” is a
double?)
– Choose good variable names (“state_code”, “beds_total”,
“region_code”)
– Make the values intuitive (open heart surgery should be 0/1 dummy
variable, not either 1 or 5, where 5 means “hospital performs open
heart surgery”)
• Stata details:
– String data types can be up to 244 characters (why 244?)
– Decimal variables are “float” by default, NOT “double”
• “float” variables have ~7 decimal places of accuracy while “double”
variables have ~15 decimal places of accuracy (floats are 4 bytes of
data, doubles are 8 bytes of data).
• When is this important? MLE, GMM. Variables that are used as
“tolerances” of search routines should always be double. We will
revisit this in lecture 3. In general, though, this distinction is not
important.
• If you are paranoid (like me!), can place “set type double,
permanently” at top of your file and all decimals will be “double” by
default (instead of “float”)
Summarizing data
• Why only 6420 observations for “fyr” variable? 0 observations for
“id” variable?
• Are there any missing “id” variables? How could we tell?
• How many observations are in the data set?
Missing data in Stata
(Disclaimer: In my opinion, this is one the worst “features” of Stata. It is counter-intuitive and
error-prone. But if you use Stata you are stuck with their bad programming language design. So
learn the details!)
• Missing values in Stata
– Missing numeric values are represented as a “.” (a period). Missing string values are “” (an
empty string of length 0)
– Best way to think about “.” value: it is “+/- infinity” (it is an unattainably large or an
unattainably small number that no valid real number can equal).
generate c = log(0) produces only missing values
– What might be wrong with following code?
drop if weeks_worked < 40
regress log_wages is_female is_black age education_years
• Missing values in Stata, new “feature” starting in version 9.1: 27 missing values!
– Now missing values can be “.”, “.a”, … , “.z”
– If “.” is infinity, then “.a” is infinity+1
– For example, to drop ALL possible missing values, you need to write code like this:
drop if age >= .
– Cannot be sure in recent data sets (especially government data sets that feel the need to
use new programming features) that “drop if age == .” will drop ALL missing age
values
– Best programming practice (in Stata 9):
drop if missing(age)
Detailed data summaries
clear
set mem 100m
set obs 50000
generate normal = invnormal(uniform())
generate ttail30 = invttail(30, uniform())
generate ttailX = invttail(5+floor(25*uniform()), uniform())
summ normal ttail* , detail
leptokurtic
distribution!
Tabulating data
clear
set obs 1000
generate c =
log(floor(10*uniform()))
tabulate c, missing
Two-way tables
clear
set obs 10000
generate rand = uniform()
generate cos = round( cos(0.25 * _pi * ceil(16 * rand)), 0.0001)
generate sin = round( sin(0.25 * _pi * ceil(16 * rand)), 0.0001)
tabulate cos sin, missing
Presenting data graphically
• Type “help twoway” to see what Stata has
built-in
– Scatterplot
– Line plot (connected and unconnceted)
– Histogram
– Kernel density
– Bar plot
– Range plot
Preparing data for analysis
• Key commands:
– generate – encode
– replace – assert
– if, in – count
– sort, gsort – append
– merge – collapse
– reshape – strfun
– by • length
– egen • lower
• count • proper
• diff • real
• group • regexm, regexr
• max • strpos
• mean • subinstr
• median • substr
• min • trim
• mode • upper
• pctile
• rank
• sd
• rowmean, rowmax, rowmin
De-meaning variables
clear
set obs 1000
generate variable = log(floor(10*uniform()))
summ variable
replace variable = variable - r(mean)
summ variable
NOTE: “infinity” – r(mean) = “infinity”
De-meaning variables, con’t
clear
set obs 1000
generate variable = log(floor(10*uniform()))
egen variable_mean = mean(variable)
replace variable = variable - variable_mean
summ variable
if/in commands
clear
set obs 50000
generate normal = invnormal(uniform())
list in 1/5
list in -5/-1
generate two_sigma_event = 0
replace two_sigma_event = 1 if (abs(normal)>2.00)
tabulate two_sigma_event
egen commands
Calculate denominator of logit log-likelihood function …
egen double denom = sum(exp(theta))
Calculate 90-10 log-income ratio …
egen inc90 = pctile(inc), p(90)
egen inc10 = pctile(inc), p(10)
gen log_90_10 = log(inc90) – log(inc10)
Create state id from 1..50 (why would we do this?) …
egen group_id = group(state_string)
Make sure all income sources are non-missing …
egen any_income_missing = rowmiss(inc*)
replace any_income_mising = (any_income_missing > 0)
by, sort, gsort
clear
set obs 1000
** randomly generate states (1-50)
gen state = 1+floor(uniform() * 50)
** randomly generate income N(10000,100)
** for each person
gen income = 10000+100*invnormal(uniform())
** GOAL: list top 5 states by income
** and top 5 states by population
sort state
by state: egen mean_state_income =mean(income)
by state: gen state_pop = _N
by state: keep if _n == 1
gsort -mean_state_income
list state mean_state_income state_pop in 1/5
gsort -state_pop
list state mean_state_income state_pop in 1/5
** make state population data file (only for 45 states!)
clear
set obs 45
egen state = fill(1 2)
gen state_population = 1000000*invttail(5,uniform())
save state_populations.dta
list in 1/5
** make state income data file (for all 50 states!)
clear
set obs 1000
gen state = 1+floor(uniform() * 50)
gen income = 10000 + 100*invnormal(uniform())
sort state
merge
command
save state_income.dta
list in 1/5
** created merged data set
clear NOTE:
use state_populations _merge==1, obs only in master
sort state _merge==2, obs only in using
save state_populations, replace
_merge==3, obs in both
clear
use state_income
sort state
merge state using state_populations.dta, uniqusing
tab _merge, missing
tab state if _merge == 2
keep if _merge == 3
drop _merge
save state_merged.dta
merge command, con’t
. ** make state population data file (only for 45 states!)
. Clear
. set obs 45 .
obs was 0, now 45 . ** created merged data set
. egen state = fill(1 2) . clear
. gen state_population = 1000000*invttail(5,uniform()) . use state_populations
. save state_populations.dta . sort state
file state_populations.dta saved . save state_populations, replace
. list in 1/5 file state_populations.dta saved
+-------------------+ . clear
| state state_p~n | . use state_income
|-------------------| . sort state
1. | 1 -4682021 |
2. | 2 1271717 | . merge state using state_populations.dta, uniqusing
3. | 3 -527176.7 | variable state does not uniquely identify observations
4. | 4 907596.9 | in the master data
5. | 5 1379361 |
+-------------------+ . tab _merge, missing
. _merge | Freq. Percent Cum.
. ** make state income data file (for all 50 states!) ------------+-----------------------------------
. clear 1 | 108 10.80 10.80
. set obs 1000 3 | 892 89.20 100.00
obs was 0, now 1000 ------------+-----------------------------------
. gen state = 1+floor(uniform() * 50) Total | 1,000 100.00
. gen income = 10000 + 100*invnormal(uniform())
. sort state . tab state if _merge == 2
. save state_income.dta no observations
file state_income.dta saved . keep if _merge == 3
. list in 1/5 (108 observations deleted)
+------------------+ . drop _merge
| state income | . save state_merged.dta
|------------------| file state_merged.dta saved
1. | 1 10056.04 | .
2. | 1 9999.274 | end of do-file
3. | 1 10042.95 |
4. | 1 10095.03 |
5. | 1 9913.146 |
+------------------+
reshape command
clear
set obs 1000
gen player = 1+floor(uniform() * 100)
bysort player: gen tournament = _n
gen score1 = floor(68 + invnormal(uniform()))
gen score2 = floor(68 + invnormal(uniform()))
gen score3 = floor(68 + invnormal(uniform()))
gen score4 = floor(68 + invnormal(uniform()))
list in 1/3
reshape long score, i(player tournament) j(round)
list in 1/12
reshape command, con’t
String functions (time permitting)
• Stata has support for basic string operations (length,
lowercase, trim, replace).
– Type “help strfun”
• Here is a small example using regular expressions. This
is fairly advanced but can be very helpful sometimes.
– Here is the data set …
game price sectionrow
1 90 FB4,r5
1 75 FB4-5
1 90 4-5
1 80 5FB12
2 80 Field Box 4,12
2 60 4FieldBox12
2 50 Field Box 12, Row 17
2 90 Field Box 2, Rw 5
– GOAL: Get section number and row number for each
observation
Regular expressions
clear
insheet using [Link]
replace sectionrow = subinstr(sectionrow," ", "", .)
local regex = "^([a-zA-Z,.-]*)([0-9]+)([a-zA-Z,.-]*)([0-9]+)$"
gen section = regexs(2) if regexm(sectionrow, "`regex'")
gen row = regexs(4) if regexm(sectionrow, "`regex'")
list
Back to the “standard RA project”
Recall the steps …
1. Read in data
2. Effectively summarize/tabulate data, present graphs
3. Prepare data set for analysis (generate, reshape, parse,
encode, recode)
4. Preliminary regressions and output results
To do (4) we will go through a motivating example …
QUESTION: What is the effect of winning the coin toss on
the probability of winning a cricket match?
Data
team1 team2 toss choice outcome result date_str
West Indies A Priestley's XI A Priestley's XI decided to bat West Indies by 3 wickets 1897_2_15
Trinidad A Priestley's XI Trinidad decided to bat Trinidad by 10 wickets 1897_2_19
Trinidad A Priestley's XI Trinidad decided to bat Trinidad by 8 wickets 1897_2_25
Barbados A Priestley's XI A Priestley's XI decided to bat Barbados by an innings and 411897_1_13
runs
Barbados A Priestley's XI A Priestley's XI decided to field A Priestley's XI by 3 wickets 1897_1_18
Barbados A Priestley's XI A Priestley's XI decided to field Barbados by 136 runs 1897_1_21
Jamaica A Priestley's XI Jamaica decided to bat A Priestley's XI by an innings and 3 1897_3_13
runs
Jamaica A Priestley's XI A Priestley's XI decided to bat A Priestley's XI by 10 wickets 1897_3_16
Jamaica A Priestley's XI A Priestley's XI decided to bat A Priestley's XI by an innings and 151897_3_27
runs
New South Wales A Shaw's XI A Shaw's XI decided to field A Shaw's XI by 9 wickets 1886_12_10
New South Wales A Shaw's XI A Shaw's XI decided to field New South Wales by 122 runs 1887_2_18
New South Wales A Shaw's XI A Shaw's XI decided to bat A Shaw's XI by an innings and 371885_1_24
runs
Victoria A Shaw's XI Victoria decided to bat A Shaw's XI by 9 wickets 1887_3_4
Victoria A Shaw's XI A Shaw's XI decided to bat A Shaw's XI by 118 runs 1884_11_14
New South Wales A Shaw's XI A Shaw's XI decided to bat New South Wales by 6 wickets 1886_11_19
New South Wales A Shaw's XI New South Wales decided to bat A Shaw's XI by 4 wickets 1884_11_21
Victoria A Shrewsbury's XI Victoria decided to bat A Shrewsbury's XI by an innings and 456
1887_12_16
runs
New South Wales A Shrewsbury's XI New South Wales decided to bat A Shrewsbury's XI by 10 wickets 1887_12_9
New South Wales A Shrewsbury's XI New South Wales decided to bat New South Wales by 153 runs 1888_1_13
New South Wales A Shrewsbury's XI New South Wales decided to field New South Wales by 10 wickets 1887_11_10
Trinidad AB St Hill's XII Trinidad decided to bat Trinidad by 146 runs 1901_1_10
Victoria AC MacLaren's XI Victoria decided to bat AC MacLaren's XI by 8 wickets 1902_2_22
New South Wales AC MacLaren's XI New South Wales decided to bat AC MacLaren's XI by an innings and 128
1902_2_31
runs
South Australia AC MacLaren's XI South Australia decided to bat AC MacLaren's XI by 6 wickets 1902_3_14
Basic regressions and tables
set mem 500m
insheet using [Link], tab names
gen year = real(substr(date_str, 1, 4))
assert(year>1880 & year<2007 & floor(year)==year)
gen won_toss = (toss == team1)
gen won_match = (outcome == team1)
summarize won_toss won_match
encode team1, gen(team_id)
Basic
regressions
and tables
Global and local variables
• Scalar variables in Stata can be either local or global. Only
difference is that global variables are visible outside the current
DO file
• Syntax:
clear
local local_variable1 = "local variable"
local local_variable2 = 14170
global global_variable1 = "global variable"
global global_variable2 = 14170
local whoa_dude = "$global_variable2"
display "`local_variable1'"
display $global_variable2
display `whoa_dude'
???
Global and local variables, con’t
local var1 = "var3"
global var2 = "var3"
local var3 = 14170
di "`var1'"
di "$var2"
di "``var1''"
di "`$var2'"
NOTE: Last two lines use syntax that is
somewhat common in ADO files written
by Stata Corp.
Control structures (loops)
foreach var of varlist math reading writing history science {
regress `var’ class_size, robust
estimates store `var’
}
est table `var’
forvalues year = 1940(10)2000 {
regress log_wage is_female is_black educ_yrs exper_yrs ///
if year == `year’, robust
estimates store mincer`year’
}
local i = 1
while (`i’ < 100) {
display `i’
local i = `i’ + 1
}
forvalues i = 1/100 {
display `i’
}
Using control structures for
data preparation
EXAMPLE: Find all
1-city layover
flights given data
set of available origin dest carrier
SFO ORD Delta
flights ORD SFO Delta
ORD CMH Delta
CMH ORD Delta
SFO
ORD
CMH ORD RCA Delta
RCA ORD Delta
CHO RCA Delta
RCA RCA CHO Delta
CHO
Using control structures for
data preparation, con’t
LAYOVER BUILDER ALGORITHM
In the raw data, observations are (O, D, C, . , . ) tuple where
O = origin
D = destination
C = carrier string
and last two arguments are missing (but will be the second
carrier and layover city )
FOR each observation i from 1 to N
FOR each observation j from i+1 to N
IF D[i] == O[j] & O[i] != D[j]
CREATE new tuple (O[i], D[j], C[i], C[j], D[i])
Control structures for
Data Preparation
insheet using [Link], tab names
gen carrier2 = ""
gen layover = ""
local numobs = _N
forvalues i = 1/`numobs' {
di "doing observations `i' ..."
forvalues j = 1/`numobs' {
if (dest[`i'] == origin[`j'] & origin[`i'] != dest[`j']) {
**create new observation for layover flight
local newobs = _N + 1
set obs `newobs'
quietly {
replace origin = origin[`i'] if _n == `newobs'
replace dest = dest[`j'] if _n == `newobs'
replace carrier = carrier[`i'] if _n == `newobs'
replace carrier2 = carrier[`j'] if _n == `newobs'
replace layover = dest[`i'] if _n == `newobs'
}
}
}
}
Control structures for
Data Preparation
A diversion: Intro to Algorithms
• The runtime of this algorithm (as written) is O(N2) time, where N is the number of
observations in the data set (*)
• Whether using Matlab or C, the runtime will be asymptotically equivalent. This is
important to keep in mind when thinking about making calls to C. Most of the time
you only get a proportional increase speed, not an asymptotic increase. In some
cases, thinking harder about getting your algorithm to run in O(N*log(N)) time
instead of O(N2) is much more important than making calls to a programming
language with less overhead.
• In order to improve asymptotic runtime efficiency, need better data structures
than just arrays and/or matrices – specifically, need hash tables (**)
• Perl, Java, and C++ all have built-in hash tables. Some C implementations also
available. We will cover this in more detail in Perl class this afternoon.
• With proper functioning hash tables, the runtime in previous algorithm can be
reduced to O(N) (!!)
(*) the runtime is actually O(N3) in Stata but would be O(N2) in a standard Matlab
and/or C implementation for reasons we will discuss in the next lecture
(**) also called “Associative Arrays”, “Lookup Tables” or “Index Files”
Exercises
• Go to following URL:
[Link]
• Download each DO file
– No DTA files! All data files are loaded from the web (see
“help webuse”)
– First run DO file “AS IS” before doing anything else (this will
save data files in your local directory)
• 3 exercises (in increasing difficulty and – in my
opinion – decreasing importance)
A: Preparing a data set, running some preliminary
regressions, and outputting results
B: More with finding layover flights
C: Using regular expressions to parse data