clear
set seed 123456
local num_of_replications = 1000
* treatment effects (absolute deviations)
local narrative_effect             = 0.00
local nondiagnostic_effect         = 0.00
local narr_and_nondiagnostic_effect = 2.15
* baseline outcome distribution (truncated normal on [0,20])
local mu = 8.9
local sd = 5
local L  = 0
local U  = 20
local subjects = 184
local T = 15


set matsize `num_of_replications'
mat estimates = J(`num_of_replications',1,.)


* ---- within-subject correlation ----
* Random-intercept ICC target: ICC = su^2 / (su^2 + se^2)
local ICC   = 0.50          // set to 0..1 (0.40 gives strong same-person similarity)
local totSD = 3             // SD of the combined {random intercept + idiosyncratic} part
local su = sqrt(`ICC' * `totSD'^2)
local se = sqrt((1-`ICC') * `totSD'^2)

quietly forvalues j = 1/`num_of_replications' {
    clear
    set obs `subjects'
    gen id = _n

    * subject-level treatment (1/4 per cell)
    gen str14 treatment = "Control_no"
    replace treatment="Control_yes"    if id > `subjects'/4   & id <= `subjects'/2
    replace treatment="Narrative_no"   if id > `subjects'/2   & id <= `subjects'*3/4
    replace treatment="Narrative_yes"  if id > `subjects'*3/4

    * panel: 15 rows per subject
    expand `T'
    bysort id: gen t = _n

    gen byte narrative     = inlist(treatment,"Narrative_no","Narrative_yes")
    gen byte nondiagnostic = inlist(treatment,"Control_yes","Narrative_yes")

    * ---- baseline price per row (truncated normal on [0,20]) ----
    scalar pL = normal((`L' - `mu')/`sd')
    scalar pU = normal((`U' - `mu')/`sd')
    gen double u_draw = runiform(pL,pU)
    gen double price  = `mu' + `sd'*invnormal(u_draw)
    drop u_draw

    * ---- subject random intercept (creates ICC) ----
    bysort id: gen double u_i = .
    bysort id: replace u_i = rnormal(0, `su') if _n==1
    bysort id: replace u_i = u_i[1]

    * ---- idiosyncratic error: i.i.d.
	gen double e = rnormal(0, `se')

    * ---- apply treatment shifts ----
    replace price = price + `narrative_effect'                  if treatment=="Narrative_no"
    replace price = price + `nondiagnostic_effect'              if treatment=="Control_yes"
    replace price = price + `narr_and_nondiagnostic_effect'     if treatment=="Narrative_yes"

    * add correlated person-level and serial components
    replace price = price + u_i + e

    * regress with clustering at subject level
	reg price i.nondiagnostic##i.narrative, vce(cluster id)
    lincom 1.nondiagnostic#1.narrative
    scalar pvalue = r(p)
	matrix estimates[`j',1] = pvalue 
}

svmat estimates, names(pvalues)
gen significant = (pvalues < 0.05) if pvalues != .
ci means significant