Week7-1

Week7/build_week7_report.py

@@ -0,0 +1,275 @@
+"""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))