*! Chi2 goodness of fit test; v.2.0 by Stas Kolenikov
program define chi2fit, rclass
  version 6.0
  syntax varlist (numeric min=2 max=2) [if] [in] [fw pw/], [BINS(int 0) Detail DF(int 0)]
* the syntax is: the first variable is the actual data,
* the second variable is the cdf of the hypothetical data

  tokenize `varlist'
  local actual `1'
  local hypoth `2'

  marksample touse
  qui {
*    tempvar hypoth
*    g double `hypoth' = `exp' if `touse'

    tempvar counts
    if "`weight'"~="" { 
                        g `counts'=`exp' if `touse' 
                        sum `counts' if `touse'
                        replace `counts' = `counts'*r(N)/r(sum)
                      }
               else   { g byte `counts' = 1 if `touse' }

    qui sum `counts' if `touse'
    local total = r(sum)
* total number of observations used

    if `total'<32 {
      di in red "Are you sure you want to run this program with so few observations?"
      exit 8
    }

    if `bins'==0    { local bins=int(log(`total')/log(2)+1) }
    if `bins'<=`df'+3 { local bins=`df'+3 }
    if "`detail'"~="" {
      noi di in gre _n "Using " in yel %3.0f `bins' in gre " bins for grouping" _n
      noi di in gre "  Bins |  From  |   To   | Under H0 |  Actual"
      noi di in gre "-------+--------+--------+----------+----------"
    }

    local sorted : sortedby
    if "`sorted'"=="" {
      tempvar tosort
      g long `tosort'=_n
      local sorted `tosort'
    }

    sort `touse' `hypoth'
    cap assert `actual' <= `actual'[_n+1] if `touse'
    if _rc==9 {
       di in red "Are you sure `exp' is a good enough cdf for `actual' ?"
       exit 9
    }
    tempvar binno lower upper Nact Nhyp chi2
    g byte `binno'=int(`hypoth'*`bins'+1) if `touse'

    local i=1
    g double `lower'=.
    g double `upper'=.
    g long `Nact'=.
    g long `Nhyp'=.
    while `i'<=`bins' {
      if `i'==1       {
                        replace `lower'=0             in `i'
                        local lbndry "-infty"
                      }
           else       {
                        replace `lower'=`upper'[_n-1] in `i'
                        local lbndry `rbndry'
                      }
      if `i'==`bins'  {
                        replace `upper'=1             in `i'
                        local rbndry="+infty"
                      }
           else       {
                        sum `hypoth' if `binno'==`i', meanonly
                        local left=r(max)
                        sum `hypoth' if `binno'==`i'+1, meanonly
                        local right=r(min)
                        replace `upper'=.5*(`right'+`left') in `i'
                        sum `actual' if `binno'==`i', meanonly
                        local left=r(max)
                        sum `actual' if `binno'==`i'+1, meanonly
                        local right=r(min)
                        local rbndry=string(.5*(`right'+`left'),"%8.3f")
                      }
      sum `counts' if `binno'==`i'
      cap assert r(N)>5
      if _rc==9 { /* we should rebin our stuff */ }
      replace `Nact'=r(sum) in `i'
      replace `Nhyp'=(`upper'-`lower')*`total' in `i'
      if "`detail'"~="" {
        #delimit ;
         noi di _col(5) %-3.0f in yel `i' in gre _col(8) "| " in yel
                "`lbndry'" _col(17) in gre "| " in yel "`rbndry'"
                _col(26) in gre "| " in yel %7.2f `Nhyp'[`i']
                in gre _col(37) "| " in yel %5.0g `Nact'[`i']
        ;
        #delimit cr
      }
      local i=`i'+1
    } /* end of while bin */
    g double `chi2'=(`Nact'-`Nhyp')^2/`Nhyp' in 1/`bins'

    sum `chi2', meanonly
    return scalar df     = `bins'-1-`df'
    return scalar chi2   = r(sum)
    return scalar pvalue = chiprob(return(df),return(chi2))
    sort `sorted'
  } /* end of quietly */
  if "`detail'"~="" { di in gre "-------+--------+--------+----------+----------" }
  #delimit ;
  di in gre _n "       Goodness of fit chi-square test:" _n
               "       comparing `actual' with `exp'"
     _n(2) _col(18) "Chi2 (" in yel %4.0f return(df)
     in gre ") = " in yel %-9.3f return(chi2)
     _n in gre _col(18) "Prob > Chi2 = " in yel %6.4f return(pvalue)
  ;
  #delimit cr
end

/* History:
  1.2 ???
  2.0 syntax changed, [pw] added

*/