/*  Part I: Using Do-files to Manage & Analyze Data*/

use http://www.ats.ucla.edu/stat/stata/stat130/hsberr

save hsberr, replace

/* begin hsbtest.do
clear
use `1'
nmissing
summarize
assert id>=1 & id<=200
assert (gender>=1 & gender<=2)     | gender==.
assert (race>=1 & race<=4)         | race==.
assert (ses>=1 & ses<=3)           | ses==.
assert (schtyp>=1 & schtyp<=2)     | schtyp==.
assert (prog>=1 & prog<=3)         | prog==.
assert (read>=1   & read<=100)     | read==.
assert (write>=1   & write<=100)   | write==.
assert (math>=1    & math<=100)    | math==.
assert (science>=1 & science<=100) | science==.
assert (socst>=1   & socst<=100)   | socst==. 
/* end hsbtest.do */

/* begin hsbfix */
clear
use `1'
replace id=193     if id==1193
replace read=47    if read==147
replace science=61 if science==-61
replace gender=.   if gender==5
replace race=.     if race<1 | race>4
replace ses=.      if ses<1 | ses>3
replace schtyp=.   if schtyp<1 | schtyp>2
replace prog=.     if prog<1 | prog>3
generate female = gender
recode female 1=0 2=1
label define fem 0 "male" 1 "female"
label value female fem
save `1'2, replace 
/* end hsbfix */

do hsbtest hsberr

do hsbfix hsberr

do hsbtest hsberr2 

/* begin hsberr2an.do */
log using hsberr2.log, text replace
use hsberr2, clear 
summarize read write math science 
univar read write math science 
tabstat write, stat(n mean sd p25 p50 p75) by(female)
tabstat write, stat(n mean sd p25 p50 p75) by(prog)
ttest write, by(gender)
hist write, normal start(30) width(5)
more
kdensity write, normal
more
kdensity math, normal
more
pnorm write
more
qnorm write 
more
tab1 female ses prog
tab prog ses, all 
correlate write read science female
more
graph matrix read science write, half 
regress write read female
rvfplot
more
indexplot rstu, scatter
more
indexplot cooksd, scatter mlabel(id)
log close
/* end hsberr2an */ 

do hsberr2an.do

/*  Part II: Programing an Ado-file */

/* Median Program -- Version #1    */
/*    basic program with one variable        */

