MSimulate Command

CLT is a powerful property

* Stata has a wonderfully effective simulate function that allows users to easily simulate data and analysis in a very rapid fashion.

* The only drawback is that when you run it, it will replace the data in memory with the simulated data automatically.

* Which is not a big problem if you stick a preserve in front of your simulate command.

* However, you may want to run sequential simulates and keep the data form all of the simulations together rather than only temporarily accessed.

* Fortunately we can accomplished this task by writing a small program.

cap program drop msim
program define msim

  * Gettoken will split the arguments fed into msim into those before colon and those after.
  gettoken before after : 0 , parse(":")

  * I really like this feature of Stata!

  * First let's strip the colon.  The 1 is important since we want to make sure to only remove the first colon.
  local simulation = subinstr("`after'", ":", "", 1)

  * Now what I propose is that the argument in `before' is used as an extension for names of variables created by the simulate command.

  * First let's save the current data set.


  * Generate an id that we will later use to merge in more data
  cap gen id = _n

  * Save the current data to a temporary location
  tempfile tempsave
  save `tempsave'

  * Now we run the simulation which wipes out the current data.
  `simulation'

  * First we will rename all of the variables to have an extension equal to the first argument
  foreach v of varlist * {
    cap rename `v' `v'`before'
  }

  * Now we need to generate the ID to merge into

  cap gen id = _n

  merge 1:1 id using `tempsave'
  * Get rid of the _merge variable generated from the above command.
  drop _merge
end


* Let's write a nice little program that we would like to simulate.
cap program drop simCLT
program define simCLT

  clear
  set obs `1'
  * 1 is defined as the first argument of the program sim

  * Let's say we would like to see how many observations we need for the central limit theorem (CLT) to make the means of a bernoulli distribution look normal.  Remember, so long as the mean and variance is defined the generally central limit theorem will eventually force any random distribution of means to approximate a normal distribution as the number of observations gets large.
  gen x = rbinomial(1,.25)

  sum x
end

* So let's see first how the simulate command works initially

simulate, rep(200): simCLT 100

* The simulate command will automatically save the returns from the sum command as variables (at least in version 12)

hist mean, kden
* The mean is looking good but not normal

* Now normally what we need to do would be to run simulate again with a different argument.

* But instead let's try our new command with 200!


* But instead let's try our new command!
clear
* Clear out the old results
msim 100: simulate, rep(200): simCLT 100

msim 200: simulate, rep(200): simCLT 200

* Looks good!

msim 400: simulate, rep(200): simCLT 400

msim 1000: simulate, rep(200): simCLT 1000

msim 10000: simulate, rep(200): simCLT 10000

msim 100000: simulate, rep(200): simCLT 100000

msim 1000000: simulate, rep(200): simCLT 100000

* The next two commands can take a little while.
msim 10000000: simulate, rep(200): simCLT 1000000

msim 100000000: simulate, rep(200): simCLT 10000000

sum mean*
* We can see that the standard deviations are getting smaller with a larger sample size.

* How is the histograms looking?
foreach v in 100 200 400 100 1000 10000 100000 1000000 10000000 {
  hist mean`v', name(s`v', replace)  nodraw  title(`v') kden

}

graph combine s100 s200 s400 s1000 s1000 s10000 s100000 s1000000 s10000000 ///
  , title("CLT between 100 and 10,000,000 observations")

* We can see that the distribution of means approximates the normal distribution as the number of draws in each sample gets large.

* This is one of the fundamental findings of statistics and pretty awesome if you think about it.

Non-Parametric PDF Fit Test

* This is an idea that I decided to explore before inspecting how others have addressed the problem.

* As noted by my previous post we cannot use standard independence based draw reasoning in order to test model fit.

* The following command will simulate random draws in the distribution being tested and see how closely they fit to exact pdf draws.

* It will then use that information to see if we can reject the null that the observed distribution is the same as the null distribution.

cap: program drop dist_check
program define dist_check, rclass

  * The syntax command parses the input into useful macros used for later calculations.
  syntax  varlist, Tgenerate(numlist >= 100) [Dist(string)]
 
 
  if "`dist'"=="" local dist="normal"
 
  if "`dist'"=="normal" di "Testing if `varlist' is normally distributed"
  if "`dist'"=="uniform" di "Testing if `varlist' is uniformly distributed"
  if "`dist'"=="poisson" di "Testing if `varlist' has a poisson distributed"

  di "Simulate `tgenerate' draws"
 
  * Construct a t-vector in Mata to store estimation results.
  mata: tdist = J(`tgenerate', 1, .)
 
  *******
  * This will randomly draw from the distribution being tested a number of times equal to tgenerate.
  *******
 
  preserve
  qui: drop if `varlist' == .
 
  forv i=1(1)`tgenerate' {
 
    tempvar x zx z p d

    if "`dist'"=="normal" qui: gen `x' = rnormal()

if "`dist'"=="poisson" {
      qui: sum `varlist'
      qui: gen `x' = rpoisson(r(mean))
    }

if "`dist'"=="uniform" qui: gen `x' = runiform()

    sort `x'
    qui: gen `p' = (_n-.5)/_N

if "`dist'"=="normal" {
      qui: egen `zx' = std(`x')
      qui: gen `z' = invnormal(`p')
      qui: gen `d' = (`z'-`zx')^2
   }

if "`dist'"=="poisson" {
      qui: sum `x'
      qui: gen `z' = round(invpoisson(r(mean), 1-`p'))
 * The invpoisson distribution in Stata is misspecified.  1-p is neccessary.
      qui: gen `d' = (`z'-`x')^2
    }

