"""Verify Stata outputs against Excel and regenerate the two offline reports.
Requires: numpy, scipy, openpyxl. Run after week7_analysis.do.
"""
from pathlib import Path
import base64
import csv
import hashlib
import html
import json
import re
import xml.etree.ElementTree as ET
import numpy as np
from scipy import stats
from openpyxl import load_workbook

ROOT = Path(__file__).resolve().parent

def read_csv(name):
    with (ROOT / name).open(encoding='utf-8-sig', newline='') as f:
        return list(csv.DictReader(f))

def numeric(row):
    return {k: float(v) for k, v in row.items()}

model = numeric(read_csv('model_statistics.csv')[0])
coeff = {r['term']: {k: float(v) for k, v in r.items() if k != 'term'} for r in read_csv('regression_coefficients.csv')}
robust = {r['term']: {k: float(v) for k, v in r.items() if k != 'term'} for r in read_csv('robust_coefficients.csv')}
groups = [numeric(r) for r in read_csv('income_group_summary.csv')]
audit = [numeric(r) for r in read_csv('source_totals_audit.csv')]
fitted = [numeric(r) for r in read_csv('households_fitted.csv')]

# Independent numerical validation: derive everything directly from source cells.
workbook = load_workbook(ROOT / 'Eg2_1_1/Eg2_1_1.xlsx', data_only=True)
sheet = workbook['Sheet1']
source = []
for col in range(3, 13):
    for row in range(4, 17):
        value = sheet.cell(row, col).value
        if value is not None:
            source.append((row, float(sheet.cell(3, col).value), float(value)))
source_array = np.array(source)
assert len(source) == 98
assert all((r['excel_row'], r['income'], r['consume']) == item for r, item in zip(fitted, source))
x, y = source_array[:, 1], source_array[:, 2]
X = np.column_stack([np.ones(len(x)), x])
beta = np.linalg.lstsq(X, y, rcond=None)[0]
yhat = X @ beta
resid = y - yhat
n, k = X.shape
inv = np.linalg.inv(X.T @ X)
rss = resid @ resid
se = np.sqrt(np.diag((rss / (n-k)) * inv))
hc1 = inv @ (X.T @ (X * (resid ** 2)[:, None])) @ inv * n / (n-k)
se_hc1 = np.sqrt(np.diag(hc1))
np.testing.assert_allclose([model['intercept'], model['slope']], beta, rtol=1e-10, atol=1e-9)
np.testing.assert_allclose([r['fitted'] for r in fitted], yhat, rtol=1e-10, atol=1e-8)
np.testing.assert_allclose([r['residual'] for r in fitted], resid, rtol=1e-10, atol=1e-8)
for i, term in enumerate(['_cons', 'income']):
    np.testing.assert_allclose(coeff[term]['se'], se[i], rtol=1e-10)
    np.testing.assert_allclose(robust[term]['se'], se_hc1[i], rtol=1e-10)
    for output, standard_error in [(coeff[term], se[i]), (robust[term], se_hc1[i])]:
        p = 2 * stats.t.sf(abs(beta[i] / standard_error), n-k)
        np.testing.assert_allclose(output['p'], p, rtol=1e-9, atol=0)
np.testing.assert_allclose(model['r_squared'], 1-rss/np.sum((y-y.mean())**2), rtol=1e-12)
np.testing.assert_allclose(model['root_mse'], np.sqrt(rss/(n-k)), rtol=1e-12)
for g, a in zip(groups, audit):
    yy = y[x == g['income']]
    assert g['n'] == len(yy) == a['n']
    assert a['detail_total'] == yy.sum()
    np.testing.assert_allclose(g['mean_consume'], yy.mean(), rtol=1e-12)
    assert a['difference'] == a['detail_total'] - a['reported_total']
