*! 1.00 Ariel Linden 14Dec2016

program drmmws_bs, rclass

version 13

	syntax varlist [if] [in] [,						///
			Ovars(string) Pvars(string)				///
			NSTRata(string) 						///
			COMMon									///
			ATT										///
			Family(string)							///
			Link(string)							///
			MEDian									///
			SEED(string)							///
			REPS(string)]							//


		quietly {

			local orig `0'
			tokenize `varlist'
			local outcome `1'
			macro shift
			local treat `1'
			macro shift
			local preds `*'

			marksample touse
			count if `touse'
			if r(N) == 0 error 2000
			local N = r(N)
			replace `touse' = -`touse'


			// Give values to options that were not set in the outcome or pscore models
			if "`ovars'" == "" {
				local ovars `preds'
			}
			if "`pvars'" == "" {
				local pvars `preds'
			}
			if "`family'" == "" {
				local family gaussian
			}
			if "`link'" == "" {
				local link   identity
			}

			// Fit propensity score model
			tempvar pscore

			tabulate `treat' if `touse'
				if r(r) != 2 {
				di as err "With a binary treatment, `treat' must have exactly two values (coded 0 or 1)."
			exit 420
			}
			else if r(r) == 2 {
				capture assert inlist(`treat', 0, 1) if `touse'
				if _rc {
				di as err "With a binary treatment, `treat' must be coded as either 0 or 1."
			exit 450
				}
			}
			logit `treat' `pvars' if `touse'
			predict `pscore' if `touse'


			mmws `treat' if `touse', pscore(`pscore') nstrata(`nstrata') `att' `common' replace


			// Fit outcome model
			tempvar pom1 pom0 pomdiff

			// median treatment effects (MTE)
			if "`median'" != "" {

				qreg `outcome' `ovars' [pw = _mmws] if `treat' == 1 & `touse'

				if "`common'" != "" {
					predict `pom1' if _support==1 & `touse'
				}
				else {
					predict `pom1' if `touse'
				}

				sum `pom1', meanonly
				return scalar poms1 = r(mean)

				qreg `outcome' `ovars' [pw = _mmws] if `treat' == 0 & `touse'

				if "`common'" != "" {
					predict `pom0' if _support==1 & `touse'
				}
				else {
					predict `pom0' if `touse'
				}
			} // end median

			// average treatment effects
			else {
				glm `outcome' `ovars' [pw = _mmws] if `treat' == 1 & `touse', link(`link') family(`family')

				if "`common'" != "" {
					predict `pom1' if _support==1 & `touse'
				}
				else {
					predict `pom1' if `touse'
				}

				sum `pom1', meanonly
				return scalar poms1 = r(mean)

				glm `outcome' `ovars' [pw = _mmws] if `treat' == 0 & `touse', link(`link') family(`family')

				if "`common'" != "" {
					predict `pom0' if _support==1 & `touse'
					}
				else {
					predict `pom0' if `touse'
				}
			} // end else (glm)

			sum `pom0', meanonly
			return scalar poms0 = r(mean)

			gen `pomdiff' = `pom1' - `pom0'

			if "`att'" != "" {
				sum `pomdiff' if `treat' == 1, meanonly
			}
			else {
				sum `pomdiff', meanonly
			}
			return scalar drmmws = r(mean)

		} //end quietly

end