version 10.0
capture log close
log using chap04, replace
use food, clear

* add observation 
edit
set obs 41
replace income=20 in 41

* estimate regression
regress food_exp income
predict yhat
predict ehat, residuals
predict sef, stdf

* compute t-critical value
scalar define tc = invttail(38,.025)
di "t critical value 97.5 percentile = "  tc
gen lb = yhat - tc*sef
gen ub = yhat + tc*sef
list income lb yhat ub in 41
drop in 41

* R2
pwcorr food_exp income yhat

* examine distribution of food_exp
summarize food_exp, detail

* create $ income
gen inc_dollar = income*100
reg  food_exp inc_dollar

* create reciprocal of income and regress
gen rincome = 1/income
reg food_exp rincome

summarize income
return list
scalar xbar = r(mean)
di xbar
scalar xbar2 = xbar^2
lincom -rincome/xbar2

predict ryhat
twoway (scatter food_exp income) (connected ryhat income)
more
twoway (scatter food_exp income) (connected yhat income) (connected ryhat income)
more

* create log(income) and regress
gen lincome = log(income)
reg food_exp lincome
lincom lincome/xbar
predict linlog
twoway (scatter food_exp income) (connected yhat income) (connected linlog income)
more

* analyze residuals from original equation
histogram ehat, percent
more

summarize ehat, detail
return list
scalar jb = (r(N)/6)*( (r(skewness)^2) + ((r(kurtosis)-3)^2)/4 )
di "Jarque-Bera Statistic = " jb
scalar chic = invchi2tail(2,.05)
di "Chi-square(2) 95th percentile = " chic

scalar pvalue = chi2tail(2,jb)
di "Jarque-Bera p-value = " pvalue

*------------------------
* another example
*------------------------

use wa-wheat, clear

* examine data
summarize
rename  greenough yield
twoway (scatter yield time)
more

* regression
reg yield time
predict yhat

* analyze residuals
predict ehat, residuals
histogram ehat, percent
more

summarize ehat, detail
scalar jb = (r(N)/6)*( (r(skewness)^2) + ((r(kurtosis)-3)^2)/4 )
di "Jarque-Bera Statistic = " jb
scalar chic = invchi2tail(2,.05)
di "Chi-square(2) 95th percentile = " chic
scalar pvalue = chi2tail(2,jb)
di "Jarque-Bera p-value = " pvalue
twoway (scatter yield time) (connect yhat time)
more
twoway (scatter yield time) (lfit yield time)
more
twoway connected ehat time, yline(0)
more
rvpplot time, recast(bar) yline(0)
more

*------------------------
* cubic equation for yield
*------------------------

generate timecube = (time^3)/1000000
reg yield timecube
predict yhat2

* analyze residuals
predict ehat2, residuals
histogram ehat2, percent
more
summarize ehat2, detail
scalar jb = (r(N)/6)*( (r(skewness)^2) + ((r(kurtosis)-3)^2)/4 )
di "Jarque-Bera Statistic = " jb
scalar pvalue = chi2tail(2,jb)
di "Jarque-Bera p-value = " pvalue

twoway (scatter yield time) (connect yhat2 time)
more
rvpplot time, recast(bar) yline(0)
more

*------------------------
* Wage equation
*------------------------

use cps_small, clear

* summarize and plot
summarize
twoway (scatter wage educ) 
more

* add one observation
edit 
set obs 1001
replace educ=12 in 1001

* create log(wage) and plot
gen lwage = log(wage)
twoway (scatter lwage educ) 
more

* log-linear regression
reg lwage educ
predict lwagehat
predict ehat, residuals
predict sef, stdf

* calculate sigma-hat^2
ereturn list
scalar sig2 = e(rss)/e(df_r)
di "sigma-hat squared = " sig2

* Analyze resdiduals

histogram ehat, percent
more
summarize ehat, detail
scalar jb = (r(N)/6)*( (r(skewness)^2) + ((r(kurtosis)-3)^2)/4 )
di "Jarque-Bera Statistic = " jb
scalar chic = invchi2tail(2,.05)
di "Chi-square(2) 95th percentile = " chic
scalar pvalue = chi2tail(2,jb)
di "Jarque-Bera p-value = " pvalue
rvpplot educ, yline(0)
more

* compute natural and corrected predictor and plot
gen yhatn = exp(lwagehat)
gen yhatc = yhatn*exp(sig2/2)
twoway (scatter wage educ, sort) (connected yhatn educ, sort) (connected yhatc educ, sort)
more

* list predicted values
list educ yhatn yhatc in 1001
summarize wage if educ==12 in 1/1000

* R^2
correlate wage yhatn yhatc

* prediction interval
scalar tc = invttail(998,.025)
gen lb_lwage = lwagehat - tc*sef
gen ub_lwage = lwagehat + tc*sef
gen lb_wage = exp(lb_lwage)
gen ub_wage = exp(ub_lwage)

* list and plot

list lb_wage ub_wage in 1001
twoway (scatter wage educ, sort) (connect ub_wage educ, sort) (connect lb_wage educ, sort) (connect yhatn educ, sort)
more

log close
translate chap04.smcl chap04.txt, replace