********************************************************************************
*  Empirical Exercise 5.1 - Earnings and Height: Hypothesis Testing
*  Stock & Watson - Introduction to Econometrics (4th Edition)
*
*  Code: Fatih Kansoy - www.fatih.ai
*  Course: Financial Econometrics, SJTU Summer School 2025
********************************************************************************

clear all
set more off

* Set the working directory to the folder that holds the data before running,
* for example:  cd "path/to/data"

* Load data
import delimited "Earnings_and_Height.csv", clear

* Quick data inspection
describe
summarize

********************************************************************************
* DATA OVERVIEW
********************************************************************************
* Dataset: Random sample of U.S. workers (n = 17,870)
* Key variables:
*   - earnings: Annual earnings ($)
*   - height: Height in inches
*   - sex: 0 = Female, 1 = Male
*   - occupation: Occupation category (1-15)

* Check sample sizes by gender
tab sex

********************************************************************************
* KEY FORMULAS FOR HYPOTHESIS TESTING
********************************************************************************
* t-statistic: t = (β̂₁ - β₁,₀) / SE(β̂₁)
* Under H₀: β₁ = 0, this becomes t = β̂₁ / SE(β̂₁)
*
* Critical values (two-sided):
*   - 10% level: |t| > 1.645
*   - 5% level:  |t| > 1.96
*   - 1% level:  |t| > 2.576
*
* 95% CI: β̂₁ ± 1.96 × SE(β̂₁)
* 90% CI: β̂₁ ± 1.645 × SE(β̂₁)

********************************************************************************
* (a) REGRESSION: ALL WORKERS
********************************************************************************

di "=============================================="
di "(a) REGRESSION: ALL WORKERS"
di "=============================================="

regress earnings height

