*! version 3.0   03/26/92  (STB10: sbe9)
program define brier
	version 3.0
	local varlist "req min(2) max(2)"
	local if "opt"
	local in "opt"
	local options "Group(integer 10)"
	parse "`*'"
	parse "`varlist'", parse(" ")
	cap assert `1'==0 | `1'==1
	if (_rc) { error 149 }

	tempvar FMD GID DJ INSAMP NGROUP FJ DJ LOG

	quietly {
		sort `2'
		gen byte `INSAMP' = 1 `if' `in'
		replace `INSAMP' = 0 if `1'==. | `2'==.
		gen float `FMD' = sum(`INSAMP')
		local totsamp = `FMD'[_N]
		replace `FMD' = (`FMD'-1)/`FMD'[_N]        /* relative percentile */
		gen int `GID' = int(`group' *  `FMD') + 1  /* group ID */
		sort `GID'
		summ `1' if `INSAMP'
		local dbar = _result(3)
		summ `2' if `INSAMP'
		local fbar = _result(3)
		local fvar = _result(4) * (_result(1)-1)/_result(1)

		summ `2' if `INSAMP' & `1'==1
		local f1 = _result(3)
		replace `FMD' = `f1' if `1'==1
		summ `2' if `INSAMP' & `1'==0
		local f0 = _result(3)
		replace `FMD' = `f0' if `1'==0
		summ `FMD' if `INSAMP'
		local sfmin = _result(4) * (_result(1)-1)/_result(1)

		by `GID': gen long `NGROUP' = sum(`INSAMP')
		by `GID': replace `NGROUP' = `NGROUP'[_N]
		by `GID': gen float `FJ' = sum(`2')/`NGROUP'
		by `GID': replace `FJ' = `FJ'[_N]
		by `GID': gen float `DJ' = sum(`1')/`NGROUP'
		by `GID': replace `DJ' = `DJ'[_N]
		by `GID': gen byte `LOG' = _n == _N

		replace `FMD' = (`2'-`1')^2                /* (f(i)-d(i))^2 */
		summ `FMD' if `INSAMP'
		local brier = _result(3)

		replace `FMD' = (`FJ'-`1')^2
		summ `FMD' if `INSAMP'
		local brier1 = _result(3)

		replace `FMD' = `NGROUP'*`DJ'*(1-`DJ')
		summ `FMD' if `LOG'
		local sanders = _result(1)*_result(3)/`totsamp'

		replace `FMD' = `NGROUP'*(`FJ'-`DJ')^2
		summ `FMD' if `LOG'
		local relinsm = _result(1)*_result(3)/`totsamp'

		local oiv = `dbar'*(1-`dbar')

		replace `FMD' = `NGROUP'*(`DJ'-`dbar')^2
		summ `FMD' if `LOG'
		local murphy = _result(1)*_result(3)/`totsamp'
		local relinla = (`fbar'-`dbar')^2
		local sfd1 = (`f1'-`f0')*`oiv'

		replace `FMD' = (`2'-`FJ')^2
		summ `FMD' if `INSAMP'
		local gerr = _result(3)

	}

	di in gr "Brier score" _col(30) in ye %7.4f `brier'
	di in gr "Sanders-modified Brier score" _col(30) in ye %7.4f `brier1'
	di in gr "Sanders resolution" _col(30) in ye %7.4f `sanders'
	di in gr "Outcome index variance" _col(30) in ye %7.4f `oiv'
	di in gr "Murphy resolution" _col(30) in ye %7.4f `murphy'
	di in gr "Reliability-in-the-small" _col(30) in ye %7.4f `relinsm'
	di in gr "Forecast variance" _col(30) in ye %7.4f `fvar'
	di in gr "Excess forecast variance" _col(30) in ye %7.4f `fvar'-`sfmin'
	di in gr "Minimum forecast variance" _col(30) in ye %7.4f `sfmin'
	di in gr "Reliability-in-the-large" _col(30) in ye %7.4f `relinla'
	di in gr "2*Forecast-Outcome-Covar" _col(30) in ye %7.4f 2*`sfd1'
	di in gr "Grouping forecast error" _col(30) in ye %7.4f `gerr'
/*

	di `brier'-(`oiv'+`fvar'+`relinla'-2*`sfd1')
	di `brier1'-`sanders'-`relinsm'
	di `sanders'-`oiv'+`murphy'
	di `relinsm'-`fvar'-`relinla'+2*`sfd1'-`murphy' - `gerr'

	assert abs(`brier1'-`sanders'-`relinsm')<.0001
	assert abs(`sanders'-`oiv'+`murphy')<.0001
	assert abs(`relinsm'-`fvar'-`relinla'+2*`sfd1'-`murphy')<.0001

The last assertion does not hold because there is a difference between
brier and brier1 conceptions.  Yates statements hold for the original
Brier score, but the Sanders and Murphy decompositions hold for the
brier1 score.  The difference between brier and brier1 is the
grouping error.  It is probably somewhat more complicated than I
computed here.

*/
end