*! version 1  22Nov2000 by Duolao Wang       (STB-59: sbe42)
* hernes -- a command that fits the Hernes (ASR 1972) model of
* proportion of a cohort every-married by age.
* The command takes two arguments: "age" for age in years and
* "rate" observed age-specific cumulative proportions married

program define hernes,rclass
	version 6
	syntax varlist(min=2 max=2) [if] [in] [, method(string)]
	tokenize `varlist'
	local Age `1'
	local P `2'
	tempvar touse
	mark `touse' `if' `in'
	qui keep if `touse'
	gen t=`Age'-`Age'[1]

	tempname n m K 
	scalar `n'=_N
	scalar `K'=`n'-int(`n'/3)*3
	preserve

	tempvar logg s1logg1 s2logg2 s3logg3
	
	qui  {
	summarize `P'
	if r(M) == 0 { noisily error 2000}
	if r(min) == 0 {
		di in red `"Cumulative Marriage Rate Is Zero"'
		exit 499
		}
	return clear
	}

	* transform data according to formula (13)
	gen `logg'=log(`P'/(1 - `P'))		 

	* expand if the P-vector does not have a length evenly divisible by 3

	qui expand 3 if _n==`n' &`K'==1
	qui expand 2 if _n==`n' &`K'==2
	scalar `m'=_N/3	

	* calculate sums of each third
	gen `s1logg1'=sum(`logg'*(_n<=`m'))
	gen `s2logg2'=sum(`logg'*((_n>=(`m' + 1))&(_n<=(2 * `m'))))
	gen `s3logg3'=sum(`logg'*((_n>=(2 * `m' + 1))&(_n<=_N)))


	tempname s1logg s2logg s3logg bT b loga logk k a A F0

	scalar `s1logg'=`s1logg1'[_N]
	scalar `s2logg'=`s2logg2'[_N]
	scalar `s3logg'=`s3logg3'[_N]

	* get b from formula (14)
	scalar `bT'=(`s3logg' - `s2logg')/(`s2logg' - `s1logg')
	scalar `b'=`bT'^(1/`m')
  
	* get log(a) from formula (15)
	scalar `loga'=((`s2logg' - `s1logg') * (`b' - 1))/(`bT' - 1)^2

	* get log(k) from formula (16)
	scalar `logk'=(1/`m') * (`s1logg' - ((`bT' - 1)/(`b' - 1) * `loga'))
	scalar `k'=exp(`logk')

	* get A ("average initial marriageability")

	** since log.a = A / log(b) from equation (6)
	scalar `a'=exp(`loga')
	scalar `A'=`loga' * log(`b')	
	scalar `F0'= 1/(1 + (1/`k'))


	**method=hernes
	if "`method'"=="hernes"|"`method'"=="he"|"`method'"=="HE"|"`method'"=="Hernes"|"`method'"=="HERNES" {
	restore

	gen rate_hat= 1/(1 + (1/(`k' * `a'^(`b'^(t)))))	
	gen diff=`P'-rate_hat
	format rate_hat  diff %9.3f
	label var rate_hat "Estimated Cumulative First Marriage Rate"
	label var diff "Difference between Observed and Estimated Rate"
	drop t 

	return scalar  A      = `A'
	return scalar  a      = `a'
	return scalar  b      = `b'
	return scalar  k      = `k'
	return scalar  F0     = `F0'
	exit
	}


	**method=nl
	if "`method'"==""|"`method'"=="nl"|"`method'"=="NL" {

	global a0=`a'
	global b0=`b'
	global k0=`k'

	restore,preserve
	nl hernes `P'
	restore

	gen rate_hat= 1/(1 + (1/($k * $a^($b^(t)))))	
	gen diff=`P'-rate_hat
	format rate_hat  diff %9.3f
	label var rate_hat "Estimated Cumulative First Marriage Rate"
	label var diff "Difference between Observed and Estimated Rate"
	drop t

	return scalar  A      = log($a)*log($b)
	return scalar  a      = $a
	return scalar  b      = $b
	return scalar  k      = $k
	return scalar  F0     = 1/(1 + (1/$k))
	exit
	}

end