log = (ROOT / 'week7_stata.log').read_text(encoding='utf-8', errors='replace')
assert 'WEEK7_ANALYSIS_COMPLETED' in log
assert not re.search(r'^r\(\d+\);', log, flags=re.M)
assert (ROOT / 'consumption_fit.png').read_bytes().startswith(b'\x89PNG')
assert (ROOT / 'consumption_fit.pdf').read_bytes().startswith(b'%PDF')
ET.parse(ROOT / 'consumption_fit.svg')
assert '不同可支配收入水平' in (ROOT / 'week7_analysis.do').read_text(encoding='utf-8')
verification = {
    'status': 'passed', 'date': '2026-10-23', 'observations': n,
    'source_sha256': hashlib.sha256((ROOT/'Eg2_1_1/Eg2_1_1.xlsx').read_bytes()).hexdigest(),
    'independent_ols_intercept': float(beta[0]), 'independent_ols_slope': float(beta[1]),
    'max_abs_fitted_difference': float(np.max(np.abs(np.array([r['fitted'] for r in fitted])-yhat))),
    'checks': ['all 98 Excel source records match Stata data, including worksheet row',
               'OLS coefficients, fitted values, residuals, R-squared, RMSE',
               'conventional and HC1 standard errors and two-sided t-test p-values',
               'all income-group counts, means, and source-total differences',
               'Stata log completion and no r() errors', 'PNG, SVG, and PDF file structure']
}
(ROOT / 'verification.json').write_text(json.dumps(verification, ensure_ascii=False, indent=2), encoding='utf-8')

md, body = [], []
sections = []

def inline(text):
    text = html.escape(text)
    text = re.sub(r'`([^`]+)`', r'<code>\1</code>', text)
    text = re.sub(r'\*\*([^*]+)\*\*', r'<strong>\1</strong>', text)
    return text

def heading(title, anchor):
    md.append('## ' + title)
    body.append(f'<h2 id="{anchor}">{html.escape(title)}</h2>')
    sections.append((title, anchor))

def paragraph(text, cls=''):
    md.append(text)
    body.append(f'<p class="{cls}">{inline(text)}</p>')

def equation(text):
    md.append('> ' + text)
    body.append(f'<div class="equation">{html.escape(text)}</div>')

def table(headers, rows):
    md.append('| ' + ' | '.join(headers) + ' |\n| ' + ' | '.join(['---']*len(headers)) + ' |\n' + '\n'.join('| ' + ' | '.join(str(c).replace('|', r'\|') for c in r) + ' |' for r in rows))
    body.append('<div class="table-wrap"><table><thead><tr>' + ''.join('<th scope="col">'+inline(h)+'</th>' for h in headers) + '</tr></thead><tbody>' + ''.join('<tr>'+''.join('<td>'+inline(str(c))+'</td>' for c in row)+'</tr>' for row in rows) + '</tbody></table></div>')

def code(text, lang='stata'):
    md.append('```'+lang+'\n'+text+'\n```')
    body.append('<pre><code>'+html.escape(text)+'</code></pre>')

def image(path, alt):
    md.append(f'![{alt}]({path})')
    encoded = base64.b64encode((ROOT/path).read_bytes()).decode()
    body.append(f'<figure><img src="data:image/png;base64,{encoded}" alt="{html.escape(alt)}"><figcaption>{html.escape(alt)}</figcaption></figure>')

def asset_link(path, label):
    return f'<a href="{html.escape(path, quote=True)}" download>{html.escape(label)}</a>'

md.append('# Week 7：随机扰动项、残差与家庭消费回归')
paragraph('研究日期：2026 年 10 月 23 日。数据：Eg2_1_1/Eg2_1_1.xlsx（Sheet1）。实际执行软件：Stata/SE 15.1，路径 D:/Stata15/StataSE-64.exe；脚本使用 version 15.0 语法。')
paragraph('本研究使用原表 98 条家庭消费明细，估计每月可支配收入与每月消费支出的线性关系，并参照给定图片完成 Stata 作图。所有金额单位均为元。')
equation(f'Ŷ = {model["intercept"]:.4f} + {model["slope"]:.6f}X；n = 98；R² = {model["r_squared"]:.6f}')