if "`dist'"=="uniform" {
      qui: gen `z' = _n/_N
      qui: gen `d' = (`z'-`x')^2
    }
   
qui: reg `d'
    local t = _b[_cons]/_se[_cons]
    mata: tdist[`i', 1]=`t'

drop `x' `z' `p' `d'
cap drop `zx'
  }
 
  * From the above loop we should have a vector of t values.
  * We can use that vector to construct a Confidence Interval used observed when true t values as cutoff points.
 
  mata: tsorted = sort(tdist, 1)

  * One tailed test
  local c90 = floor(`tgenerate'*.90)+1
  mata: st_local("t90", strofreal(tsorted[`c90']))
 
  local c95 = floor(`tgenerate'*.95)+1
  mata: st_local("t95", strofreal(tsorted[`c95']))

  local c99 = floor(`tgenerate'*.99)+1
  mata: st_local("t99", strofreal(tsorted[`c99']))
 
  * Two Tailed
    * 90% CI
  local c90U = floor((`tgenerate'+.5)*.95)+1
  mata: st_local("t90U", strofreal(tsorted[`c90U']))
 
  local c90L = ceil((`tgenerate'+.5)*.05)-1
  mata: st_local("t90L", strofreal(tsorted[`c90L']))
 
    * 95% CI
  local c95U = floor((`tgenerate'+.5)*.975)+1
  mata: st_local("t95U", strofreal(tsorted[`c95U']))

  local c95L = ceil((`tgenerate'+.5)*.025)-1
  mata: st_local("t95L", strofreal(tsorted[`c95L']))

    * 99% CI
  local c99U = floor((`tgenerate'+.5)*.99)+1
  mata: st_local("t99U", strofreal(tsorted[`c99U']))

  local c99L = ceil((`tgenerate'+.5)*.01)-1
  mata: st_local("t99L", strofreal(tsorted[`c99L']))

  ** Now we do the estimation
    tempvar x zx z p d

qui: gen `x' = `varlist'

    sort `x'
    qui: gen `p' = (_n-.5)/_N

  * We transform the data in different ways depending upon what distribution we are assuming.
if "`dist'"=="normal" {
      qui: egen `zx' = std(`x')
      qui: gen `z' = invnormal(`p')
      qui: gen `d' = (`z'-`zx')^2
   }

if "`dist'"=="poisson" {
      qui: sum `x'
      qui: gen `z' = round(invpoisson(r(mean), 1-`p'))
      qui: gen `d' = (`z'-`x')^2
    }

if "`dist'"=="uniform" {
      qui: sum `x'
      qui: gen `z' = (`x'-r(min))/(r(max)-r(min))
      qui: gen `d' = (`z'-`x')^2
    }

* This is the regression of interest.
    qui: reg `d'
    local t = _b[_cons]/_se[_cons]

* Now we compare our t with the ts that were drawn when random variables were drawn from our distribution.
di _newline "Estimated t:  `: di %9.3f `t''" _newline

di "One-way Analysis"
di "  CI   (1%) :    0.000   to `: di %9.3f `t99''"
di "  CI   (5%) :    0.000   to `: di %9.3f `t95''"
di "  CI   (10%):    0.000   to `: di %9.3f `t90''"

    mata: psearch = abs(tsorted:-`t')
mata: a = 1::`tgenerate'
mata: b = psearch :== min(psearch)
    mata: st_local("position", strofreal(sum(b:*a)))
local p1 = 1-`position'/`tgenerate'
local p2 = min(`position'/`tgenerate',1-`position'/`tgenerate')*2

di "One-sided  p:`: di %9.3f `p1''   (`position' out of `tgenerate')"

di _newline "Two-way Analysis"
di "  CI   (1%): `: di %9.3f `t99L''    to `: di %9.3f `t99U''"
di "  CI   (5%): `: di %9.3f `t95L''    to `: di %9.3f `t95U''"
di "  CI  (10%): `: di %9.3f `t90L''    to `: di %9.3f `t90U''"
di "Two-sided p: `: di %9.3f `p2''"

    return scalar p1 = `p1'
    return scalar p2 = `p2'

  restore
end

* Let's see how this works on a sample draw.
clear
set obs 1000
gen x = rnormal()+1
dist_check x, t(100) dist(poisson)
* Interestingly, some data distributions seem to fit to fit "better" pdfs from different distributions.
* Thus the estimated t is smaller for some distributions.
* For this reason I have included a two sided confidence interval.

* The following small program will be used to construct a Monte Carlo simulation to see how well the dist_check command is working.
cap program drop checker
program define checker
  syntax [anything], obs(numlist > 0) rdraw(string) reps(numlist > 0) dist(string)
  clear
  set obs `obs'
  gen x = `rdraw'
  dist_check x, t(`reps') dist(`dist')
end

* Easily generates data
checker , obs(1000) rdraw(runiform()*30) reps(200) dist(normal)

* Now let's see it in action
simulate p1 = r(p1) p2 = r(p2) , reps(100): ///
  checker , obs(50) rdraw(rnormal()*30) reps(100) dist(normal)
* With the above simulation almost by definition we know that it is sized propertly but just as a reference.
 
* We should reject at the 10% level for both one sided and two sided confidence intervals.
gen r1_10 = p1<= .1
gen r2_10 = p2<= .1

sum
* mean r1_10 and  mean r2_10 are rejection rates for the one-tailed and two-tailed test respectively.

* Becasue the distribution is truelly the null we are rejecting the null only 10% of the time at the 10% level.

* This looks great.

* Now let's see about the test's power to reject the null when the true distribution is not the null.

checker , obs(1000) rdraw(rnormal()) reps(100) dist(uniform)
  * This distribution is almost uniform
hist x

simulate p1 = r(p1) p2 = r(p2) , reps(100): ///
  checker , obs(1000) rdraw(rnormal()) reps(100) dist(uniform)
* With the above simulation almost by definition we know that it is sized propertly but just as a reference.
 
gen r1_10 = p1<= .1
gen r2_10 = p2<= .1

sum

checker , obs(100) rdraw(rnormal()+runiform()*50) reps(100) dist(uniform)
  * This distribution is almost uniform
hist x

simulate p1 = r(p1) p2 = r(p2) , reps(100): ///
  checker , obs(100) rdraw(rnormal()+runiform()*50) reps(100) dist(uniform)
* With the above simulation almost by definition we know that it is sized propertly but just as a reference.
 
gen r1_10 = p1<= .1
gen r2_10 = p2<= .1

sum

* I think there might be some slight errors.  I hope this post is useful though I would not recommend using this particular command as there are more effective and better established commands already developed for testing distribution fit assumptions.

Macro Parsing

Do file

* In Stata when programming commands we often want to be able to be able to develop our own syntax for our own commands.

* The command "syntax" is often at providing a method for parsing command inputs into easily managed forms.

* However, sometimes we may want syntax options not permitted by the command "syntax".

* If we were to program our own instrumental variable regression that looked like Stata's these limitations might cause problems.

* myivreg y x1 x2 (w1 w2 w3 = z1 z2 z3)

* This is one method how this kind of parsing could be accomplished using the gettoken command.

* Gettoken parses or breaks a macro when a certain character or set of characters is input.

* Let's see it in action

local example1 abc def ghi,jkl
di "`example1'"

* We can see that the local example1 is a disorganized string

* We can seperate it into before and after the space with the following commands.

gettoken before after: example1

di "`before'"

di "`after'"

* The default parse character is " "

* We can specify an alternative parse character or set of characters using the parse option.

gettoken before after: example1, parse(",")

di "`before'"

di "`after'"

* Notice that the parse character is included in the after parse

* Let's see gettoken in: myivreg y x1 x2 (w1 w2 w3 = z1 z2 z3)

cap program drop myivreg
program define myivreg

  di "The local 0 contains the entire command line: `0'"

  gettoken main options: 0 , parse(",")

  gettoken dep_exog brackets: 0, parse("(")

  gettoken dep_var exog : dep_exog

  di "Dependent Variables: `dep_var'"
  di "Exogenous Variables: `exog'"

  gettoken cntr: brackets, parse(")")

  local cntr=subinstr("`cntr'", "(", "", .)

  di "Inside brackets: `cntr'"

  gettoken endog inst: cntr, parse("=")
  local inst=subinstr("`inst'", "=", "", .)

  di "Endogenous Variables: `endog'"
  di "Instrumental Variables: `inst'"

  * At this point we have all of the main pieces of elements we would need to run our ivreg.

  ivreg `dep_var' `exog' (`endog' = `inst')

end

* Let's generate some dummy data
clear
set obs 100

* This will generate some strange data:
*  x2 -> z3 -> w0 -> z2 -> w1 -> z1 -> w2 -> x1 -> y
* Where -> means causes.  This is not a valid iv regression data
foreach v in  x2 z3 w0 z2 w1 z1 w2 x1 y {
  gen `v' = rnormal() `corr_vars'
  local corr_vars ="+`v'"
}

myivreg y (w0=z1)

myivreg y x1 x2 (w0 w1 w2 = z1 z2 z3)

* Everthing seems to be working well.
myivreg y x1 x2 (w0 w1 w2 = z1 z2 z3), test

* However, adding the option "test" does not cause a problem.  Is this a problem.  Not neccessarily.

* We just have not specified any code to ensure that unused syntax is accounted for.

* Combining the gettoken command with the syntax command can accomplish this.


cap program drop myivreg2
program define myivreg2

  di "The local 0 contains the entire command line: `0'"

  * This is a new bit of code that ensures that no options will be included
  syntax anything

  gettoken dep_exog brackets: anything, parse("(")

  gettoken dep_var exog : dep_exog

  di "Dependent Variables: `dep_var'"
  di "Exogenous Variables: `exog'"

  gettoken cntr: brackets, parse(")")

  local cntr=subinstr("`cntr'", "(", "", .)

  di "Inside brackets: `cntr'"

  gettoken endog inst: cntr, parse("=")
  local inst=subinstr("`inst'", "=", "", .)

  di "Endogenous Variables: `endog'"
  di "Instrumental Variables: `inst'"

  * At this point we have all of the main pieces of elements we would need to run our ivreg.

  ivreg `dep_var' `exog' (`endog' = `inst')

end

myivreg2 y x1 x2 (w0 w1 w2 = z1 z2 z3)

myivreg2 y x1 x2 (w0 w1 w2 = z1 z2 z3), test
* Now does not work.

US Daily Gun Deaths R Animation - Sandy Hook

R script

# Listenning to NPR about gun deaths in the US got me thinking.

# Find the article here:
# http://www.slate.com/articles/news_and_politics/crime/2012/12/gun_death_tally_every_american_gun_death_since_newtown_sandy_hook_shooting.html

# Let's create an animation of gun deaths since Dec 14, 2012
gun.deaths <- getcsv.php="" gun-deaths="" http:="" p="" read.csv="" slate-interactives-prod.elasticbeanstalk.com="">
gun.deaths$victimID

# I first need to change the dates from days to days for later reference
  # First I will make a table of the dates wich will include a count
  deaths.tab = table(gun.deaths$date)
    # Calculate Cumulative Dead
    cum.deaths = deaths.tab
    for (i in 1:(length(cum.deaths)-1)) cum.deaths[i+1] = cum.deaths[i]+deaths.tab[i+1]

  plot(deaths.tab, main="Daily Total US Deaths by Gun", ylab="Death count")

  # Number of days in our data (constantly updated every time we run the code
  ndays = length(deaths.tab)

  # This complicated bit of code will force the dates which are currently factor variables into string variables
  gun.deaths$dates = t(data.frame(lapply(gun.deaths$date, as.character), stringsAsFactors=FALSE))

  # Create an empty factor to be filled
  gun.deaths$day = NA
  # This command loops through all of the days and checks if each individual entry in the data from is from that day.
  # If it is, then it assigns that day to the day entry.
  for (i in 1:ndays) gun.deaths$day[gun.deaths$dates==names(deaths.tab)[i]] = i

# We will cut the data into different age categories

# Some individuals have ages missing.  We will code these as category 0.
  gun.deaths$age[is.na(gun.deaths$age)] <- 999="" p="">
  gun.deaths$age.cat = ""
  gun.deaths$age.cat[gun.deaths$age<12 children="" p="">  gun.deaths$age.cat[gun.deaths$age>11 & gun.deaths$age<20 p="" teens="">  gun.deaths$age.cat[gun.deaths$age>=20 & gun.deaths$age<30 adults20="" p="">  gun.deaths$age.cat[gun.deaths$age>=30 & gun.deaths$age<40 adults30="" p="">  gun.deaths$age.cat[gun.deaths$age>=40 & gun.deaths$age<65 madults="" p="">  gun.deaths$age.cat[gun.deaths$age>65 & gun.deaths$age<999 nbsp="" oadults="" p="">
# Adjust the latitude and logitude variables to account for a rescaling of the graph later
# as well as some noise which will help identify when there are multiple deaths in the same city.
  nll = length(gun.deaths$lng)
  gun.deaths$x = ((gun.deaths$lng+125)/(60))*ndays+rnorm(nll)*.006
  # For the graph that will be produced the 20 year olds have the highest likelihood of death.
  # Thus they will provide the y upper limit.
    ymax = ceiling(sum(gun.deaths$age.cat=="adults20")/50)*50
  gun.deaths$y = ((gun.deaths$lat-24)/(27))*ymax+rnorm(nll)*.06

# Generate subsets of the data.
  children         = gun.deaths[gun.deaths$age.cat == "children",]
  teens            = gun.deaths[gun.deaths$age.cat == "teens",]
  adults20         = gun.deaths[gun.deaths$age.cat == "adults20",]
  adults30         = gun.deaths[gun.deaths$age.cat == "adults30",]
  madults          = gun.deaths[gun.deaths$age.cat == "madults",]
  oadults          = gun.deaths[gun.deaths$age.cat == "oadults",]

# This will count cumulative deaths by data subset
children.tab = table(children$date)
  cum.children = children.tab
  for (i in 1:(length(cum.children)-1)) cum.children[i+1] = cum.children[i]+children.tab[i+1]

teens.tab = table(teens$date)
  cum.teens = teens.tab
  for (i in 1:(length(cum.teens)-1)) cum.teens[i+1] = cum.teens[i]+teens.tab[i+1]

adults20.tab = table(adults20$date)
  cum.adults20 = adults20.tab
  for (i in 1:(length(cum.adults20)-1)) cum.adults20[i+1] = cum.adults20[i]+adults20.tab[i+1]

adults30.tab = table(adults30$date)
  cum.adults30 = adults30.tab
  for (i in 1:(length(cum.adults30)-1)) cum.adults30[i+1] = cum.adults30[i]+adults30.tab[i+1]

madults.tab = table(madults$date)
  cum.madults = madults.tab
  for (i in 1:(length(cum.madults)-1)) cum.madults[i+1] = cum.madults[i]+madults.tab[i+1]

oadults.tab = table(oadults$date)
  cum.oadults = oadults.tab
  for (i in 1:(length(cum.oadults)-1)) cum.oadults[i+1] = cum.oadults[i]+oadults.tab[i+1]

# This counts the total deaths
cum.total = cum.adults20+cum.adults30+cum.teens+cum.children+cum.madults+cum.oadults

# In order to get a map of the US we will need to install the library maps
  library(maps)

# Rescale the US map to fit in our data
  map.us = map('state', plot=F)
  map.us$x = ((map('state', plot=F)$x+125)/(60))*ndays
  map.us$y = ((map('state', plot=F)$y-24)/(27)*ymax)

# Static Plot - Final image
i=ndays
dev.new(width=15, height=10)
plot(cum.adults20, type="n", ylim=c(0,ymax),
     ylab="Cumulative Deaths by Age Group" ,
     main="US gun Deaths Since Dec 14, 2012",  cex.main=1.5,
     xlab=paste(names(deaths.tab)[i],"day",toString(i),
           "   ", toString(deaths.tab[i]), "Killed",
           "  ", toString(cum.deaths[i]), "Dead"))
   
  grid(lwd = 2) # grid only in y-direction

  # Draw the US state map
  lines(map.us, type="l")

  # Place dots for every death
  lines(adults20$x, adults20$y, type="p", col="blue",pch=16, cex=.5)
  lines(adults30$x, adults30$y, type="p", col="purple",pch=16, cex=.5)
  lines(madults$x, madults$y, type="p", col="orange",pch=16, cex=.5)
  lines(oadults$x, oadults$y, type="p", col=colors()[641],pch=16, cex=.5)
  lines(teens$x, teens$y, type="p", col="red",pch=16, cex=.5)
  lines(children$x, children$y, type="p", col=colors()[258],pch=16, cex=.5)

  lines(cum.teens,    type="l", col=gray(.9), lwd=2)
  lines(cum.children, type="l", col=gray(.9), lwd=2)
  lines(cum.adults20, type="l", col=gray(.9), lwd=2)
  lines(cum.adults30, type="l", col=gray(.9), lwd=2)
  lines(cum.madults,  type="l", col=gray(.9), lwd=2)
  lines(cum.oadults,  type="l", col=gray(.9), lwd=2)

  lines(cum.teens,    type="l", col="red", lwd=1, cex=.5)
  lines(cum.children, type="l", col=colors()[258], lwd=1, cex=.5)
  lines(cum.adults20, type="l", col="blue", lwd=1, cex=.5)
  lines(cum.adults30, type="l", col="purple", lwd=1, cex=.5)
  lines(cum.madults,  type="l", col="orange", lwd=1, cex=.5)
  lines(cum.oadults,  type="l", col=colors()[641], lwd=1, cex=.5)

  lines(c(ndays,ndays), c(-15,ymax-27), lwd=2, lty=2)

  legend("topright", "Today", cex=1.5, bty="n")

legend(0, ymax+20, expression(Children, Teens, "Adults 20s", "Adults 30s",
   "Adults 40-65", "Adults 65+"), lty = 1:1, col = c(colors()[258], "red",
   "blue", "purple", "orange", colors()[641]),  adj = c(0, 0.6), ncol = 6,
   cex=.75, lwd=2
   )

# Sequential Graphic - Animation

# Animation package must be installed
library(animation)

deaths.animation = function() {

for (i in c(1:ndays, rep(ndays,10))) {
  # Adding the rep(ndays,10)) causes the animation to wait for the equivalent of 10 days after ending.

plot(cum.adults20, type="n", ylim=c(0,ymax),
     ylab="Cumulative Deaths by Age Group" ,
     main="US gun Deaths Since Dec 14, 2012", cex.main=2,
     xlab=paste(names(deaths.tab)[i],"day",toString(i), "   ", toString(deaths.tab[i]), "Killed", "  ", toString(cum.deaths[i]), "Dead"))
   
  # Draw the US state map
  grid(lwd = 2) # grid only in y-direction

  # Place dots for every death
  lines(map.us, type="l")

  t.adults20 = adults20[adults20$day  t.adults30 = adults30[adults30$day  t.madults  = madults[madults$day  t.oadults = oadults[oadults$day  t.teens = teens[teens$day  t.children = children[children$day
  lines(t.adults20$x, t.adults20$y, type="p", col="blue",pch=16, cex=1.5)
  lines(t.adults30$x, t.adults30$y, type="p", col="purple",pch=16, cex=1.5)
  lines(t.madults$x, t.madults$y, type="p", col="orange",pch=16, cex=1.5)
  lines(t.oadults$x, t.oadults$y, type="p", col=colors()[641],pch=16, cex=1.5)
  lines(t.teens$x, t.teens$y, type="p", col="red",pch=16, cex=1.5)
  lines(t.children$x, t.children$y, type="p", col=colors()[258],pch=16, cex=1.5)

  p.adults20 = adults20[adults20$day==i,]
  p.adults30 = adults30[adults30$day==i,]
  p.madults  = madults[madults$day==i,]
  p.oadults = oadults[oadults$day==i,]
  p.teens = teens[teens$day==i,]
  p.children = children[children$day==i,]

  lines(p.adults20$x, p.adults20$y, type="p", col="blue",pch=1, cex=5)
  lines(p.adults30$x, p.adults30$y, type="p", col="purple",pch=1, cex=5)
  lines(p.madults$x, p.madults$y, type="p", col="orange",pch=1, cex=5)
  lines(p.oadults$x, p.oadults$y, type="p", col=colors()[641],pch=1, cex=5)
  lines(p.teens$x, p.teens$y, type="p", col="red",pch=1, cex=5)
  lines(p.children$x, p.children$y, type="p", col=colors()[258],pch=1, cex=5)

  lines(cum.teens[1:i],    type="b", col=gray(.9), lwd=7)
  lines(cum.children[1:i], type="b", col=gray(.9), lwd=7)
  lines(cum.adults20[1:i], type="b", col=gray(.9), lwd=7)
  lines(cum.adults30[1:i], type="b", col=gray(.9), lwd=7)
  lines(cum.madults[1:i],  type="b", col=gray(.9), lwd=7)
  lines(cum.oadults[1:i],  type="b", col=gray(.9), lwd=7)

  lines(cum.teens[1:i],    type="b", col="red", lwd=2)
  lines(cum.children[1:i], type="b", col=colors()[258], lwd=2)
  lines(cum.adults20[1:i], type="b", col="blue", lwd=2)
  lines(cum.adults30[1:i], type="b", col="purple", lwd=2)
  lines(cum.madults[1:i],  type="b", col="orange", lwd=2)
  lines(cum.oadults[1:i],  type="b", col=colors()[641], lwd=2)


  lines(c(ndays,ndays), c(-15,ymax-27), lwd=2, lty=2)

  legend("topright", "Today", bty="n", cex=2.5)

legend(0, ymax+20, expression(Children, Teens, "Adults 20s", "Adults 30s",
   "Adults 40-65", "Adults 65+"), lty = 1:1, col = c(colors()[258], "red",
   "blue", "purple", "orange", colors()[641]),  adj = c(0, 0.6), ncol = 6,
   cex=1.5, lwd=1
   )
 
  ani.pause()
}
}
# deaths.animation()

ani.options(ani.width=1200, ani.height=900)
saveGIF(deaths.animation())
# In order for the saveGIF function to work you need install Image Magic Display (http://www.imagemagick.org/script/index.php)

# Map Updated Jan-31-2012

(OLD MAP BELOW)

Return Text

do file

* The following program will create a tiny program that will send to the display input text.
cap program drop rt
program define rt
    local caller : di _caller()
    * This will parse the `0' text so that there is a before the colon and an after the colon
capt _on_colon_parse `0'

* This will save the after the colon in the macro s(after)
local command = s(after)

* Now we will display the command input
di as input `"`command'"'

* Next we will run the command.
`command'
end

rt: di "adsf"
rt: clear
rt: set obs 1000
rt: gen a = "324"
rt: sum

* Without rt the following commands will see only the output.
forv i = 1/5 {
 clear
 set obs `i'000
 gen a = rnormal()
 sum
}

* With rt prefix we will see both the command as well as the output.
forv i = 1/5 {
 rt: clear
 rt: set obs `i'000
 rt: gen a = rnormal()
 rt: sum
}

Stata Syntax Command's Syntax

do file

* When programming Stata "programs" (loosely known as functions in R) it is often times very useful to use the "syntax" command to parse arguments.

* This command takes the arguments input into the command and uses them to control the function of the command.

* For example:

cap program drop example1
program define example1, rclass

  syntax [anything]
  * The [] tell stata that this argument is optional
  * The anything tells Stata that this argument can be of any form excluding a comma (,)

  di "`anything'"

  * The local function 1 returns the first argument in the command
  di "`1'"

  * While 2 returns the second
  di "`2'"

end

example1 hello world
example1 *
example1 2123 123213 hello world

* We can see the display changes as a result of how the arguments are structured

* Programmers can force the user to input only inputs of a particular structure.

* For instance:


cap program drop example2
program define example2, rclass

  syntax [varlist]
  * The [] tell stata that this argument is optional
  * The anything tells Stata that this argument can be of any form excluding a comma (,)

  di "`varlist'"

  * The local function 1 returns the first argument in the command
  di "`1'"

  * While 2 returns the second
  di "`2'"

  * Variables targetted in this manner can be modified:
  foreach v of varlist `varlist' {
    rename `v' r_`v'_rename
  }

end

clear

gen hello = "content of variable hello"
gen world = "content of variable world"

example2 hello world
* Which works because we created the variables hello and world

example2 hello world
* This does not work because the above command already renamed the variables hello and world

* There is a lot of trickiness involved in effectively programming the syntax command and unfortunately it can be very frustrating because the only feedback stata gives is "invalid syntax" when something goes wrong.

* Also, unfortunately stata does not have many actual coding examples in which syntax is displayed.

* However, I think the following example will go a long way to answering many questions.

* This is the start of a command that will simulate data and test different estimators on that data.

* It will also set default parameters if the user does not input parameters into the simulation directly.

* Before executing the command, the command will display to the user the parameter set to be estimated

cap program drop example3
program define example3, rclass

*
  syntax [anything] [, NGrp(numlist >0 integer) NInd(numlist >0 integer) ///
                    Rgrp(numlist integer) Udist(string) NEXogenous NENdogenous]
  * The first few letters being capital allows the user to abreviate the option to just those two or three letters.

  * Set defaults if the user has not defined any.
  if "`ngrp'"==""   local ngrp=50
  if "`nind'"==""   local nind=20
  if "`rgrp'"==""   local rgrp=0
  if "`udist'"==""  local udist="rnormal()*5"

  * This sets the local parameter exogenous to be equal to 1 by default unless the user specifies ontherwise
  local exogenous=1
  if "`nexogenous'"!="" local exogenous=0

  local endogenous=1
  if "`nendogenous'"!="" local endogenous=0

  di as text _newline "Simulation Paramaters"
  di "Number of groups: `ngrp'"
  di "Average number of individuals in each group: `nind'"
  di "Group size random (0 no, 1 yes): `rgrp'"
  di "Error distribution: `udist'"

  di _newline _c "Models estimated: "
  if `exogenous' == 1 di _c "EXOGENOUS"
  if `endogenous' == 1 di " ENDOGENOUS"
  if  `exogenous' == 0 & `endogenous' == 0 di "No model estimated!"

  ********************

  * The actual simulation

  ********************

  * End of the command
end

  example3 , ngrp(454)
  example3 , ng(454)
  * See how the abreviation is possible

  example3 , nex  udist(runifrom()*3)
  example3 , nex  udist(runifrom()*3) r(1)

  example3 , nex nen

* It is very useful to check out the "help syntax" file.

User written functions in R

R Script

# Writing your own functions can make your programs work much more efficiently, decreasing the lines of codes required to accomplish the desired results as well as simultaneously reducing the chance of error through repetative code.

# It is easy to take whole blocks of arguments from one function to the next with the ... syntax.  This is a sort of catch all syntax that can take any number of arguments.  Thus if you do not like the default choices for some functions such as the default of paste to use " " to conjoin different pieces then a simple function that I often include is the following:

# Concise paste function that joins elements automatically together.
p = function(...) paste(..., sep="")
  x=21
  p("For example x=",x)
  # This function is particularly useful when using the get or assign commands since they take text identifiers of object names.
  a21 = 230
  get(p("a",x))

# Print is a useful function for returning feeback to the user.
# Unfortunately this function only takes one argument.  By conjoining it with my new paste function I can easily make it take multiple arguments.
# Concise paste and print function
pp = function(...) print(p(...))
  # Print displays text.
  for (i in seq(1,17,3)) pp("i is equal to ",i)

# Round to nearest x
# Functions can also take any number of specific arguments.
# These arguments are either targeted by use of the order arguments are placed in or by specific references.
# If an argument is left blank then the default for that argument is used when available or an error message is returned.
round2 = function(num, x=1) x*round(num/x)
  # This rounding function is programmed to function similar to Stata's round function which is more flexible than the built in round command in R.
  # In R you specify the number of digits to round to.
  # In Stata however you specify a number to round.
  # Thus in Stata round(z,.1) is the same in R as round(z,1).
  # However, the Stata command is more general than the R one since Stata for instance could round to the nearest .25 while the R command would need to be manipulated to accomplish the same goal.
  # I have therefore written round2 which rounds instead to the nearest x.
  round(1.123, 2)
  round2(1.123, .01)
  # Yeild the same result. Yet, round2 will work with the following values.
  # Order is not neccessary.  The following produces identical results.
  round2(x=.01,num=1.123)
  round2(123, 20)
  # The round2 has a default x of 1 so round2 can be used quickly to round to the nearest whole number.
  round2(23.1)
  # The original round has the same default.

# Using modular arithmatic is often times quite helpful in generating data.
# The following function reduces a number to the lowest positive integer in the modular arithmatic.
mod = function(num, modulo) num - floor(num/modulo)*modulo
  # This syntax is programmed in a similar manner to Stata's mod command.
  mod(15,12)
  # This kind of a 12 modular system is the kind of system used for hours.
  # Thus mod(15,12) is the same as
  mod(3,12)
  # Or even:
  mod(-9,12)
  # There is a built in operator that does the same thing as this function.
  -9 %% 12

# Check for duplicate adjacent rows.  First row of a duplicate set is ignored but subsequent rows are not.
radjdup = function(...) c(F, apply(...[2:(dim(...)[1]),]==...[1:dim(...)[1]-1,],1,sum)==dim(...)[2])

# Half the data is duplicated.
  example.data = data.frame(id=1:100, x=mod(1:100, 5), y=rep(1:10,each=10))

# The data needs to be sorted for radjdup to work.
  example.order = example.data[order(example.data$x, example.data$y),]

cbind(example.order, duplicate=radjdup(example.order[,2:3]))

# Half the observations should be duplicated. We must sort our data for the adjected row duplicate command to be of any use with the current data.
  example.data = example.data[order(example.data[,1], example.data[,2]),]

# We expect to see half the observations are duplicate flagged.  This is because the first instance of any duplicate is not flagged.
  cbind(example.order, duplicate=radjdup(example.order[,2:3]))

Easy Monte Carlo Sampler Command

R script

# In R I am often times unsure about the easiest way to run a flexible Monte Carlo simulation.
# Ideally, I would like to be able to feed a data generating function into a command with an estimator and get back out interesting values like mean estimates and the standard deviation of estimates.
# Stata has a very easly command to use called "Simulate".
# I am sure others have created these types of command previous but I don't know of them yet so I will write my own.

MC_easy = function(dgp, estimator, obs, reps, save_data_index=0) {
  # dgp (Data Generating Proccess) is the name of a function that creates the data that will be passed to the estimator
  # estimator is the name of a function that will process the data and return a vector of results to be analyzed
  # obs is the number of observations to be created in each spawned data set
  # reps is the number of times to loop through the dgp and estimation routine
  # save data index is a vector of data repetitions that you would like saved

  # I typically start looping programming by manually going through the first loop value before constructing the loop for both debugging and demonstrative purposes.

  # This command will create a a starting sample data set.
  start_data = get(dgp)(obs)
  # This command runs the estimate values on that starting data set.
  MC_results = get(estimator)(start_data)

  # Create an empty list of save data values
  save_data = list()
  # If the first data set generated has been specified as a data set to save then save it as the first entry in the save_data list.
  if (sum(1==save_data_index)>0) save_data[[paste("1",1,sep="")]] = start_data

  for (i in 2:(reps-1)) {
    temp_data = get(dgp)(obs)
 
    # If this repetition is in the save_data_index then save this data to the list save_data
    if (sum(i==save_data_index)>0) save_data[[paste("i",i,sep="")]] = temp_data
 
    MC_results=rbind(MC_results, get(estimator)(temp_data))
  }
  # Display the results of the estimations
  MC_results = as.data.frame(MC_results)
  print(rbind(mean=mean(MC_results), sd= sd(MC_results), var= sd(MC_results)^2))

  # If the number of items for which the data was saved is equal to zero then only return the results of the estimation.
  if (sum(save_data_index>0)==0) return(MC_results)

  # If the number of items is greater than zero also return the saved data.
  if (sum(save_data_index>0)>0)  return(list(MC_results,save_data))
}

#### Done with the Monte Carlo easy simulation command.

# Let's try see it in action.


# Let's define some sample data generating program
endog_OLS_data = function(nobs) {
  # Let's imagine that we are wondering if OLS is unbiased when we have a causal relationship between an unobserved variable z and the observed variables x2 and x3 and the unobserved error e.
  # x4 is an observed variable that has no relationship with any of the other variables.
  z = rnorm(nobs)
  x1 = rnorm(nobs)
  x2 = rnorm(nobs) + z
  x3 = rnorm(nobs) - z
  x4 = rnorm(nobs)

  e = rnorm(nobs)*2 + z

  y = 3 + x1 + x2 + x3 + 0*x4 + e
  data.frame(y, x1, x2, x3, x4)
}

sample_data = endog_OLS_data(10)
sample_data
# Everything appears to be working well

# Now let's define a sample estimation
OLS_est = function(temp_data) {
  # Summary grabs important values from the lm command and coef grabs the coefficients from the summary command.
  lm_coef = coef(summary(lm(y~x1+x2+x3+x4, data=temp_data)))

  # I want to reduce lm_coef to a single vector with values for each b, each se, and each t
  lm_b = c(lm_coef[,1])
    names(lm_b) = paste("b",0:(length(lm_b)-1),sep="")
  lm_se = c(lm_coef[,2])
    names(lm_se) = paste("se",0:(length(lm_se)-1),sep="")
  lm_se2 = c(lm_coef[,2])^2
    names(lm_se2) = paste("se2_",0:(length(lm_se2)-1),sep="")

  lm_rej = c(lm_coef[,4])<.1
    names(lm_rej) = paste("rej",0:(length(lm_rej)-1),sep="")
  c(lm_b,lm_se, lm_se2,lm_rej)
}
OLS_est(sample_data)

# From this is it not obvious which estimates if any are biased.
# Of course our sample size is only 10 so this is no surprise.
# However, even at such a small sample size we can identify bias if we repeat the estimation many times.

# Let's try our MC_easy command.

MC_est = MC_easy(dgp="endog_OLS_data", estimator="OLS_est", obs=10, reps=500, save_data_index=1:5)
# We can see that our mean estimate on the coefficient is close to 3 which is correct.
# While b1 is close to 1 (unbiased hopefully) while b2 is too large and b3 too small and b4 is close to zero which is the true value.
# What is not possible to observe with 500 repetitions but capable of being observed with 5000 is that the standard error is a downward biased estimate of the standard deviation.
# This is a well known fact.  I think even there is a proof that there is no unbaised estimator for the standard error for the OLS estimator.
# However, se^2 is an unbiased estimator of the variance of the coefficients.
# Though the variance of se^2 is large, making it neccessary to have a large number of repetitions before this becomes apparent.
# The final term of interest is the rejection rates.  We should expect that for the coefficients b0-b3 that the rejections rates should be above random chance 10% (since there trully exists an effect) which is what they all ended up being though with b1 and b3 the rejection rate is only slightly above 20% which indicates that our power given our sample size given the effect size of x1 and x3 is pretty small.
# The final note is that though we do not have much power we are at least not overpowered when it comes to overrejecting the null when in fact the null is true.  For b4 the mean rejection rate is close to 10% which is the correct power.

# We can recreate the results table by accessing elements of the MC_est list returned from the function.
rbind(mean=mean(MC_est[[1]]),sd=sd(MC_est[[1]]),var=sd(MC_est[[1]])^2)

# We may also be curious about the medians and maxes.
summary(MC_est[[1]])

# Or we may like to check if our data was generated correctly
MC_est[[2]]

Estimating Person Characteristics from IRT Data - 3PL Model

Original Code

# One of the basic tasks in item response theory is estimating the ability of a test taker from the responses to a series of items (questions).

# Let's draw the same pool of items that we have used on several previous posts:

# First let's imagine we have a pool of 100 3PL items.
set.seed(101)

npool = 500

pool = as.data.frame(cbind(item=1:npool, a=abs(rnorm(npool)*.25+1.1), b=rnorm(npool), c=abs(rnorm(npool)/7.5+.1)))
summary(pool)

# Drawing on code from a previous post we can calculate several useful functions:

# (http://www.econometricsbysimulation.com/2012/11/matching-item-information.html)

# Each item has a item characteristic curve (ICC) of:
PL3 = function(theta,a, b, c) c+(1-c)*exp(a*(theta-b))/(1+exp(a*(theta-b)))

# Let's imagine that we have a single test taker and that test taker has taken the first 15 items from the pool.

items.count = 15
items.taken = pool[1:items.count,]

# And that the person has a latent theta ability of 1
theta = 1.3

# Let's calculate the cut points for each of the items.
items.cut = PL3(theta, items.taken$a, items.taken$b, items.taken$c)

# We can see how the cut point works by graphing
plot(0,0, type="n", xlim=c(-3,3),ylim=c(0,1), xlab=~theta,
           ylab="Probability of Correct Response", yaxs = "i", xaxs = "i" , main="Item Characteristics Curves and Ability Level")

for(i in 1:items.count) {
  lines(seq(-3,3,.1),PL3(seq(-3,3,.1), items.taken$a[i], items.taken$b[i], items.taken$c[i]), lwd=2)
  abline(h=items.cut[i], col="blue")
}
abline(v=theta,col="red", lwd=3)

# Now let's draw a uniform draw that we will use to calculate whether each item as passed.

rdraw = runif(items.count)

# Finally, we will calculate item responses
item.responses = 0
item.responses = items.cut > rdraw

###############################################
# Done with Simulation - Time for Estimation

# We want to use the information we know about the items (the item parameters and the responses) in order to estimate a best guess at the true ability of the test taker.

# First we must check if the person got all of the items either correct or incorrect.

sum(item.responses)
# If this is either a 0 or a number equal to the number of items then we cannot esimate an interior maximum without additional assumptions.

# We will attempt to recover our theta value using the r command optim

# First we need to specify the function to optimize over.

MLE = function(theta) sum(log((item.responses==T)*PL3(theta, items.taken$a, items.taken$b, items.taken$c) +
                              (item.responses==F)*(1-PL3(theta, items.taken$a, items.taken$b, items.taken$c))))
# The optimization function takes as its argument the choice variables to be optimized (theta).
# The way the above optimization works is that you specify the probility of each response piecewise.
# If the response is correct, then you count the CDF of theta up to that point as contributing to the probability of observing a correct outcome.  If the response is negative, then you count it as contributing the the probability of a incorrect outcome.  You then choose the theta that produces the greatest total probability.

MLEval = 0
theta.range = seq(-3,3,.1)
for(i in 1:length(theta.range)) MLEval[i] =MLE(theta.range[i])

plot(theta.range, MLEval, type="l", main="Maximim Likelihood function", xlab= ~theta, ylab="Sum of Log Likelihood")
abline(v=theta, col="blue")
# We can visually see that the maximum of the slope will not be at the true value though it will be close.

optim(0,MLE, method="Brent", lower=-6, upper=6, control=list(fnscale = -1))
abline(v=optim(0,MLE, method="Brent", lower=-6, upper=6, control=list(fnscale = -1))$par, col="red")

# We can see that we can estimate theta reasonably well with just 15 items from a paper test (red line estimate, blue line true).  However, looking at the graph of the ICCs, we can see that for most of the items, the steepest point (where they have the most discriminating power) is at an ability set lower than the test taker's ability.  Thus, this test provides the most information about a person who has a lower ability than the person with a theta=1.2.

# We can use R's optim function to find the ideal theta that would maximize the information from this test.
# Item information is:
PL3.info = function(theta, a, b, c) a^2 *(PL3(theta,a,b,c)-c)^2/(1-c)^2 * (1-PL3(theta,a,b,c))/PL3(theta,a,b,c)

# Notice, this is not the best way of defining the test information function since the items are not arguments.
test.info = function(theta) sum(PL3.info(theta, items.taken$a, items.taken$b, items.taken$c))

# Construct a vector to hold the test information
info = 0
for(i in 1:length(theta.range)) info[i]=test.info(theta.range[i])

plot(theta.range, info, type="l", main="Information Peaks Slightly Above 0", xlab= ~theta, ylab="Information")
abline(v=theta, col="blue")
# But we want to know about the test taker at theta

optim(0,test.info, method="Brent", lower=-6, upper=6, control=list(fnscale = -1))
# The person this test would be best suited to evaluate would have an ability rating of .19
abline(v=optim(0,test.info, method="Brent", lower=-6, upper=6, control=list(fnscale = -1))$par, col="red")

Stata - Write your own "fast" bootstrap

  * Imagine we would like to bootstrap the standard errors of a command using a bootstrap routine.

  * I created a previous post demonstrating how to write a bootstrap command.  This is a similar post however the bootstrap is much faster than the previous one.
 
  * First let's generate some data.
 
  clear
  set obs 1000
  gen x1 = rnormal()
  gen x2 = 2*rnormal()
  gen u = 5*rnormal()
 
  gen y = 5 + x1 + x2 + u
 
  local dependent_var x1 x2
 
  local command_bs reg y `dependent_var'
 
  * First let's see how the base command works directly.
  `command_bs'
 
  * As a matter of comparison this is the built in bootstrap command.
  bs, rep(100): `command_bs'
 
  * The following code is yet another user written bootstrap alternative.
 
  * I wrote this
 
  * Specify the number of bootstrap iterations
  local bs = 100
 
  * Save the number of observations to be drawn
  local N_obs = _N
 
  * Number of terms plus one for the constant of coefficients to be saved
  local ncol = wordcount("`dependent_var'")+1
 
  mata theta = J(`bs',`ncol',0)
  forv i = 1(1)`bs' {
    * Preserve the initial form of the data
preserve

    * Draw the indicators of the resample
mata draw = ceil(runiform(`N_obs',1):*`N_obs')

* Create an empty matrix to hold the number of items to expand
mata: expander = J(`N_obs',1, 0)

* Count the number of items per observation to generate.
    qui mata: for (i=1 ; i <= `N_obs'; i++) expander[i]=sum(draw:==i) ; "draws complete"

* Pull the expander value into stata
getmata expander=expander

* Drop unnessessary data
qui drop if expander == 0

* Expand data, expand==1 does nothing
    qui expand expander

* Run a regression
    qui `command_bs'

* Send to mata the matrix of results
mata theta[`i',] = st_matrix("e(b)")
   
* Configure the visual display
di _c "."
    if (int(`i'/10)==`i'/10) di _c "|"
    if (int(`i'/50)==`i'/50) di " " `i'

restore
  }
 
  * The estimates of the coefficients have been saved into a theta matrix.
  mata theta
 
  * Now let's calculate the standard deviations.
  mata bs_var = variance(theta)
  mata bs_var
  mata bs_se = diag(bs_var):^.5
  mata bs_se

Maximum Likelihood Using R

# The short tutorial by Professor Marco Steenbergen is very helpful in learning to program maximum likelihood in R.

# http://www.unc.edu/~monogan/computing/r/MLE_in_R.pdf

# I will attempt my first R MLE with a probit.

# Let's first generate our data

# Declare the number of observations to create.
nobs=1000

x1 = rnorm(nobs)*.15^.5
x2 = rnorm(nobs)*.35^.5
u = rnorm(nobs)*(.5)^.5

# I have generated the variables so that the inside of the normal CDF has a standard deviation of 1.  This should help the coefficients be more easily recovered.

inside = 2*x1+2*x2+u

# Using basic
# Var(inside) = (.15*2)^2*var(N(0,1)) + (.35*2)^2*var(N(0,1)) + (.5^.5)^2*var(N(0,1)) = .15+.35+.5 =1

# The probability of have a success is equal the cumulative density function at the inside value.
p = pnorm(inside)

y = rbinom(nobs,1,p)

mydata = data.frame(y=y,x1=x1,x2=x2)

# The estimates we would like to duplicate are:
summary(glm(y~x1+x2, family=binomial(link="probit"), data=mydata))

# First we must specify the function to maximize
# All we need to do is to specify a function that takes choice values and maximizes over


myprobit <- function(theta) {
  beta0 <- theta[1]
  beta1 <- theta[2]
  beta2 <- theta[3]

  # Calculate the log likelihood values for a set of betas
  log.lik <- sum(log(pnorm((-1)^(1-mydata$y)*(beta0 + beta1*mydata$x1 + beta2*mydata$x2))))
    # This might look quite confusing.  However, it basically says that if y=1 then count beta0+beta1*x1+beta2*x2 as positively contributing to the probability.  However, if y=0 then we need to count not the contribution but in effect the negative contribution of the inside term.  That is when y is positive count the inside term as predicting a positive outcome.  When it is negative count the inside term as falsely predicting a positive outcome.  I am not sure how much sense this all makes.  I am still trying to wrap my head around it all.

  # Return that log likelihood value.
  return(log.lik)
}

# Let's test out our likelihood function using a generic set of parameters
myprobit(c(1,1,1))

# Let's see what the probit log likelihoods look like as we vary B1 from -1 to 2
b.range <- seq(-1,2.5,.1)
b0.loglik <- b1.loglik <- b2.loglik <- 0*b.range
for (i in 1:length(b.range)) {
  b0.loglik[i] <- myprobit(c(b.range[i],1,1))
  b1.loglik[i] <- myprobit(c(1,b.range[i],1))
  b2.loglik[i] <- myprobit(c(1,1,b.range[i]))
}

# I will define a quick min max function to help setting up the range of plots
minmax <- function(x) c(min(x),max(x))

plot(minmax(b.range), minmax(c(b0.loglik, b1.loglik, b2.loglik)),
      type="n", main="Log-Likelihood Around the point (1,1,1)",
      ylab = "log likelihoods", xlab= ~ beta )

lines(c(1,1),minmax(c(b0.loglik, b1.loglik, b2.loglik)), col="green", lwd=2)
     
lines(b.range, b0.loglik, col="black", lwd=2)
lines(b.range, b1.loglik, col="blue", lwd=2)
lines(b.range, b2.loglik, col="purple", lwd=2)

# The graph displays how the log likelyhood changes with respect to different coefficients around the point (1,1,1).  If we were using the Newton-Ralphsony algorithm then B0 (green) would step left, while B1 and B2 would step right.  This is because the alorithm use the determinate at the point to approximate the next place to step.

# Looking good.  By specifying the control fnscale as =-1 we turn the optimization problem from the default minimize to maximize.
optim(c(1,1,1), myprobit, control=list(fnscale = -1))
# Great! Though, it is worth noting that the coefficients are not identical.  This is not neccessarily a problem since we would expect different optimization algorithms to produce different results.

# The vector c(1,1,1) is the starting values of for the coefficients to be estimated.

# Notice that this has worked fine for this specific model prob(p)=pnorm(B0+B1*x1+B2*x2).
# Clearly, a more sophisticated level of programming would be neccessary to duplicate the flexibility of the probit command.  That is, a model that allows for more variables to be specified.

Program Your own Bootstrap Command #R#

# A similar bootstrap command is written in Stata: see http://www.econometricsbysimulation.com/2012/07/write-your-own-bootstrap-command.html

# The bootstrap command is an extremely powerful command that resamples from you sampling distribution in order to estimate a standard error for an estimator.  Often, using the bootstrap resampling technique can yeild estimates of standard errors on estimators that are otherwise extremely difficult to estimate (using the Delta method, the predominant method for estimating standard errors for non-linear functions).

# This post will develop a simple function that takes a data source and an estimator and generates standard errors.

# First let's simulate some data.

nobs = 200

x1 = rnorm(nobs)
x2 = rnorm(nobs)/2
u = rnorm(nobs)

y = 10-x1+2*x2+u*2

# Now that we have some vectors let's join those vectors together as a data.frame
mydata = data.frame(y=y, x1=x1, x2=x2)

# A simple OLS regression will generate good standard errors.

summary(lm(y~x1+x2, data=mydata))
# We can see that x2 has a standard deviation of about twice that of x1, this is due to x2 having half the standard deviation of x1.

# As a check, let's say we would like to estimate the standard deviation of our coefficients using bootstrapping instead.

# Let's define our bootstrapping function
data.boot <- function(Input.data, command, reps, est=F, hist=F) {
  # data is the source data that will be passed into the boot strapping command
  # command is the function that will do the operation that we desire and return a vector
  # reps is the number of repetitions in the simulation

  # Let us first get our base estimates
  original.est <- get(command)(Input.data)

  # Calculate the number of estimates needed
  nest = length(original.est)

  # Make a holder matrix for our bootstrap results
  holder <- matrix(NA, ncol = nest , nrow = reps)
    # Note, the length of our original.est is the number of coefficients we will estimate standard errors for.

  # Calculate the number of observations to resample from
  nobs <- nrow(Input.data)

  # Now let's start our bootstrapping resampling loop
  for (i in 1:reps) {
    # Draw nobs uniform integers draws from 1 to nobs
    posdraws <- ceiling(runif(nobs)*nobs)
 
    # Sample a random sample from our original data set
    draw <- Input.data[posdraws,]
      # Now we should have a new data frame called draw which has nobs draws.  Observing draw it should be easy to see there there are some observations which are repeated.
   
    # Apply our command to the new data set and save it as a row in holder
    holder[i,] <- get(command)(draw)
  }

  # We have completed our bootstrap rountine.  Now we need only perform whatever statistics we want on the results.

  # Calculate the standard deviation for each column
  sds <- apply(holder,2,sd)

  if (hist==T) {
    mfrow=c(1,1)
    frame()
    if (nest<= 3) par(mfrow=c(nest,1))
    if ((nest> 3)&(nest<= 6)) par(mfrow=c(3,2))
    for (j in 1:nest) hist(holder[,j], main=paste("Estimates of ", names(original.est)[j]), xlab="")
  }
  print(rbind(original.est,sds))
  return(list(estimates=holder, sds=sds))
}

# We will define our first function that we would like to bootstrap the standard errors
example1 <- function(BSdata) lm(y~x1+x2, data=BSdata)[[1]]
  # The key is that we need the function to return a vector to the boostrap command.  Thus the subscript [[1]] restricts the returned values to only be the estimates.  It my take some manipulation to get a single vector out of an estimate.

example1.res <- data.boot(mydata, "example1", 200, hist=T)
  # This will run a bootstrap of the estimates (200 times) and make a histogram

# We can see that the standard deviation of our bootstrap estimates are similar to those of our linear model.
summary(lm(y~x1+x2, data=mydata))

# Let's generate some more complex data to test our bootstrapping on
x1 = runif(nobs)
x2 = rnorm(nobs)
x3 = exp(rnorm(nobs))
x4 = log(runif(nobs))
u = rnorm(nobs)

y = -10+2*x1+2*x2+.2*x3+1.4*x4+u*3

mydata2 <- data.frame(y=y, x1=x1, x2=x2, x3=x3, x4=x4)

# Let's try something a bit more complex, first we will be estimating more coefficients
example2 <- function(BSdata) {
  # First save the results of the OLS regression into a
  a<- lm(y~x1+x2+x3+x4, data=BSdata)
  # Next save the coefficients
  coef <- a[[1]]
  # Save the r squared
  r2 <- summary(a)$r.squared

  # Specify a vector to return
  return(c(coef,r2=r2))
}
example2(mydata2)

example2.results <- data.boot(mydata, "example2", 200, hist=T)

We can see that we are able to estimate not only standard coefficients but also the standard deviation of any estimator generated from sampling.  The extreme draws in r2 is probably due to some of the sampling distributions having very low values (or high values) of the error term by chance in those extreme draws.

Nested functions in R and Stata

#----------------- Starting in R -------------------------#

# Both R and Stata have built in limit of how many nested functions each package will allow.
a <- 0

recursive_call <- function() {
  a <<- a+1
  print(a)
  recursive_call()
}

recursive_call()
# R gives up at 2494 calls.

#-----------------Switching to Stata---------------------*

* This is an identical program (function) in Stata

global a=0

cap program drop recursive_call
program recursive_call

  global a=$a + 1
  di $a
  recursive_call

end


* Stata on the other hand throws in the towel at 63.  This looks bad but look at the response by Alan Riley!


* The Following Code is code from a email comment by Alan RileyStataCorp LP Vice President, Software Development

/* --------------------------- And In Mata------------------------------ */
set more off // don't pause every screen full of output
mata:

a = 0

function recursive_call() {
external a
a = a+1
printf("%f\n", a)
  recursive_call()
}

recursive_call()

end

/*

Alan: The above will error out at some point when it hits its recursion limit,
and Mata will print out a (long) traceback log intended to help in
debugging when a Mata function hits some system limit. Due to such deep
recursion, this traceback log will be long and so the output from
recursive_call() will scroll off the top of your screen. In any case,
after running the above, you can then re-enter Mata and see what the value
of "a" is:

mata:
a
end



*/

You should see 32766.




As an aside, you will notice that

* I imposed an upper limit of recursive calls at 10000, if left uncontrolled it will keep going till it gives an error. I am impressed with the speed and power of Mata!
recursive_call

Linear models with heteroskedastic errors

* Imagine we have a sample of hotels and prices of rooms in those hotels
clear
set obs 200
set seed 101
* This sets the random seed at 101.  Thus every time this simulation is run it will product identical results.

gen star = ceil(runiform()*7)/2+.5
  label var star "Number of stars of hotel"

gen id = _n

gen het = abs(rnormal())
  label var het "Unobserved heterogeneity."

* Generate a the number of reservations observed per hotel.
gen num_reservations = rpoisson(150)
expand num_reservations

gen seasonality = abs(rnormal())
  label var seasonality "Seasonal demand"

gen v = abs(rnormal())
  label var v "Unobserved variance in the variance term"

gen u = rnormal()*(5+het*15+7.5*star+22.5*seasonality+v*10)
  label var u "Error term"

gen p = 175 + 20*star + 15*seasonality + u + het

sum p
* There are some p values which are less than 0 but we can think of those as special deals, coupons, refunds, or other situations that might result in the effective price being less than 0 dollars.
* If we were to eliminate the less than 0 prices then we would in be enforcing left censoring which is a different problem.  See "tobit".  This blog has several posts touching on the use of the tobit.

* Estimate the price through direct OLS
reg p star seasonality

* Set the panel level observation
xtset id

* Though we have panel data we cannot effectively use fixed effect or random effects approaches to identify the price effect having one more star has on prices.
xtreg p star seasonality, fe
xtreg p star seasonality, re

* Now let's attempt a two step method to more efficient identify the error and reestimate the OLS.

reg p star seasonality
* The OLS regression looks pretty good.  Let's see if we can improve on it.  Note, the 95% confidence interval did not capture the true coefficient of 20 but that is not necessarily a problem.

reg p star seasonality, robust
reg p star seasonality, cluster(id)
* Using robust and cluster robust estimates of the standard error does not change the variance so much that the 95% confidence interval encloses the 20.

predict uhat, resid

gen uhat_abs = abs(uhat)
  label var uhat_abs "Abs of OLS residual"

two  (scatter uhat_abs seasonality) (scatter uhat_abs star), ///
              legend(label(1 "Seasonality") label(2 "Stars"))

* Since E(u)==0 var(u) is equal to u squared Var(u)=(u-E(u))^2
* Likewise, u^2 is approximated by u^2
* Similarly sd(u) should be approximated by (u^2)^.5=abs(u)
reg uhat_abs star seasonality

predict uhat_abs_hat, xb
  * I might be doing this wrong.  I am trying aweights which say something about weighting by the inverse of the variance of the observation.)

gen uhat2 = 1/uhat_abs_hat^2

reg p star seasonality
* Let's see how the unweighted estimate performs

reg p star seasonality [aweight = uhat2]
* Using variance weights does not seem to improve the estimate

* Alternatively we can use the MLE estimator allowing the conditional standard deviation of the error as well as the conditional mean to vary linearly.

cap program drop mle_ols
program mle_ols
  args log_like xb sigma
  qui replace `log_like' = ln(normalden($ML_y1-`xb',0,`sigma'))
end

ml model lf mle_ols (price: p = star seasonality) (sigma: star seasonality)
ml maximize

* Using the MLE estimator we seem to have gained precision in the 3rd decimal place of the coefficients.  The 95% CI still does not enclose the 20 but it is closer.  Still, this is not indicative of a problem.  Perhaps if we simulated this 1000 times and substantially more than 50 times the CI did not enclose the 20 then we might be worried.

3 ways to write a probit command

* Probit is a commonly used command to model binary outcome variables.

* The maximum likelihood sytax in Stata is specific
* I will first define three different equivalent ml programs.  Then we will solve them to see how they perform.

cap program drop myprobit1
program define myprobit1
  * Tell Stata that this code is written in version 11
  version 11
  * Define the arguments input in the maximum likelihood function
  * Stata automatically generates the first argument as a temporary variable name
  * that is used to store the natural log of the likelihood value of the likelihood function.
  args ln_likehood xb
  * The probit tries to directly maximize the probability of seeing the outcome (either y==1 or y==0)
  * If y==1 then we maximize the probability of seeing a positive outcome
  qui replace `ln_likehood' = ln(  normal(`xb')) if $ML_y==1
  * If y==0 then we want to maximize the probability of seeing a negative outcome
  qui replace `ln_likehood' = ln(1-normal(`xb')) if $ML_y==0
end

cap program drop myprobit2
program define myprobit2
  version 11
  args ln_likehood xb
  * Because of the symetry of the normal CDF the following set of expressions is equivalent.
  qui replace `ln_likehood' = ln(normal( `xb' )) if $ML_y==1
  qui replace `ln_likehood' = ln(normal(-`xb' )) if $ML_y==0
end

cap program drop myprobit3
program define myprobit3
  version 11
  args ln_likehood xb
  * This is a somewhat opeque way of writing the same thing as above.
  * I do not prefer it to the first two forms except that it is likely to run slightly faster
  * since it is a single command.
  qui replace `ln_likehood' = ln(normal((-1)^(1-$ML_y)*`xb'))
end


clear
set obs 1000

gen x1 = rnormal()*.5
gen x2 = rnormal()*.5
gen u = rnormal()*(.5^.5)

gen p=normal(x1+x2+u)

gen y=rbinomial(1,p)

probit y x1 x2
* This is the benchmark

ml model lf myprobit1 (y=x1 x2)
ml maximize
* Basically works the same except there is no pseudo r2 reported.

timer on 1
ml model lf myprobit2 (y=x1 x2)
ml maximize
timer off 1

timer on 2
ml model lf myprobit3 (y=x1 x2)
ml maximize
timer off 2

timer list

di "myprobit2 is " r(t1)/r(t2) " times slower than myprobit3"

* The following is a brief exploration attempting to run two probits at once.
* It shows how ml can optimize two equations jointly
cap program drop mytwoprobit1
program define mytwoprobit1
  version 11
  args ln_likehood xb1 xb2
  qui replace `ln_likehood' = ln(normal((-1)^(1-$ML_y1)*`xb1')) ///
                            + ln(normal((-1)^(1-$ML_y2)*`xb2'))
end

* This is equivalent to the previous except that it uses the binormal distribution rather than two seperate normals.  I tried specifying the correlation as a third equation but that did not work.  I will need to play around a bit more to figure out how to program a biprobit.
cap program drop mytwoprobit2
program define mytwoprobit2
  version 11
  args ln_likehood xb1 xb2
  qui replace `ln_likehood' = ln(binormal((-1)^(1-$ML_y1)*`xb1', ///
                                          (-1)^(1-$ML_y2)*`xb2', ///
 0))

end

clear
set obs 1000

gen x1 = rnormal()*.5
gen x2 = rnormal()*.5

gen p1=normal(x1-x2)
gen p2=normal(-x1+x2)

gen y1=rbinomial(1,p1)
gen y2=rbinomial(1,p2)

ml model lf mytwoprobit1 (y1:y1=x1 x2) (y2:y2=x1 x2)
ml maximize

ml model lf mytwoprobit2 (y1:y1=x1 x2) (y2:y2=x1 x2)
ml maximize

Draw n values from a defined set

# This is a simple function (program in Stata) that takes an input vector of values (d) and randomly draws from them a random subset set of values of length n.

# We want a to program our own function to draw x draws from a set of possible draws
pull <- function(n,d) {
  # Create an empty vector to store the output
  v <- c()
  # Loop from 1 to n and draw 1 observation from 1 randomly for each of those
  for(i in 1:n) v=append(v, d[ceiling(runif(1)*length(d))])
    # Okay so I know this looks pretty intractable.  The for loop allows a singe line of code to follow on the same line.  So it is just saying loop through i from 1 to n.  Then the vector v has one value added to the end of it.  That value is a random draw from the vector d.
  # Return the compiled vector of random draws v
  v
}

# Now let's see how well our little code works.
pull(20,c("Agent Bob","Agent Sue","008","Agent Henrietta"))

# This will draw five numbers randomly between 1 and 10
pull(5,1:10)

# One can do this same kind of thing in Stata.  However, it would be much less elegant.

###########
# After writing this function I realized my original solution was not nearly as elegant as it could have been.
###########


# We want a to program our own function to draw x draws from a set of possible draws
pull <- function(n,d) d[ceiling(runif(n)*length(d))]

# Now let's see how well our little code works.
pull(20,c("Agent Bob","Agent Sue","008","Agent Henrietta"))

# This will draw five numbers randomly between 1 and 10
pull(5,1:10)


# This alternative function specification uses the preferred vector forms in R's language.
# To get a hint of what it is doing play around with this.
vect1 <- c("a","b","c","d","e")
# To see vect1 just type
vect1

# To get a single element of vect1
vect1[2]

# To get a subset of vect1
vect1[c(1,2,3)]

# To get a subset with repetition
vect1[c(1,2,1,3,1)]

#Now
ceiling(runif(20)*length(vect1))

# This is a nice way of figuring out how code is working in R.  Simply take it apart and see what it does by itself.