
use [DF_7],clear


capture program drop plot_mgr_w162_densities
program define plot_mgr_w162_densities
    version 16.0
    syntax , MGRID(string) ///
        [ PRIORvar(real 2.43) PRIorsd(real -1) ///
          NGRID(integer 2001) XRANGE(numlist min=2 max=2) BW(real -1) ]

    preserve
        quietly {
            * Keep one manager (mgrID is string)
            keep if mgrID == `"`mgrid'"'
			
			
local MID   = mgrID[1]
local MAP   = MAP_w162[1]
local HDIL  = HDI_lower_w162[1]
local HDIU  = HDI_upper_w162[1]
local DMU   = data_162[1]
local DSE   = data_162se[1]

			
			quietly summarize MAP_w162, meanonly
scalar MAP = r(mean)

quietly summarize HDI_lower_w162, meanonly
scalar L = r(mean)

quietly summarize HDI_upper_w162, meanonly
scalar U = r(mean)

quietly summarize data_162, meanonly
scalar MU = r(mean)

quietly summarize data_162se, meanonly
scalar SE = r(mean)

			
            count
            local Ndraws = r(N)
            if `Ndraws' == 0 {
                di as err "No observations found for mgrID = `mgrid'"
                restore
                exit 2000
            }

            * Data mean and SE
            summarize data_162, meanonly
            scalar MU = r(mean)

            summarize data_162se, meanonly
            scalar SE = r(mean)

            if missing(MU) | missing(SE) {
                di as err "data_162 and/or data_162se missing for mgrID = `mgrid'"
                restore
                exit 198
            }
            if (SE <= 0) {
                di as err "data_162se must be > 0 for mgrID = `mgrid'"
                restore
                exit 198
            }

            * Prior SD: default interpret N(0,2.43) (can be modified if desired!)
            scalar PRIOR_SD = cond(`priorsd' > 0, `priorsd', sqrt(`priorvar'))

            *  x-range
            summarize mw162, meanonly
            scalar draw_min = r(min)
            scalar draw_max = r(max)
            scalar draw_sd  = r(sd)

            if "`xrange'" != "" {
                scalar LO = real(word("`xrange'",1))
                scalar HI = real(word("`xrange'",2))
            }
            else {
                scalar SPAN = 6 * max(PRIOR_SD, SE, draw_sd)
               scalar LO = -10
scalar HI = 15
            }

            * Bandwidth 
            local bwopt ""
            if (`bw' > 0) {
                local bwopt "bwidth(`bw')"
            }
            else {
                summarize mw162, detail
                scalar __sd  = r(sd)
                scalar __iqr = r(p75) - r(p25)
                scalar __s0  = min(__sd, __iqr/1.34)
                scalar __h   = 0.9 * __s0 * (`Ndraws')^(-0.2)
                if (!missing(__h) & __h > 0) local bwopt "bwidth(`=__h')"
            }

            * Append grid rows at bottom so draws stay intact
            local N0 = _N
            set obs `=`N0' + `ngrid''

            gen byte isgrid = 0
            replace isgrid = 1 in `=`N0'+1'/`=`N0'+`ngrid''

            gen double x = .
            replace x = LO + (HI-LO)*(_n-(`N0'+1))/(`ngrid'-1) if isgrid

            gen double prior = .
            gen double like  = .
            replace prior = normalden(x, 0, PRIOR_SD) if isgrid
            replace like  = normalden(x, MU, SE)      if isgrid

            * IMPORTANT: do NOT pre-create post; let kdensity create it
            capture drop post
        }

        * kdensity 
        capture noisily kdensity mw162, at(x) generate(post) nograph bw(.35)
        if _rc {
            * fallback syntax for some builds
            capture drop post
            kdensity mw162, at(x) gen(post) nograph `bwopt'
        }

        keep if isgrid

   twoway ///
    (line prior x, lwidth(medthick) lcolor(red) lpattern(dash)) ///
    (line like  x, lwidth(medthick) lcolor(blue) lpattern(dash)) ///
    (line post  x, lwidth(medthick) lcolor(black) lpattern(solid)) ///
    , ///
    legend(off) ///
    graphregion(color(white)) plotregion(color(white)) ///
    xlabel(-10(5)15, labsize(8pt) nogrid) ///
    ylabel(, nogrid) ///
    ytitle("") yscale(off) xtitle("") ///
    xsize(5.5) ysize(2.5)
		
		
		noi di as txt "mgrID: " as res "`MID'"
noi di as txt "MAP_w162: " as res %9.3f `MAP'
noi di as txt "HDI95:   " as res %9.3f `HDIL' as txt " , " as res %9.3f `HDIU'
noi di as txt "data_162: " as res %9.3f `DMU'  as txt " , " as res %9.3f `DSE'
noi di as txt "prior: N(0, 2.43)"



    restore
end


plot_mgr_w162_densities, mgrid([SELECT DESIRED MGR ID])


* Vector exports (pick what you need)
graph export "bayes_illustrate.pdf", replace      
graph export "bayes_illustrate.eps", replace      
graph export "bayes_illustrate.svg", replace      
graph export "bayes_illustrate.emf", replace      
restore