heading('1. 随机扰动项与残差项', 'concepts')
paragraph('**随机扰动项** uᵢ 是实际观测值与真实总体条件均值之间的偏差。写作 Yᵢ = m(Xᵢ) + uᵢ，其中 m(Xᵢ) = E(Yᵢ | Xᵢ)。它反映给定收入后，家庭偏好、人口规模、其他未纳入因素及随机变化造成的消费差异。真实总体条件均值通常未知，因此 uᵢ 通常不可直接观测。')
equation('uᵢ = Yᵢ − E(Yᵢ | Xᵢ)；若 m(Xᵢ) = β₀ + β₁Xᵢ，则 uᵢ = Yᵢ − β₀ − β₁Xᵢ。')
paragraph('**残差项** eᵢ 是实际观测值与拟合值之间的偏差。它在估计回归系数后可以计算，是对相应随机扰动的估计。')
equation('eᵢ = Yᵢ − Ŷᵢ = Yᵢ − β̂₀ − β̂₁Xᵢ。')
table(['比较维度', '随机扰动项 uᵢ', '残差项 eᵢ'], [
    ['参照对象', '真实总体回归函数 E(Y | X)', '利用数据估计的回归函数 Ŷ'],
    ['参数', '真实且通常未知的 β₀、β₁', '可由数据算出的 β̂₀、β̂₁'],
    ['能否直接计算', '通常不能', '拟合后可以'],
    ['零均值性质', '以真实条件均值定义时 E(u | X) = 0；设定线性模型时需相应条件成立', '含截距 OLS 的样本残差满足 Σeᵢ = 0（数值计算有舍入误差）'],
    ['含义', '实际值相对真实条件均值的偏离', '实际值相对估计条件均值的偏离']
])
equation('eᵢ = uᵢ − (β̂₀ − β₀) − (β̂₁ − β₁)Xᵢ。')
paragraph('因此，残差与扰动项并非完全相同。含截距 OLS 还满足 ΣXᵢeᵢ = 0；这些是拟合所得的样本代数性质，不能据此证明总体扰动与解释变量无关。若线性条件均值设定不成立，拟合残差还会包含函数形式误设的影响。')
paragraph('**关于“总体回归曲线”**：总体回归函数是 E(Y | X)，本题将它设定为 β₀ + β₁X。工作簿没有说明所列家庭是否覆盖某个明确总体的全部成员，因此本文将拟合直线称为“总体回归函数的估计（样本回归线）”。各收入组的实际平均消费是组内经验均值，也不直接当作已知的真实总体条件均值。', 'note')

heading('2. 数据整理与原表核对', 'data')
paragraph('工作表采用横向分组布局：C3:L3 为 10 个收入水平（800、1100、…、3500），C4:L16 为家庭消费明细，C17:L17 为原表合计。标题文字在源工作簿中已有乱码，数字单元格可正常读取；变量含义依据给定示例图和表格结构确认。')
paragraph('整理过程：读取 C3:L17；核对收入表头；分别计算每列明细和并与末行合计比较；移除表头及合计行；将消费列转换为长表；删除空白明细；保留收入、消费和原 Excel 行号。空白不视为 0；相同金额的重复记录保留，因为可能对应不同家庭。每条明细等权，共得到 98 条记录。')
paragraph('**发现 4 组明细与原表合计不一致。** 差额按“明细之和 − 原表合计”定义。以下差异不足以判断究竟是明细录入错误、合计错误还是遗漏，因此不补造记录，也不按合计反推消费值；回归仅依据现有明细。', 'warning')
table(['收入 X', '记录数', '明细之和', '原表合计', '差额'], [[f'{a["income"]:.0f}', f'{a["n"]:.0f}', f'{a["detail_total"]:.0f}', f'{a["reported_total"]:.0f}', f'{a["difference"]:+.0f}'] for a in audit])
paragraph('特别是收入 2300 元组，明细之和为 21868 元，而原表合计为 23870 元，差额为 −2002 元。应先向数据提供者核实，再考虑更正；本次结果对应当前文件的可见明细。')