program define median1
  version 7
  sort `1'
  quietly count if `1' ~= .
  local n = r(N)
  local mid = int(`n'/2)
  local odd = mod(`n',2)
  
  if `odd' {
    local median = `1'[`mid'+1]
  }
  else {
    local median = (`1'[`mid'] + `1'[`mid'+1])/2
  }
  
  display "Median of `1' = `median'"
end 

/* end median1.ado  */

median1 write

/* Median Program -- Version #2                        */
 /*   multiple variables and saves results in return list     */

program define median2, rclass
  version 7
  syntax varlist
  
  preserve /* preserve data */
  
  display
  display as text " Variable        N    Median"
  display as text "----------------------------"
  
  foreach var of local varlist {
    quietly count if `var' ~= .
    local n = r(N)
    local mid = int(`n'/2)
    local odd = mod(`n',2)
    sort `var'
    
    if `odd' {
      local median = `var'[`mid'+1]
    }
    else {
      local median = (`var'[`mid'] + `var'[`mid'+1]) / 2
    }
    display as res %9s "`var'" %9.0f `n' %10.2f `median'
  }
  return local Mdn = `median'
  return local N = `n'
  
  restore /* restore data */
end 

/* end median2.ado  */

median2 read write math

/* Median Program -- Version #3    */ 
/*   Allows for if and in                           */

program define median3, rclass
  version 7
  syntax varlist [if] [in]        
  
  display
  display as text " Variable        N    Median"
  display as text "----------------------------"
  
  foreach var of local varlist {
    preserve 
    if "`if'"~="" | "`in'"~="" {
      quietly keep `if'  `in' 
    }
    quietly keep if ~missing(`var')
    quietly count 
    local n = r(N)
    local mid = int(`n'/2)
    local odd = mod(`n',2)
    sort `var'
    
    if `odd' {
      local median = `var'[`mid'+1]
    }
    else {
      local median = (`var'[`mid'] + `var'[`mid'+1])/2
    }
    display as res %9s "`var'" %9.0f `n' %10.2f `median'
    restore
  }

  return local Mdn = `median'
  return local N = `n'
end
/* end median3.ado  */

median3 read write math if gender==1

/* Median Program -- Version #4        */ 
/*   Adds option for mean              */

program define median4, rclass
  version 7
  syntax varlist [if] [in], [mean]       
  
  display
  if "`mean'" == "" {
    display as text " Variable        N    Median"
    display as text "----------------------------"
  }
  else {
    display as text " Variable        N    Median      Mean"
    display as text "--------------------------------------"
  }
  
  foreach var of local varlist {
    preserve 
    if "`if'"~="" | "`in'"~="" {
      quietly keep `if'  `in' 
    }
    quietly keep if ~missing(`var')
    quietly count 
    local n = r(N)
    local mid = int(`n'/2)
    local odd = mod(`n',2)
    sort `var'
    
    if `odd' {
      local median = `var'[`mid'+1]
    }
    else {
      local median = (`var'[`mid'] + `var'[`mid'+1])/2
    }
    if "`mean'" == "" {
      display as res %9s "`var'" %9.0f `n' %10.2f `median'
    }
    else {
      quietly summarize `var'
      local mean = r(mean)
      display as res %9s "`var'" %9.0f `n' %10.2f `median' %10.2f `mean'
    }
    restore
  }

  return local Mdn = `median'
  return local N = `n'
end
/*  end median4.ado  */

median4 read write math if gender==1

median4 read write math if gender==1, mean

/*  Part III: Matrix Programming  */

/*  begin matreg.ado  */
program define matreg
  version 6.0
 
    syntax varlist(min=2 numeric) [if] [in] [, Level(integer $S_level)]
    marksample touse                       /* mark cases in the sample */
    tokenize "`varlist'"

    quietly matrix accum sscp = `varlist' if `touse'
    local nobs = r(N)
    local df = `nobs' - (rowsof(sscp) - 1) /* df residual */

    matrix XX = sscp[2...,2...]            /* X'X */
    matrix Xy = sscp[1,2...]               /* X'y */

    matrix b = Xy * syminv(XX)             /* (X'X)-1X'y */
    local k = colsof(b)                    /* number of coefs */
    matrix hat = Xy * b'
    matrix V = syminv(XX) * (sscp[1,1] - hat[1,1])/`df'
    matrix C = corr(V)
    matrix seb = vecdiag(V)
    matrix seb = seb[1, 1...]
    matrix t = J(1,`k',0)
    matrix p = t

    local i = 1
    while `i' <= `k' {
      matrix seb[1,`i'] = sqrt(seb[1,`i'])
      matrix t[1,`i'] = b[1,`i']/seb[1,`i']
      matrix p[1,`i'] = tprob(`df',t[1,`i'])
      local i = `i' + 1
    }

  display
  display "Dependent variable: `1'"
  display
  display "Regression coefficients"
  matrix list b
  display
  display "Standard error of coefficients"
  matrix list seb
  display
  display "Values of t"
  matrix list t
  display
  display "P values for t"
  matrix list p
  display
  display "Covariance of the regression coefficients"
  matrix list V
  display
  display "Correlation of the regression coefficients"
  matrix list C

  matrix drop sscp XX Xy b hat V seb t p C 
end 
/*  end matreg.ado  */

use http://www.ats.ucla.edu/stat/stata/notes/hsb2 

regress write read female 

matreg write read female 

/*  beging matreg2.ado  */
program define matreg2, eclass
  version 6.0
 
    syntax varlist(min=2 numeric) [if] [in] [, Level(integer $S_level)]
    marksample touse                       /* mark cases in the sample */
    tokenize "`varlist'"

    quietly matrix accum sscp = `varlist' if `touse'
    local nobs = r(N)
    local df = `nobs' - (rowsof(sscp) - 1) /* df residual */

    matrix XX = sscp[2...,2...]            /* X'X */
    matrix Xy = sscp[1,2...]               /* X'y */

    matrix b = Xy * syminv(XX)             /* (X'X)-1X'y */
    local k = colsof(b)                    /* number of coefs */
    matrix hat = Xy * b'
    matrix V = syminv(XX) * (sscp[1,1] - hat[1,1])/`df'
   
    estimates post b V, dof(`df') obs(`nobs') depname(`1') /*
      */ esample(`touse')
    est local depvar "`1'"
    est local cmd "matreg"
    
    display
    estimates display, level(`level')
 
  matrix drop sscp XX Xy hat
end 
/*  end matreg2.ado  */

regress write read female 

matreg2 write read female 

vce 

vce, corr 

test female