import numpy as np
import pandas as pd
import statsmodels.api as sm

# 生成模拟数据
np.random.seed(123)
n = 100
X = np.random.randn(n, 3)
X = sm.add_constant(X)  # 添加常数项
beta_true = np.array([1.5, 2.0, -1.0, 0.5])
epsilon = np.random.randn(n) * 0.5
y = X @ beta_true + epsilon

# 方法1：矩阵直接计算
def ols_matrix_method(X, y):
    """使用矩阵公式计算OLS系数"""
    XTX_inv = np.linalg.inv(X.T @ X)
    beta_hat = XTX_inv @ X.T @ y
    return beta_hat

# 方法2：使用statsmodels
def ols_statsmodels(X, y):
    """使用statsmodels进行OLS回归"""
    model = sm.OLS(y, X)
    results = model.fit()
    return results

# 比较结果
beta_matrix = ols_matrix_method(X, y)
model_sm = ols_statsmodels(X, y)
beta_sm = model_sm.params

print("真实参数:", beta_true)
print("矩阵法估计:", beta_matrix)
print("Statsmodels估计:", beta_sm.values)
print("\n系数比较（最大绝对差）:", np.max(np.abs(beta_matrix - beta_sm)))

# 加载包
library(lmtest)
library(sandwich)
library(ggplot2)

# 示例数据
set.seed(123)
n <- 100
x <- runif(n, 1, 10)
# 异方差数据生成：方差随x增大而增大
y <- 2 + 3*x + rnorm(n, 0, 0.5*x)

# 1. OLS回归
model <- lm(y ~ x)
summary(model)

# 2. 异方差检验
# White检验（一般形式）
white_test <- bptest(model, ~ x + I(x^2))
cat("White检验p值:", white_test$p.value, "\n")

# BP检验
bp_test <- bptest(model)
cat("BP检验p值:", bp_test$p.value, "\n")

# 残差图诊断
residuals_df <- data.frame(
  fitted = fitted(model),
  residuals = residuals(model),
  x = x
)

ggplot(residuals_df, aes(x = fitted, y = residuals)) +
  geom_point(alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "残差 vs 拟合值图", 
       subtitle = "漏斗形状提示异方差") +
  theme_minimal()

# 3. 异方差稳健标准误
# HC0-HC4类型稳健标准误
coeftest(model, vcov = vcovHC(model, type = "HC0"))
coeftest(model, vcov = vcovHC(model, type = "HC3"))  # 小样本推荐

# 4. 加权最小二乘法(WLS)
# 假设方差与x成比例
w <- 1/x^2  # 权重与方差成反比
model_wls <- lm(y ~ x, weights = w)
summary(model_wls)

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.sandbox.regression.gmm import IV2SLS
import linearmodels.iv as iv
import warnings
warnings.filterwarnings('ignore')

# 生成模拟数据（存在内生性）
np.random.seed(123)
n = 500

# 工具变量z，与内生变量x相关，但与误差项不相关
z = np.random.normal(0, 1, n)

# 内生变量x（受z和误差项v影响）
v = np.random.normal(0, 1, n)  # x的误差项
x = 0.5 + 0.8*z + v  # 第一阶段：z影响x

# 结构方程误差项，与v相关（导致内生性）
epsilon = 0.7*v + np.random.normal(0, 1, n)  # 内生性来源：Cov(x, epsilon) ≠ 0

# 因变量y
y = 2 + 3*x + epsilon

# 创建DataFrame
df = pd.DataFrame({
    'y': y,
    'x': x,
    'z': z,
    'const': 1
})

# 1. OLS估计（有偏）
ols_model = sm.OLS(df['y'], df[['const', 'x']]).fit()
print("=== OLS估计结果（有偏）===")
print(ols_model.summary().tables[1])

# 2. 2SLS估计（使用statsmodels）
# 方法1：使用IV2SLS
from statsmodels.sandbox.regression.gmm import IV2SLS
iv_model_sm = IV2SLS(df['y'], df[['const']], df[['x']], df[['const', 'z']]).fit()
print("\n=== 2SLS估计结果（statsmodels）===")
print(iv_model_sm.summary())

# 方法2：使用linearmodels（更专业的计量经济学包）
from linearmodels.iv import IV2SLS
formula = 'y ~ 1 + [x ~ z]'  # 1表示常数项，[x ~ z]表示x是内生变量，z是工具变量
iv_model_lm = IV2SLS.from_formula(formula, df).fit(cov_type='robust')
print("\n=== 2SLS估计结果（linearmodels）===")
print(iv_model_lm.summary)

# 3. 第一阶段回归
first_stage = sm.OLS(df['x'], df[['const', 'z']]).fit()
print("\n=== 第一阶段回归 ===")
print("系数估计:", first_stage.params)
print("F统计量:", first_stage.fvalue)
print("R-squared:", first_stage.rsquared)

# 4. 弱工具变量检验（经验法则：F > 10）
print("\n=== 弱工具变量检验 ===")
print(f"第一阶段F统计量: {first_stage.fvalue:.2f}")
print("判断：F > 10" if first_stage.fvalue > 10 else "警告：可能存在弱工具变量问题")

# 5. Hausman内生性检验
# 比较OLS和2SLS系数差异
from scipy.stats import chi2
beta_ols = ols_model.params['x']
beta_iv = iv_model_lm.params['x']
var_ols = ols_model.cov_params().loc['x', 'x']
var_iv = iv_model_lm.cov.loc['x', 'x']

hausman_stat = (beta_iv - beta_ols)**2 / (var_iv - var_ols)
p_value = 1 - chi2.cdf(hausman_stat, 1)
print(f"\n=== Hausman内生性检验 ===")
print(f"统计量: {hausman_stat:.4f}, p值: {p_value:.4f}")
print("结论：拒绝原假设（存在内生性）" if p_value < 0.05 else "结论：不能拒绝原假设")