heading('3. 回归设定、过程与估计结果', 'results')
paragraph('设 X 为家庭每月可支配收入，Y 为每月消费支出。假定总体条件均值为线性形式：E(Y | X) = β₀ + β₁X；使用普通最小二乘法（OLS）选择截距和斜率，使 98 条明细的残差平方和最小。[1]')
equation('min Σ(Yᵢ − β̂₀ − β̂₁Xᵢ)²；β̂₁ = Σ(Xᵢ − X̄)(Yᵢ − Ȳ) / Σ(Xᵢ − X̄)²；β̂₀ = Ȳ − β̂₁X̄。')
paragraph(f'样本平均收入 X̄ = {x.mean():.4f} 元，平均消费 Ȳ = {y.mean():.4f} 元。Stata 命令为 `regress consume income`，并使用 `predict ..., xb` 与 `predict ..., residuals` 保存拟合值及残差。[1][2]')
equation(f'Ŷ = {model["intercept"]:.6f} + {model["slope"]:.9f}X')
table(['变量', '系数', '常规标准误', 't 值', 'p 值', '常规 95% 置信区间'], [
    ['收入 X', f'{coeff["income"]["estimate"]:.6f}', f'{coeff["income"]["se"]:.6f}', f'{coeff["income"]["t"]:.4f}', f'{coeff["income"]["p"]:.3e}', f'[{coeff["income"]["ci_low"]:.6f}, {coeff["income"]["ci_high"]:.6f}]'],
    ['截距', f'{coeff["_cons"]["estimate"]:.6f}', f'{coeff["_cons"]["se"]:.6f}', f'{coeff["_cons"]["t"]:.4f}', f'{coeff["_cons"]["p"]:.6f}', f'[{coeff["_cons"]["ci_low"]:.6f}, {coeff["_cons"]["ci_high"]:.6f}]']
])
table(['指标', '结果'], [
    ['观测数 / 残差自由度', '98 / 96'],
    ['R² / 调整 R²', f'{model["r_squared"]:.6f} / {1-(1-model["r_squared"])*(n-1)/(n-k):.6f}'],
    ['Root MSE（元）', f'{model["root_mse"]:.6f}'],
    ['残差平方和（元²）', f'{model["rss"]:.6f}'],
    ['常规 F(1,96) / p 值', f'{model["F"]:.6f} / {model["p_F"]:.3e}']
])
paragraph('**经济含义**：在这组数据的线性拟合中，收入每高 1 元，拟合的平均消费约高 0.7288 元；每高 100 元，约高 72.88 元。可将斜率描述为估计的边际消费倾向，但单变量描述性回归不能独自识别收入对消费的因果作用。')
paragraph('R² = 0.9307 表示收入的线性拟合解释了本数据中约 93.07% 的消费总离差平方和。它不是预测准确率，也不能证明不存在遗漏变量。截距对应收入为 0 时的拟合消费，但 0 超出观测收入范围 800—3500 元，因而不作实质经济解释。')
paragraph(f'**补充稳健性核对**：`estat hettest` 的 Breusch–Pagan / Cook–Weisberg 检验为 χ²(1) = {model["bp_chi2"]:.4f}，p = {model["bp_p"]:.6f}，在其检验条件下拒绝同方差原假设。另运行 `regress consume income, vce(robust)` 计算 HC1 稳健标准误；回归系数和拟合直线保持一致。[1][2]')
table(['变量', 'HC1 稳健标准误', '稳健 t 值', '稳健 p 值', '稳健 95% 置信区间'], [
    ['收入 X', f'{robust["income"]["se"]:.6f}', f'{robust["income"]["t"]:.4f}', f'{robust["income"]["p"]:.3e}', f'[{robust["income"]["ci_low"]:.6f}, {robust["income"]["ci_high"]:.6f}]'],
    ['截距', f'{robust["_cons"]["se"]:.6f}', f'{robust["_cons"]["t"]:.4f}', f'{robust["_cons"]["p"]:.6f}', f'[{robust["_cons"]["ci_low"]:.6f}, {robust["_cons"]["ci_high"]:.6f}]']
])
paragraph('采用稳健标准误后，收入系数仍显著为正。稳健推断仍依赖相应抽样、独立性及大样本近似条件；数据没有提供抽样设计信息，且存在上述合计差异，所以结果宜作本组数据的教学性、描述性分析。')

