clear
set more off
set scheme sj

sjlog using plreg_ex1, replace
use plreg_example
plreg tc wage pcap puc kwh life lf kmwire, nlf(cust) generate(func) bw(1) nodraw
sjlog close, replace

sjlog using plreg_ex2, replace
regress tc cust custsq wage pcap puc kwh life lf kmwire
sjlog close, replace

use plreg_example, clear
summarize tc, meanonly
local mean_tc = r(mean)
summarize cust, meanonly
local mean_cust = r(mean)
summarize custsq, meanonly
local mean_custsq = r(mean)
regress tc cust custsq wage pcap puc kwh life lf kmwire
generate quad_func = `mean_tc'+_b[cust]*(cust-`mean_cust')+_b[custsq]*(custsq-`mean_custsq')
plreg tc wage pcap puc kwh life lf kmwire, nlf(cust) generate(func) bw(1) nodraw
matrix coeff = e(b)
matrix score xb_diff = coeff
summarize xb_diff, meanonly
local mean_xb_diff=r(mean)
generate tc_hat = tc-(xb_diff-`mean_xb_diff') // note normalization
twoway scatter tc_hat cust, msymbol(Oh) || ///
line func quad_func cust, lpattern(solid shortdash) ///
ytitle(Log total cost per year) ylabel(5(0.2)6) ///
xtitle(Log customers) ///
legend(order(2 3) label(2 "Nonparametric function") label(3 "Quadratic function")) 
graph export plreg1.eps, replace