* Store results for display
local b1 = _b[height]
local se1 = _se[height]
local t1 = _b[height] / _se[height]
local p1 = 2 * ttail(e(df_r), abs(`t1'))
local r2 = e(r2)
local n = e(N)

* Calculate 95% CI manually
local ci95_lo = `b1' - 1.96 * `se1'
local ci95_hi = `b1' + 1.96 * `se1'

di ""
di "RESULTS:"
di "  Slope (β̂₁):        " %10.2f `b1'
di "  Std. Error:        " %10.2f `se1'
di "  t-statistic:       " %10.2f `t1'
di "  p-value:           " %10.4f `p1'
di ""
di "(i) Is the slope statistically significant?"
di "    |t| = " %5.2f abs(`t1') " > 1.96 → Significant at 5% level"
di "    |t| = " %5.2f abs(`t1') " > 2.576 → Significant at 1% level"
di "    p-value ≈ 0 < 0.01 → REJECT H₀: β₁ = 0"
di ""
di "(ii) 95% Confidence Interval:"
di "     " %8.2f `ci95_lo' " to " %8.2f `ci95_hi'
di "     (Does not include 0 → Significant)"

********************************************************************************
* (b) REGRESSION: WOMEN ONLY (sex == 0)
********************************************************************************

di ""
di "=============================================="
di "(b) REGRESSION: WOMEN ONLY"
di "=============================================="

regress earnings height if sex == 0

* Store results
local b1_f = _b[height]
local se1_f = _se[height]
local t1_f = _b[height] / _se[height]
local p1_f = 2 * ttail(e(df_r), abs(`t1_f'))
local n_f = e(N)

* Calculate 95% CI
local ci95_lo_f = `b1_f' - 1.96 * `se1_f'
local ci95_hi_f = `b1_f' + 1.96 * `se1_f'

di ""
di "RESULTS:"
di "  Sample size:       " %10.0f `n_f'
di "  Slope (β̂₁):        " %10.2f `b1_f'
di "  Std. Error:        " %10.2f `se1_f'
di "  t-statistic:       " %10.2f `t1_f'
di "  p-value:           " %10.4f `p1_f'
di ""
di "(i) Is the slope statistically significant?"
di "    |t| = " %5.2f abs(`t1_f') " > 1.96 → Significant at 5% level"
di "    p-value < 0.01 → REJECT H₀: β₁ = 0"
di ""
di "(ii) 95% Confidence Interval:"
di "     " %8.2f `ci95_lo_f' " to " %8.2f `ci95_hi_f'

********************************************************************************
* (c) REGRESSION: MEN ONLY (sex == 1)
********************************************************************************

di ""
di "=============================================="
di "(c) REGRESSION: MEN ONLY"
di "=============================================="

regress earnings height if sex == 1

* Store results
local b1_m = _b[height]
local se1_m = _se[height]
local t1_m = _b[height] / _se[height]
local p1_m = 2 * ttail(e(df_r), abs(`t1_m'))
local n_m = e(N)

* Calculate 95% CI
local ci95_lo_m = `b1_m' - 1.96 * `se1_m'
local ci95_hi_m = `b1_m' + 1.96 * `se1_m'

di ""
di "RESULTS:"
di "  Sample size:       " %10.0f `n_m'
di "  Slope (β̂₁):        " %10.2f `b1_m'
di "  Std. Error:        " %10.2f `se1_m'
di "  t-statistic:       " %10.2f `t1_m'
di "  p-value:           " %10.6f `p1_m'
di ""
di "(i) Is the slope statistically significant?"
di "    |t| = " %5.2f abs(`t1_m') " > 2.576 → Significant at 1% level"
di "    p-value ≈ 0 → REJECT H₀: β₁ = 0"
di ""
di "(ii) 95% Confidence Interval:"
di "     " %8.2f `ci95_lo_m' " to " %8.2f `ci95_hi_m'

********************************************************************************
* COMPARISON TABLE: All vs Women vs Men
********************************************************************************

di ""
di "=============================================="
di "COMPARISON: Regression Results by Gender"
di "=============================================="
di ""
di "                    All         Women        Men"
di "------------------------------------------------------"
di "Sample size      " %8.0f `n' "     " %8.0f `n_f' "    " %8.0f `n_m'
di "Slope ($/inch)   " %8.2f `b1' "     " %8.2f `b1_f' "    " %8.2f `b1_m'
di "Std. Error       " %8.2f `se1' "     " %8.2f `se1_f' "    " %8.2f `se1_m'
di "t-statistic      " %8.2f `t1' "     " %8.2f `t1_f' "    " %8.2f `t1_m'
di "------------------------------------------------------"
di "95% CI Lower     " %8.2f `ci95_lo' "     " %8.2f `ci95_lo_f' "    " %8.2f `ci95_lo_m'
di "95% CI Upper     " %8.2f `ci95_hi' "     " %8.2f `ci95_hi_f' "    " %8.2f `ci95_hi_m'
di "------------------------------------------------------"
di "Significant?     " "    YES          YES         YES"
di "=============================================="

********************************************************************************
* (d) TEST: Is the Effect the Same for Men and Women?
********************************************************************************

di ""
di "=============================================="
di "(d) TEST: Is the Height Effect the Same for Men and Women?"
di "=============================================="
di ""
di "Method: Interaction Model"
di "Model: Earnings = β₀ + β₁×Height + β₂×Male + β₃×(Height×Male) + u"
di ""
di "H₀: β₃ = 0 (effect of height is the same for men and women)"
di "H₁: β₃ ≠ 0 (effect of height differs by gender)"

* Create interaction variables
gen male = (sex == 1)
gen height_male = height * male

* Run interaction regression
regress earnings height male height_male

* Store coefficient on interaction term
local b3 = _b[height_male]
local se3 = _se[height_male]
local t3 = _b[height_male] / _se[height_male]
local p3 = 2 * ttail(e(df_r), abs(`t3'))

* Calculate 95% CI for interaction
local ci95_lo_int = `b3' - 1.96 * `se3'
local ci95_hi_int = `b3' + 1.96 * `se3'

di ""
di "RESULTS:"
di "  Coefficient on Height×Male (β̂₃): " %10.2f `b3'
di "  Std. Error:                       " %10.2f `se3'
di "  t-statistic:                      " %10.2f `t3'
di "  p-value:                          " %10.4f `p3'
di ""
di "  95% CI for β₃: (" %8.2f `ci95_lo_int' ", " %8.2f `ci95_hi_int' ")"
di ""
di "INTERPRETATION:"
di "  β₁ = " %6.2f _b[height] " (effect for women: base group)"
di "  β₁ + β₃ = " %6.2f _b[height] + `b3' " (effect for men)"
di ""
if `p3' < 0.05 {
    di "CONCLUSION: REJECT H₀ at 5% level"
    di "  The effect of height on earnings is significantly DIFFERENT"
    di "  for men and women. Men gain $" %6.2f `b3' " more per inch than women."
}
else {
    di "CONCLUSION: FAIL TO REJECT H₀ at 5% level"
    di "  Cannot conclude the effect differs by gender."
}

********************************************************************************
* (e) OCCUPATIONS WHERE STRENGTH IS UNLIKELY IMPORTANT
********************************************************************************

di ""
di "=============================================="
di "(e) Does Height Effect Disappear in Non-Physical Occupations?"
di "=============================================="

* First, check occupation variable
tab occupation

* Run regressions by occupation category
di ""
di "Regression results by occupation:"
di ""
di "Occup.  Slope      SE       t-stat    p-value   n"
di "------------------------------------------------------"

forvalues occ = 1/15 {
    quietly reg earnings height if occupation == `occ'
    local b_occ = _b[height]
    local se_occ = _se[height]
    local t_occ = _b[height] / _se[height]
    local p_occ = 2 * ttail(e(df_r), abs(`t_occ'))
    local n_occ = e(N)
    di %4.0f `occ' "    " %8.2f `b_occ' "  " %8.2f `se_occ' "  " %7.2f `t_occ' "  " %7.4f `p_occ' "  " %6.0f `n_occ'
}

* Group occupations
* Typically: 1-5 might be more physical/manual, 6+ might be professional/office
di ""
di "------------------------------------------------------"
di "GROUPED ANALYSIS:"
di ""

* Professional/Office occupations (where strength less likely to matter)
* This depends on the actual occupation coding - adjust as needed
di "Professional/Service occupations (e.g., occupation >= 6):"
quietly reg earnings height if occupation >= 6
di "  Slope: $" %6.2f _b[height] "/inch, t = " %5.2f _b[height]/_se[height] ", n = " e(N)

di ""
di "Manual/Physical occupations (e.g., occupation <= 5):"
quietly reg earnings height if occupation <= 5
di "  Slope: $" %6.2f _b[height] "/inch, t = " %5.2f _b[height]/_se[height] ", n = " e(N)

* Run full regression for non-physical occupations for detailed output
di ""
di "Detailed regression for occupations where strength unlikely important:"
reg earnings height if occupation >= 6

********************************************************************************
* END OF E5.1
********************************************************************************