heading('4. 仿照示例图作图与残差解释', 'figure')
paragraph('Stata 使用 `twoway (scatter consume income) (lfit consume income)` 作图：深蓝色圆点表示家庭消费明细，红线表示 OLS 拟合的条件均值。横轴为每月可支配收入，纵轴为每月消费支出；淡蓝背景、白色绘图区、中文标题及底部图例均参照给定图片。[3]')
image('consumption_fit.png', '图 1｜Stata 生成的家庭消费散点与线性拟合图')
paragraph('每个固定收入水平上的竖直点列展示该收入组的经验消费分布。若收入、消费均相同，点会重合，所以图上可分辨的点数可能小于 98；这些观测在回归中全部保留。红线估计条件均值，不能代表每户家庭的实际消费。')
paragraph('按收入组得到的明细均值与拟合均值如下。回归直接使用 98 条家庭明细，未对 10 个组均值作不加权回归；各组家庭数量不同，不能随意替换分析单位。')
table(['收入 X', '家庭数', '组内平均消费', '拟合均值 Ŷ', '平均残差'], [[f'{g["income"]:.0f}', f'{g["n"]:.0f}', f'{g["mean_consume"]:.3f}', f'{g["fitted_mean"]:.3f}', f'{g["mean_residual"]:+.3f}'] for g in groups])
paragraph(f'**残差例子**：Excel 单元格 C4 对应收入 800 元、消费 561 元。其拟合消费为 {yhat[0]:.4f} 元，所以 e = 561 − {yhat[0]:.4f} = {resid[0]:.4f} 元。负号表示这户家庭的实际消费低于拟合条件均值；这不是已知的真实随机扰动 u。')
paragraph(f'数值核对：Σeᵢ = {model["residual_sum"]:.3e}，ΣXᵢeᵢ = {model["income_residual_sum"]:.3e}。两者在浮点精度下接近 0，与含截距 OLS 的正规方程一致；Root MSE = √[Σeᵢ²/(98−2)] = {model["root_mse"]:.4f} 元。')
body.append(f'''<div class="prediction"><label for="income-slider">在观测范围内查看拟合均值</label><div class="prediction-values"><span>每月收入 <output id="income-value">2000</output> 元</span><strong>拟合平均消费 <output id="pred-value">{model['intercept']+model['slope']*2000:.2f}</output> 元</strong></div><input id="income-slider" type="range" min="800" max="3500" step="10" value="2000"><p>中间收入值是线性插值；结果表示拟合条件均值，并非对某户实际消费的保证。</p></div>''')
body.append('<details><summary>查看题目提供的原始参考图</summary>')
# The original image remains available in Markdown as a regular reference link.
md.append('[原始参考图](Eg2_1_1/Eg2_1_1.png)')
original_img = base64.b64encode((ROOT/'Eg2_1_1/Eg2_1_1.png').read_bytes()).decode()
body.append(f'<img src="data:image/png;base64,{original_img}" alt="题目提供的原始参考图"></details>')

heading('5. 复现方法与 Stata 指令', 'reproduce')
paragraph('完整脚本 `week7_analysis.do` 从原始 Excel 开始，依次完成合计审计、长表转换、OLS 回归、拟合值与残差生成、图形导出、稳健标准误及异方差检验。脚本无需安装 Stata 扩展包。在 Stata 命令窗口运行：')
code('do "C:/Users/cattom/Documents/econometrics/Week7/week7_analysis.do"')
paragraph('也可以在 Windows PowerShell 中后台执行。以下命令会等待 Stata 完成；重跑会更新同名分析产物：')
code('Start-Process -FilePath "D:\\Stata15\\StataSE-64.exe" `\n  -ArgumentList \'/e do "C:\\Users\\cattom\\Documents\\econometrics\\Week7\\week7_analysis.do"\' `\n  -WorkingDirectory "C:\\Users\\cattom\\Documents\\econometrics\\Week7" `\n  -WindowStyle Hidden -Wait', 'powershell')
paragraph('下列为主要估计与绘图指令摘录；前置数据整理和完整图形样式以完整 do 文件为准：')
code('use "C:/Users/cattom/Documents/econometrics/Week7/households_clean.dta", clear\nregress consume income\npredict double fitted, xb\npredict double residual, residuals\ntwoway (scatter consume income) (lfit consume income), ///\n    title("不同可支配收入水平组家庭消费支出的条件分布图") ///\n    xtitle("每月可支配收入（元）") ytitle("每月消费支出（元）")\nregress consume income, vce(robust)')
paragraph('Stata 运行后，在工作目录执行以下命令可独立核对数值并重建两份报告。此辅助脚本需要 Python、NumPy、SciPy 和 openpyxl；它读取 Stata 输出，不替代 Stata 回归或作图。')
code('python build_week7_report.py', 'powershell')
paragraph('**已完成验证**：独立从 Excel 明细构建设计矩阵，用 NumPy 复算 OLS，并核对 98 条记录、系数、拟合值、残差、R²、Root MSE、常规及 HC1 标准误和双侧 t 检验 p 值；核对分组汇总和合计差异；检查 Stata 日志成功结束及图形文件结构。验证结果保存于 `verification.json`。')

