drop _all
set mem 60m
sysuse auto 

gen wt = weight / 1000

gen wt2  = wt^2
gen mpg2 = mpg^2

logistic for mpg mpg2 wt wt2

sum wt mpg

drop _all
range wt   1.76  4.84 10
range mpg 12.00 41.00 10
fillin wt mpg

gen wt2  = wt^2
gen mpg2 = mpg^2

predict prob_foreign

xtile qt = prob_foreign , nq(4)

keep wt mpg qt

local dx = (1/2) * (4.84 - 1.76) / 9
local dy = (1/2) * (  41 -   12) / 9

gen id = _n

expand 6
sort id

by id: replace wt  = wt  - `dx' if _n==1
by id: replace mpg = mpg - `dy' if _n<=2

by id: replace wt  = wt  + `dx' if _n==2 | _n==3

by id: replace mpg = mpg + `dy' if _n==3 | _n==4

by id: replace wt  = wt  - `dx' if _n==4

by id: replace wt  = wt  - `dx' if _n==5
by id: replace mpg = mpg - `dy' if _n==5

by id: replace wt  = .          if _n==6
by id: replace mpg = .          if _n==6

twoway area mpg wt if qt==1 , nodropbase bc(navy)  ||		///
       area mpg wt if qt==2 , nodropbase bc(green) ||		///
       area mpg wt if qt==3 , nodropbase bc(red)   ||		///
       area mpg wt if qt==4 , nodropbase bc(yellow) ,		///
       ylabel(20 30 40) xlabel(2 3 4)				///
       plotregion(margin(zero))					///
       legend(subtitle(Predicted probability foreign)		///
              order(1 "(0,.25]"  2 "(.25,.5]" 			///
	            3 "(.5,.75]" 4 "(.75,1)" )			///
       )

twoway area mpg wt if qt==1 , nodropbase bc(navy)  ||		///
       area mpg wt if qt==2 , nodropbase bc(green) ||		///
       area mpg wt if qt==3 , nodropbase bc(red)   ||		///
       area mpg wt if qt==4 , nodropbase bc(yellow) ,		///
       ylabel(20 30 40) xlabel(2 3 4)				///
       plotregion(margin(zero))					///
       legend(subtitle(Predicted probability foreign)		///
              order(1 "( 0,.25]"  2 "(.25,.5]" 			///
	            3 "(.5,.75]"  4 "(.75,1)" )			///
       )							///
       , aspect(1)