"""教学模拟，只有显式运行才写指定成果目录，不读取研究数据、不联网或安装软件。"""
from pathlib import Path
import argparse,csv,json,platform
import numpy as np
SEED=42
def generate_sample(n: int = 30, seed: int = SEED, *, heteroskedastic: bool = False):
    """返回 x/y/真实误差；两种生成方式均保持独立观测、零条件均值和已知斜率 2。"""
    if n < 3:
        raise ValueError("至少需要3个观察，才能估计截距/斜率并计算标准误。")
    rng = np.random.default_rng(seed)
    x = rng.uniform(-2, 2, n) if heteroskedastic else rng.normal(0, 1, n)
    sigma = 0.3 + 1.2 * np.abs(x) if heteroskedastic else np.ones(n)
    error = sigma * rng.normal(0, 1, n)
    return x, 1 + 2 * x + error, error

def fit_ols(x, y):
    """含截距的单变量OLS，返回系数、普通/HC1标准误、拟合值和残差。"""
    x, y = np.asarray(x, dtype=float), np.asarray(y, dtype=float)
    if x.ndim != 1 or y.shape != x.shape or len(x) <= 2:
        raise ValueError("x/y须为等长一维数组，且观察数大于2。")
    if not np.isfinite(x).all() or not np.isfinite(y).all():
        raise ValueError("示例不接受缺失或无穷值；先明确样本处理规则。")
    design = np.column_stack([np.ones(len(x)), x])
    beta, _, rank, _ = np.linalg.lstsq(design, y, rcond=None)
    if rank != 2:
        raise ValueError("解释变量必须有变化，否则无法识别斜率。")
    fitted = design @ beta
    residual = y - fitted
    n, k = design.shape
    bread = np.linalg.inv(design.T @ design)
    classical_cov = (residual @ residual) / (n - k) * bread
    # HC1只调整协方差估计，不重估系数；n/(n-k)是有限样本修正，k包括截距。
    scores = design * residual[:, None]
    hc1_cov = n / (n - k) * bread @ (scores.T @ scores) @ bread
    return {
        "beta": beta,
        "classical_se": np.sqrt(np.diag(classical_cov)),
        "hc1_se": np.sqrt(np.diag(hc1_cov)),
        "fitted": fitted,
        "residual": residual,
    }

def run(output,seed=42):
    from scipy.stats import t
    import scipy,matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    output.mkdir(parents=True,exist_ok=True)
    x,y,error=generate_sample(500,seed,heteroskedastic=True);fit=fit_ols(x,y);df=len(x)-2;crit=float(t.ppf(.975,df))
    with (output/'simulated-data.csv').open('w',newline='') as f:
        w=csv.writer(f,lineterminator='\n');w.writerow(['x','y','error','residual']);w.writerows(zip(x,y,error,fit['residual']))
    rows=[]
    for name,key in [('conventional','classical_se'),('HC1','hc1_se')]:
        for j,term in enumerate(['intercept','x']):
            b=float(fit['beta'][j]);se=float(fit[key][j]);p=float(2*t.sf(abs(b/se),df))
            rows.append({'variance':name,'term':term,'n':len(x),'coefficient':b,'standard_error':se,'ci_lower':b-crit*se,'ci_upper':b+crit*se,'p_value':p})
    with (output/'estimates.csv').open('w',newline='') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]),lineterminator='\n');w.writeheader();w.writerows(rows)
    fig,ax=plt.subplots(1,2,figsize=(12,5),layout='constrained')
    ax[0].scatter(x,fit['residual'],s=12,alpha=.6,color='#615b49');ax[0].axhline(0,color='#333',linestyle='--');ax[0].set(title='Same simulated sample (n=500)',xlabel='x',ylabel='OLS residual')
    pos=np.arange(2)
    ax[1].bar(pos-.18,fit['classical_se'],.36,label='Conventional',color='#D4C3A1')
    ax[1].bar(pos+.18,fit['hc1_se'],.36,label='HC1',color='#615b49')
    ax[1].set(title='Same coefficients; different SE estimates',xticks=pos,xticklabels=['Intercept','Slope'],ylabel='Standard error');ax[1].legend(frameon=False)
    for a in ax:a.grid(alpha=.15)
    for ext in ['svg','png']:fig.savefig(output/f'robust-standard-errors.{ext}',dpi=160)
    plt.close(fig)
    svg=output/'robust-standard-errors.svg';svg.write_text('\n'.join(l.rstrip() for l in svg.read_text().splitlines())+'\n')
    result={'simulated':True,'seed':seed,'n':500,'model':'y=1+2x+(0.3+1.2|x|)z; x~Uniform(-2,2), z~N(0,1), independent observations','df':df,'rows':rows,'runtime':{'python':platform.python_version(),'numpy':np.__version__,'scipy':scipy.__version__,'matplotlib':matplotlib.__version__}}
    (output/'result.json').write_text(json.dumps(result,indent=2)+'\n');return result
if __name__=='__main__':
    p=argparse.ArgumentParser();p.add_argument('--seed',type=int,default=42);p.add_argument('--output-dir',type=Path,default=Path(__file__).resolve().parent/'results');a=p.parse_args()
    print(json.dumps(run(a.output_dir,a.seed),indent=2));print('ROBUST_SE_COMPLETED')