heading('6. 交付文件索引', 'files')
assets = [
    ('Week7_回归分析报告.md', 'Markdown 文档，包含研究过程、结果及图片相对引用'),
    ('Week7_回归分析报告.html', 'HTML5 页面，图形与样式内嵌，可离线阅读'),
    ('week7_analysis.do', '完整、可复现的 Stata 分析脚本'),
    ('week7_stata.log', 'Stata 纯文本分析日志，含回归、分组及检验输出'),
    ('consumption_fit.png', '仿照示例生成的 1800 像素宽图形'),
    ('consumption_fit.svg', '可缩放矢量图'),
    ('consumption_fit.pdf', '适合打印的图形'),
    ('consumption_fit.gph', 'Stata 原生图形，可用 graph use 打开'),
    ('households_clean.dta', '98 条清洗后家庭明细；同名 CSV 同时提供'),
    ('households_fitted.dta', '含拟合值、残差、残差平方的明细；同名 CSV 同时提供'),
    ('income_group_summary.csv', '收入组人数、均值、标准差、范围和拟合均值；另存 DTA'),
    ('source_totals_audit.csv', '明细合计与原表合计差异审计；另存 DTA'),
    ('regression_coefficients.csv', '常规系数、标准误、t 值、p 值及区间；另存 DTA'),
    ('robust_coefficients.csv', 'HC1 稳健系数表；另存 DTA'),
    ('model_statistics.csv', '拟合优度、残差核对与异方差检验统计量'),
    ('consumption_ols.ster', '可由 estimates use 重新载入的 OLS 估计结果'),
    ('consumption_robust.ster', '稳健标准误版本的估计结果'),
    ('build_week7_report.py', '独立验证与文档生成脚本'),
    ('verification.json', '数值验证状态及原始 Excel 的 SHA-256 指纹')
]
md.append('| 文件 | 内容 |\n| --- | --- |\n' + '\n'.join(f'| [{name}]({name}) | {desc} |' for name, desc in assets))
body.append('<div class="table-wrap"><table><thead><tr><th scope="col">文件</th><th scope="col">内容</th></tr></thead><tbody>'+''.join('<tr><td>'+asset_link(name,name)+'</td><td>'+html.escape(desc)+'</td></tr>' for name,desc in assets)+'</tbody></table></div>')
paragraph('原始 Excel 与参考 PNG 保留在 Eg2_1_1 文件夹中。Stata 批处理另会自动生成 week7_analysis.log。Markdown 需与图形文件保留相对位置；HTML 页面本身可单独打开阅读，其文件下载链接需相邻产物仍在原目录。')

heading('7. 方法依据与来源', 'sources')
paragraph('数值结果以本地原始 Excel、Stata 实际执行日志及独立数值核验为依据；Stata 指令含义核对了本机随软件附带的官方帮助文件。以下编号对应正文的方法引用：')
table(['编号', '来源', '使用内容'], [
    ['[1]', 'D:/Stata15/ado/base/r/regress.sthlp', 'OLS、vce(robust) 与 HC1 的 n/(n−k) 修正'],
    ['[2]', 'D:/Stata15/ado/base/r/regress_postestimation.sthlp', '线性预测 xb、残差 residuals 与 estat hettest'],
    ['[3]', 'D:/Stata15/ado/base/t/twoway_lfit.sthlp', '由一元线性回归预测值绘制拟合线'],
    ['[4]', 'Eg2_1_1/Eg2_1_1.xlsx 与 Eg2_1_1/Eg2_1_1.png', '原始数值、变量含义及仿图样式']
])
paragraph('随机扰动与残差的比较依据正文给出的总体条件均值定义及 OLS 正规方程推导。未将原表合计作为新增观测，也未假定已知真实总体参数。')

