version 15.0
clear all
set more off
set type double
set linesize 160
capture log close _all
cd "C:/Users/cattom/Documents/econometrics/Week7"
log using "week7_stata.log", text replace

display "Week 7: household consumption and disposable income"
display "Input: numeric cells C3:L17; row 17 contains reported totals."
import excel using "Eg2_1_1/Eg2_1_1.xlsx", sheet("Sheet1") cellrange(C3:L17) clear

tempname audit
postfile `audit' double income int n double detail_total reported_total difference using "source_totals_audit.dta", replace
local income = 800
foreach col in C D E F G H I J K L {
    assert `col'[1] == `income'
    quietly summarize `col' in 2/14
    local count = r(N)
    local total = r(sum)
    local reported = `col'[15]
    post `audit' (`income') (`count') (`total') (`reported') (`total' - `reported')
    rename `col' consume`income'
    local income = `income' + 300
}
postclose `audit'
preserve
    use "source_totals_audit.dta", clear
    list, noobs separator(0)
    quietly ds, has(type numeric)
    format `r(varlist)' %21.15g
    export delimited using "source_totals_audit.csv", datafmt replace
restore

keep in 2/14
generate int excel_row = _n + 3
reshape long consume, i(excel_row) j(income)
drop if missing(consume)
sort income excel_row
generate int household_id = _n
order household_id excel_row income consume
assert _N == 98
assert !missing(income, consume)
label variable income "每月可支配收入（元）"
label variable consume "每月消费支出（元）"
label variable excel_row "Source worksheet row"
label variable household_id "Record identifier"
save "households_clean.dta", replace
quietly ds, has(type numeric)
    format `r(varlist)' %21.15g
export delimited using "households_clean.csv", datafmt replace
summarize income consume

regress consume income
estimates store OLS
estimates save "consumption_ols.ster", replace
predict double fitted, xb
predict double residual, residuals
generate double residual_sq = residual^2
scalar slope_ols = _b[income]
scalar intercept_ols = _b[_cons]
scalar rsquared_ols = e(r2)
scalar rmse_ols = e(rmse)
scalar rss_ols = e(rss)
scalar N_ols = e(N)
scalar F_ols = e(F)
scalar df_ols = e(df_r)
scalar pF_ols = Ftail(e(df_m), e(df_r), e(F))

tempname coef
postfile `coef' str12 term double estimate double se double t double p double ci_low double ci_high using "regression_coefficients.dta", replace
foreach term in income _cons {
    scalar tval = _b[`term'] / _se[`term']
    scalar crit = invttail(e(df_r), .025)
    post `coef' ("`term'") (_b[`term']) (_se[`term']) (tval) (2*ttail(e(df_r), abs(tval))) (_b[`term']-crit*_se[`term']) (_b[`term']+crit*_se[`term'])
}
postclose `coef'
preserve
    use "regression_coefficients.dta", clear
    quietly ds, has(type numeric)
    format `r(varlist)' %21.15g
    export delimited using "regression_coefficients.csv", datafmt replace
restore

quietly summarize residual, meanonly
scalar residual_sum = r(sum)
generate double income_residual = income * residual
quietly summarize income_residual, meanonly
scalar income_residual_sum = r(sum)
display "Sum of residuals = " %21.12g scalar(residual_sum)
display "Sum of income * residual = " %21.12g scalar(income_residual_sum)
drop income_residual
save "households_fitted.dta", replace
quietly ds, has(type numeric)
    format `r(varlist)' %21.15g
export delimited using "households_fitted.csv", datafmt replace

preserve
    collapse (count) n=consume (mean) mean_consume=consume (sd) sd_consume=consume (min) min_consume=consume (max) max_consume=consume, by(income)
    recast double mean_consume sd_consume
    generate double fitted_mean = scalar(intercept_ols) + scalar(slope_ols)*income
    generate double mean_residual = mean_consume - fitted_mean
    save "income_group_summary.dta", replace
    quietly ds, has(type numeric)
    format `r(varlist)' %21.15g
    export delimited using "income_group_summary.csv", datafmt replace
    list, noobs separator(0)
restore

set scheme s2color
graph set window fontface "SimSun"
twoway (scatter consume income, mcolor("26 71 111") msymbol(O) msize(small)) (lfit consume income, lcolor("174 44 51") lwidth(medthin)), title("不同可支配收入水平组家庭消费支出的条件分布图", size(medium) color(navy)) xtitle("每月可支配收入（元）", size(medsmall)) ytitle("每月消费支出（元）", size(medsmall)) xlabel(1000(500)3500) ylabel(500(500)3000, angle(90) grid glcolor("225 236 238")) xscale(range(500 4000)) yscale(range(0 3600)) legend(order(1 "consume" 2 "Fitted values") rows(1) region(lcolor(black) fcolor(white))) graphregion(color("234 242 243")) plotregion(color(white)) xsize(9) ysize(6.55) name(consumption_fit, replace)
graph save "consumption_fit.gph", replace
graph export "consumption_fit.png", width(1800) replace
graph export "consumption_fit.svg", replace
graph export "consumption_fit.pdf", replace

* Heteroskedasticity-robust standard errors as a sensitivity check.
regress consume income, vce(robust)
estimates store Robust
estimates save "consumption_robust.ster", replace
tempname robustcoef
postfile `robustcoef' str12 term double estimate double se double t double p double ci_low double ci_high using "robust_coefficients.dta", replace
foreach term in income _cons {
    scalar tval = _b[`term'] / _se[`term']
    scalar crit = invttail(e(df_r), .025)
    post `robustcoef' ("`term'") (_b[`term']) (_se[`term']) (tval) (2*ttail(e(df_r), abs(tval))) (_b[`term']-crit*_se[`term']) (_b[`term']+crit*_se[`term'])
}
postclose `robustcoef'
preserve
    use "robust_coefficients.dta", clear
    quietly ds, has(type numeric)
    format `r(varlist)' %21.15g
    export delimited using "robust_coefficients.csv", datafmt replace
restore
estimates restore OLS
estat hettest
scalar bp_chi2 = r(chi2)
scalar bp_p = r(p)

preserve
    clear
    set obs 1
    generate N = scalar(N_ols)
    generate intercept = scalar(intercept_ols)
    generate slope = scalar(slope_ols)
    generate r_squared = scalar(rsquared_ols)
    generate root_mse = scalar(rmse_ols)
    generate rss = scalar(rss_ols)
    generate F = scalar(F_ols)
    generate df_resid = scalar(df_ols)
    generate p_F = scalar(pF_ols)
    generate residual_sum = scalar(residual_sum)
    generate income_residual_sum = scalar(income_residual_sum)
    generate bp_chi2 = scalar(bp_chi2)
    generate bp_p = scalar(bp_p)
    quietly ds, has(type numeric)
    format `r(varlist)' %21.15g
    export delimited using "model_statistics.csv", datafmt replace
restore

display "WEEK7_ANALYSIS_COMPLETED"
log close
exit, clear