*! version 3.0   03/26/92  (STB10: sbe9)
program define brier2
        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 VB EB

        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)
                gen `EB'=`2'*(1-`2')
                summ `EB' if `INSAMP'
                local ebrier = _result(3)
                gen `VB' = `2'*(1-`2')*((1-(2*`2'))^2)
                summ `VB' if `INSAMP'
                local nsq = _result(1)^2
                replace `VB'=sum(`VB')
                local vbrier=`VB'[_N]/`nsq'
                local spieg = (`brier'-`ebrier')/(`vbrier'^.5)
                qui ranksum2 `2' if `INSAMP', by(`1')
                local ROC = $S_6
*$mwu/($nn * $mm)
                local ZROC = $S_4
                qui su `1' if `INSAMP'
                local meano=_result(3)
                qui su `2' if `INSAMP'
                local meanf=_result(3)
                qui corr `1' `2' if `INSAMP'
                local corra=_result(4)

                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 bl "Mean probability of forecast" in gr _col(30) %7.4f `meanf' /*
 */_skip(3) in bl "of outcome" _col(45) in gr %7.4f `meano'
        di in bl "Correlation" _col(30) in gr %7.4f `corra'
        di ""
        di in gr "Brier score" _col(30) in ye %7.4f `brier'
        di in bl "Spiegelhalter's z-statistic" _col(30) in gr %7.4f `spieg' /*
 */in bl "  p = "in gr _col(50) %7.4f 1-normprob(`spieg')
        di in bl "ROC area" _col(30) in gr %7.4f `ROC' in bl "  p(H0<=.5)=" /*
 */in gr _col(50) %7.4f normprob(`ZROC')
        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