(ROOT/'Week7_回归分析报告.md').write_text('\n\n'.join(md)+'\n', encoding='utf-8')
style = '''
:root{color-scheme:light;--ink:#182c40;--muted:#5c6d7b;--accent:#146b72;--line:#dce5e9;--paper:#fff;--bg:#f3f6f7}*{box-sizing:border-box}html{scroll-behavior:smooth;scroll-padding-top:24px}body{margin:0;background:var(--bg);color:var(--ink);font:16px/1.85 "Microsoft YaHei","PingFang SC",system-ui,sans-serif}a{color:var(--accent);text-underline-offset:3px}header{background:#143044;color:white;padding:56px 32px 44px}.hero{max-width:1080px;margin:auto}.eyebrow{letter-spacing:3px;color:#a6d7d9;font-size:12px;font-weight:700}h1{font-size:clamp(28px,4.4vw,46px);line-height:1.4;margin:14px 0 16px;letter-spacing:-1px}.subhead{color:#c7d7e0;max-width:760px;margin:0}.metrics{display:grid;grid-template-columns:repeat(4,1fr);gap:14px;margin-top:30px}.metric{border-top:1px solid #547080;padding:14px 0 0}.metric span{display:block;font-size:12px;color:#bacfd9}.metric strong{font-size:27px;font-weight:600;font-variant-numeric:tabular-nums}.layout{max-width:1240px;margin:auto;display:grid;grid-template-columns:200px minmax(0,1fr);gap:38px;padding:34px 28px 70px}nav{position:sticky;top:24px;align-self:start;font-size:13px}nav .nav-title{color:var(--muted);margin:0 0 12px;font-size:11px;letter-spacing:2px}nav a{display:block;text-decoration:none;color:var(--muted);padding:8px 12px;border-left:2px solid var(--line)}nav a:hover,nav a:focus{border-color:var(--accent);color:var(--accent);background:#e7eff0}main{min-width:0;background:var(--paper);padding:30px 36px 42px;border:1px solid var(--line);border-radius:8px}h2{font-size:23px;line-height:1.5;margin:42px 0 18px;padding-top:24px;border-top:1px solid var(--line);scroll-margin-top:24px}p{margin:15px 0}strong{font-weight:650}.equation{font-family:"Cambria Math","Microsoft YaHei",serif;font-size:19px;line-height:1.8;background:#edf5f5;color:#154d58;border-left:3px solid #329097;padding:17px 22px;margin:22px 0;overflow-wrap:anywhere}.note,.warning{padding:17px 20px;background:#f0f5f8;border-radius:4px}.warning{background:#fff6e7;border-left:3px solid #c18a29}.table-wrap{overflow-x:auto;margin:20px 0;border:1px solid var(--line);border-radius:5px}table{width:100%;border-collapse:collapse;font-size:13px;font-variant-numeric:tabular-nums}th{text-align:left;background:#edf2f5;color:#304b5e;font-weight:600;white-space:nowrap}td,th{padding:11px 13px;border-bottom:1px solid var(--line);vertical-align:top}tbody tr:last-child td{border-bottom:0}tbody tr:nth-child(even){background:#fafcfc}td:first-child{white-space:nowrap}code{font:0.88em/1.6 Consolas,"Microsoft YaHei",monospace;background:#edf2f5;padding:2px 5px;border-radius:3px;overflow-wrap:anywhere}pre{background:#182d3d;color:#e5eff5;padding:20px;overflow-x:auto;border-radius:6px;font-size:13px;line-height:1.8}pre code{background:none;color:inherit;padding:0;white-space:pre}figure{margin:24px 0}img{display:block;width:100%;height:auto;border-radius:4px}figcaption{color:var(--muted);text-align:center;font-size:12px;margin-top:9px}details{padding:15px 18px;border:1px solid var(--line);border-radius:5px;margin:20px 0}summary{cursor:pointer;color:var(--accent);font-size:14px}details img{margin-top:20px}.prediction{background:#eff7f6;padding:20px 24px;border:1px solid #cbe3df;border-radius:6px;margin:24px 0}.prediction label{font-weight:600}.prediction-values{display:flex;justify-content:space-between;gap:16px;flex-wrap:wrap;margin-top:12px;font-size:14px}.prediction output{font-variant-numeric:tabular-nums}.prediction input{width:100%;margin:18px 0 0;accent-color:var(--accent)}.prediction p{font-size:12px;color:var(--muted);margin:6px 0 0}footer{text-align:center;padding:0 24px 36px;font-size:12px;color:var(--muted)}@media(max-width:980px){.layout{grid-template-columns:1fr;gap:18px;padding:22px 18px}nav{position:static;display:flex;flex-wrap:wrap;gap:4px}nav .nav-title{display:none}nav a{border-left:0;padding:5px 10px}main{padding:22px 25px}.metrics{gap:10px}.metric strong{font-size:22px}}@media(max-width:560px){header{padding:32px 22px}.metrics{grid-template-columns:repeat(2,1fr)}main{padding:18px 16px}.layout{padding:14px 10px 36px}h2{font-size:20px}.equation{font-size:17px;padding:14px}body{font-size:15px}td,th{padding:9px 10px}}@media print{body{background:white;color:black;font-size:10pt}header{background:white;color:black;padding:12px 0}.hero{max-width:none}.eyebrow,.subhead{color:#444}h1{font-size:24pt}.metric span,.metric strong{color:#222}.layout{display:block;padding:0;max-width:none}nav,.prediction,details,footer{display:none}main{border:0;padding:0}h2{break-after:avoid;margin-top:24px;padding-top:12px}table{font-size:9pt}tr,figure{break-inside:avoid}a{color:inherit}pre{white-space:pre-wrap;color:black;background:#f4f4f4}pre code{white-space:pre-wrap}.table-wrap{overflow:visible}}
'''
nav = '<nav aria-label="报告目录"><p class="nav-title">CONTENTS / 目录</p>'+''.join(f'<a href="#{anchor}">{html.escape(title)}</a>' for title,anchor in sections)+'</nav>'
page = f'''<!DOCTYPE html>
<html lang="zh-CN"><head><meta charset="utf-8"><meta name="viewport" content="width=device-width, initial-scale=1"><meta name="description" content="Week 7 计量经济学：随机扰动项与残差，98 条家庭消费明细的 Stata 回归、数据核对与复现报告。"><title>Week 7｜家庭消费回归分析</title><style>{style}</style></head>
<body><header><div class="hero"><div class="eyebrow">ECONOMETRICS · WEEK 07 · 2026.10.23</div><h1>随机扰动、残差<br>与家庭消费回归</h1><p class="subhead">从概念区分到数据核对，用 Stata 估计收入与消费的线性关系。完整呈现原始数据问题、拟合结果与可复现过程。</p><div class="metrics"><div class="metric"><span>有效家庭记录</span><strong>98</strong></div><div class="metric"><span>收入组数</span><strong>10</strong></div><div class="metric"><span>收入系数</span><strong>{model['slope']:.4f}</strong></div><div class="metric"><span>拟合优度 R²</span><strong>{model['r_squared']:.4f}</strong></div></div></div></header><div class="layout">{nav}<main>{''.join(body)}</main></div><footer>Stata/SE 15.1 实际运行 · Python 独立数值核对 · 页面图形与样式全部内嵌，可离线阅读</footer>
<script>
const slider = document.getElementById('income-slider');
slider.addEventListener('input', () => {{
 const income = Number(slider.value);
 document.getElementById('income-value').textContent = income.toFixed(0);
 document.getElementById('pred-value').textContent = ({model['intercept']!r} + {model['slope']!r} * income).toFixed(2);
}});
</script></body></html>'''
(ROOT/'Week7_回归分析报告.html').write_text(page, encoding='utf-8')
# Static checks for complete, portable HTML and working local artifact links.
from html.parser import HTMLParser
class ValidateHTML(HTMLParser):
    def handle_starttag(self, tag, attrs):
        attrs = dict(attrs)
        if tag == 'img':
            assert attrs['src'].startswith('data:image/png;base64,')
            assert attrs.get('alt')
        if tag == 'a' and not attrs['href'].startswith('#'):
            assert (ROOT / attrs['href']).exists(), attrs['href']
ValidateHTML().feed(page)
assert '\ufffd' not in page
assert len(re.findall(r'<h2 ', page)) == 7
print(json.dumps({'verification': 'passed', 'records': n, 'report_markdown': str(ROOT/'Week7_回归分析报告.md'), 'report_html': str(ROOT/'Week7_回归分析报告.html'), 'max_abs_fitted_difference': verification['max_abs_fitted_difference']}, ensure_ascii=False, indent=2))