program define siilin, rclass
version 11.0

syntax varlist(max=3 numeric) [if] [in] [fweight/] , [Cluster(varname numeric)] [NOGRaph]

*** Temporary variables
tempvar _use _rrank _rank _temp _sum _dec _concxw _freq _area _wi _wisum _mean _line

*** Values to 1, 2 e 3
tokenize `varlist'
preserve														// preserve original dataset
quietly {
  set more off

  *** Apply "if" and "in" restrictions and remove missings for concentration variable
  gen byte _use = 1 `if' `in'
  keep if _use==1 												// drop observations not to be used in the temporary dataset

  *** Check for variables being missing
  inspect `1'
  if r(N) == 0 {
  	di as error "Error! Rank variable [" "`1'" "] is missing for all observations."
  	error 416
  }
  inspect `2'
  if r(N) == 0 {
  	di as error "Error! Outcome variable [" "`2'" "] is missing for all observations."
  	error 416
  }
  keep if `2' ~=.

  *** Id rank variable
  local rkv `1' 												// rank variable to local macro rkv
  tabulate `rkv'
  local grouped = r(N)==r(r)									// local macro grouped = 1 if no. of observations = no. of groups
  if r(r)>20 {													// too many groups to calculate SII/RII
  	noisily display as error "Error! too many groups in ranking variable"
  	return scalar er = 1
  	error 1001													// program aborts
  }

  if `grouped' == 1 {

  	*** Create midpoints for ranked groups === GROUPED DATA
  	sort `rkv'
  	if "`exp'" == "" gen _freq=1
  	else gen _freq = `exp'
  	sum _freq
  	gen _rfi =_freq / r(sum)
  	gen _rfcum = sum(_rfi)
  	gen _rkmidp =_rfi/2 + _rfcum[_n-1]
  	replace _rkmidp = _rfi/2 in 1
  }
  else {

  	*** Create midpoints for ranked groups === MICRODATA
  	sort `rkv'
  	tab `rkv', g(_group)
  	if "`exp'" == "" gen _freq=1
  	else gen _freq = `exp'
  	sum _freq
  	gen _rfi =_freq / r(sum)
  	gen _rfcum = sum(_rfi)
  	mean _rfcum, over(`rkv')
  	mat means = e(b)
  	gen _rkmidp = .
  	forval i = 1 / `e(N_over)' {
  	  replace _rkmidp = means[1,`i'] if _group`i'==1
  	}
  }

  *** Graph to check for linearity
  if "`nograph'"=="" {
  	regress `2' i.`rkv' [pw=_freq]
  	predict _mean
  	describe, short											// This block adds two observations
  	local nobs2 = r(N)+2										// including values 0 and 1 for
  	local nobs1 = r(N)+1										// the midpoint rank
  	set obs `nobs2'												// so that the graph includes a line that goes
  	replace _rkmidp=0 in `nobs2'								// from zero to one
  	replace _rkmidp=1 in `nobs1'								// end of block
  	regress `2' _rkmidp [pw=_freq]
  	predict _line
  	noisily sc _mean _rkmidp || line _line _rkmidp, xlabel(0 (.1) 1) legend(label(1 "Observed") label(2 "Fitted")) xtitle(Fractional rank) ytitle(Mean)
  }

  *** Regression to estimate SII
  if "`cluster'" == "" capture regress `2' _rkmidp [pw=_freq]
  else capture regress `2' _rkmidp [pw=_freq], vce(cluster `cluster')
  if _rc==0 {
  	return scalar er = 0

  	noisily nlcom SII: _b[_rkmidp], noheader post
  	mat def sii=e(b)
  	mat def siise=e(V)
  	return scalar sii = sii[1,1]
  	return scalar sii_se = sqrt(siise[1,1])

  	*** Estimating RII
  	if "`cluster'" == "" regress `2' _rkmidp [pw=_freq]
  	else regress `2' _rkmidp [pw=_freq], vce(cluster `cluster')
  	noisily nlcom RII: (_b[_cons]+_b[_rkmidp])/_b[_cons], noheader post
  	mat def rii=e(b)
  	mat def riise=e(V)
  	return scalar rii = rii[1,1]
  	return scalar rii_se = sqrt(riise[1,1])
  }
  else {
  	noisily display _newline as error "Error calculating SII/RII" _newline
  	return scalar er = 1
  	error 999
  }
} // end quietly

end // program