计量经济学(二)——多元线性回归模型

多元线性回归(Multiple Linear Regression, MLR)是一种统计模型,被广泛认为是计量经济学的核心基础。多元线性回归为经济研究者提供了一种有效的方法来建模和分析多个自变量与因变量之间的线性关系。

在计量经济学中,研究者常常面临复杂的经济现象,这些现象往往受多种因素影响。通过建立多元线性回归模型,研究者能够定量分析各个自变量对因变量的影响程度,从而揭示经济关系的本质。例如,在分析工资水平时,研究者可能会考虑教育、工作经验、行业和地区等多个因素,通过多元线性回归,能够有效识别这些因素的边际贡献。
多元线性回归不仅可以用于模型构建,还可以用于假设检验,检验经济理论的有效性。计量经济工作者可以使用回归分析来验证某一理论预测是否成立,如需求弹性或生产函数的形式。通过统计检验(如t检验和F检验),研究者可以判断自变量对因变量的影响是否显著,进而为经济理论提供实证支持。
多元线性回归具有良好的统计性质,如无偏性和一致性,使得其在经济数据分析中成为一种稳健的工具。随着数据分析技术的发展,多元线性回归也可以与现代机器学习方法结合,扩展其应用范围。通过对模型进行适当的调整和扩展,研究者能够应对非线性关系、多重共线性、异方差性等问题,从而提高模型的准确性和适用性。

程序求解 多元回归几何图

一、多元线性回归模型

在面对复杂的非线性系统时,将其转化为线性系统进行处理的原因和考量主要集中在以下几个方面。
简化计算与模型求解:线性模型的计算相对简单,参数估计通常使用最小二乘法,可以直接求解。而对于非线性模型,通常需要迭代算法(如牛顿法、梯度下降法)来估计参数,这些方法计算复杂且容易陷入局部最优解。将非线性问题通过适当的变换处理成线性问题,可以利用现成的线性回归方法,更快速、稳定地得到解。
稳健性与理论基础:线性回归有着良好的统计学性质,如无偏性、有效性和一致性。这些性质使得线性模型在估计和预测中表现得更加稳健。相比之下,非线性模型因为更复杂,容易受到过拟合、参数估计不准确等问题的影响。通过线性化处理,研究者可以保留这些稳健性特征,减少模型不稳定的风险。
易于解释和理解:线性模型的结果易于解释,回归系数表示每个自变量对因变量的边际贡献,直观而且容易被理解。而非线性模型通常难以直接解释自变量与因变量之间的关系。将非线性问题转化为线性模型,可以让研究者和读者更清楚地理解变量之间的相互作用。

1.1 模型假设

多元线性回归模型的基本形式为:

\[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_k X_k + \varepsilon \]

其中:

  • $ Y$ 是因变量,表示我们希望解释或预测的结果。
  • $ X_1, X_2, \cdots, X_k$ 是自变量,表示影响因变量的因素。
  • $ \beta_0$ 是截距项,$ \beta_1, \beta_2, \cdots, \beta_k$ 是各个自变量的回归系数。
  • $ \varepsilon$ 是随机误差项,表示无法通过自变量解释的部分。

为了保证回归模型能够给出有效的估计结果,必须满足以下几个经典假设:

  • 线性关系:因变量与自变量之间的关系是线性的,即因变量是自变量的线性组合。引出非线性
  • 独立性:误差项彼此之间是相互独立的,不存在相关性。引出自相关等
  • 同方差性:误差项的方差是恒定的,不随自变量的取值变化而变化。引出异方差性
  • 正态性:误差项服从正态分布,尤其在样本量较小时这一假设尤为重要。引出非正态Logit回归等
  • 无共线性:自变量之间没有完全的线性相关性,即不存在多重共线性问题。引出多重共线性

1.2 参数估计

多元线性回归的参数估计方法通常采用普通最小二乘法(Ordinary Least Squares, OLS),这是最经典的估计方法之一,其目标是最小化回归模型中的残差平方和,即预测值与实际观测值之间差距的平方和最小化。OLS 在经典线性回归模型的框架下表现出优良的统计性质,如无偏性、最小方差性和一致性,使得它成为经济学、社会科学等领域分析变量关系的重要工具。

普通最小二乘法的原理
在多元线性回归模型中,我们假设因变量 $ Y $ 与一组自变量 $ X $ 之间的关系可以用以下线性方程表示:

\[Y = X \beta + \varepsilon \]

其中:

  • $ Y $ 是因变量的向量;
  • $ X $ 是 $ n \times k $ 维的设计矩阵,包含 $ n $ 个样本和 $ k $ 个自变量;
  • $ \beta $ 是回归系数向量;
  • $ \varepsilon $ 是误差项向量。

普通最小二乘法的核心思想是通过最小化模型预测值 $ X \beta $ 与实际观测值 $ Y $ 之间的差距,即最小化残差平方和来估计回归系数 $ \beta $。残差定义为 $ \varepsilon = Y - X \beta $,因此残差平方和 $ S(\beta) $ 为:

\[S(\beta) = \sum_{i=1}^{n} \varepsilon_i^2 = (Y - X \beta)^T (Y - X \beta) \]

最小化这个残差平方和的过程可以通过求解偏导数得到最优的 $ \beta $。该优化问题的解是:

\[\hat{\beta} = (X^T X)^{-1} X^T Y \]

其中,$ \hat{\beta} $ 是估计的回归系数向量,$ X^T X $ 是设计矩阵 $ X $ 的转置矩阵,$ (X^T X)^{-1} $ 是矩阵 $ X^T X $ 的逆矩阵。这个方程是 OLS 的闭式解,用于计算最佳的回归系数。

普通最小二乘法的性质
在经典线性回归模型的假设条件下,OLS 估计具有以下三个主要性质:

  • 无偏性:OLS 估计是无偏的,即在重复抽样的条件下,OLS 的估计值 $ \hat{\beta} $ 的期望等于真实值 $ \beta $。这意味着,尽管单个样本的估计可能有误差,但在大量重复样本中,估计结果平均来说是准确的。
  • 最小方差性:根据高斯-马尔科夫定理,在所有线性无偏估计量中,OLS 具有最小的方差。这意味着在满足模型假设的条件下,OLS 是最有效的线性估计量,它能够提供最小的估计误差。
  • 一致性:随着样本量的增加,OLS 估计量趋近于真实参数值,即估计是“收敛”的。无论样本量如何大小,OLS 估计量在大样本情况下会更接近真实的参数。
    这些性质使 OLS 成为统计和计量经济学中广泛应用的估计方法。

1.3 常见模型问题及其解决方法

异方差性
异方差性是指模型中误差项的方差不是恒定的,即残差的波动随着自变量的变化而变化。这违反了经典线性回归模型的同方差性假设,导致回归系数的标准误估计值失效,进而影响t检验和F检验的准确性。
异方差性的检验

  • 图形法:通过绘制残差与拟合值的散点图,观察是否存在系统性的模式。如果残差图中出现漏斗形或其他形态,则可能存在异方差性。
  • 布鲁斯-帕根检验(Breusch-Pagan Test):该检验通过检验残差平方是否与自变量存在关联,来判断是否存在异方差性。
  • 怀特检验(White Test):通过构建自变量及其平方项来检验模型中的异方差性。
  • Goldfeld-Quandt检验:该检验方法通过将样本按某一自变量排序并分成两组,然后比较两组的残差方差。如果方差存在显著差异,表明存在异方差。

消除异方差性的方法

  • 加权最小二乘法(WLS):通过对不同观测值赋予不同的权重,以解决误差项方差不等的问题。加权较大的值代表方差较小的观测,降低其对整体估计的影响。
  • 稳健标准误(Robust Standard Errors):若难以消除异方差性,稳健标准误可以在异方差性存在的情况下给出更为准确的标准误估计,保证显著性检验的有效性。
  • 广义最小二乘法(GLS):GLS 是一种扩展的OLS方法,能够处理误差项的异方差性问题。

多重共线性
多重共线性是指模型中的自变量之间存在较高的线性相关性,使得回归系数无法准确估计。共线性问题不会影响模型的拟合优度,但会导致回归系数的标准误增大,影响自变量的显著性检验,甚至会使某些重要变量被错误地认为是不显著的。

多重共线性的检验
方差膨胀因子
方差膨胀因子是检测多重共线性的最常用方法之一。它衡量某个自变量 $ X_i $ 与其他所有自变量之间的线性相关性程度,具体公式为:

\[\text{VIF}_i = \frac{1}{1 - R_i^2} \]

其中,$ R_i^2 $ 是将自变量 $ X_i $ 对其他自变量进行回归后的决定系数。如果 $ R_i^2 $ 较高,说明 $ X_i $ 与其他自变量之间存在较强的线性关系,VIF 值也会随之增大。

判断依据

  • 当 $ \text{VIF}_i \approx 1 $ 时,表明没有多重共线性;
  • 当 $ \text{VIF}_i > 5 $ 时,通常认为存在中等程度的多重共线性;
  • 当 $ \text{VIF}_i > 10 $ 时,认为存在严重的多重共线性。

VIF 提供了一个简单易用的工具来判断某个特定变量是否与其他自变量存在共线性问题。如果所有自变量的 VIF 值都较高,说明回归模型整体上可能存在共线性问题。

特征值分解法(Eigenvalue Decomposition)和条件数(Condition Number)
特征值分解法是通过分析设计矩阵 $ X $ 的特征值来判断多重共线性的方法。特征值表示的是矩阵中的线性独立信息的多少,特征值越小,说明矩阵中的变量越不独立,存在多重共线性。在特征值分解的基础上,条件数(Condition Number)常用于进一步判定共线性。条件数定义为设计矩阵 $ X $ 的最大特征值 $ \lambda_{\text{max}} $ 与最小特征值 $ \lambda_{\text{min}} $ 的比值:

\[\kappa(X) = \frac{\lambda_{\text{max}}}{\lambda_{\text{min}}} \]

判断依据

  • 当 $ \kappa(X) < 10 $ 时,说明没有显著的多重共线性;
  • 当 $ \kappa(X) > 30 $ 时,通常认为存在严重的多重共线性。

特征值分解法是一种更加数学化的方法,适合用于深入分析多重共线性问题,特别是在有多个自变量的情况下。

相关矩阵法(Correlation Matrix)
相关矩阵法是通过检查自变量之间的相关系数来识别共线性。计算每对自变量 $ X_i $ 和 $ X_j $ 之间的皮尔逊相关系数 $ \rho_{ij}$:

\[\rho_{ij} = \frac{\text{Cov}(X_i, X_j)}{\sigma_{X_i} \sigma_{X_j}} \]

其中,$ \text{Cov}(X_i, X_j) $ 是 $ X_i $ 和 $ X_j $ 之间的协方差,$ \sigma_{X_i} $ 和 $ \sigma_{X_j} $ 分别是自变量的标准差。

判断依据

  • 当 $ |\rho_{ij}| > 0.8 $ 时,通常认为这对自变量之间存在较强的共线性;
  • 当 $ |\rho_{ij}| > 0.9 $ 时,说明存在严重的共线性。

相关矩阵法是识别共线性的初步手段,可以帮助判断变量之间的线性关系,但它无法解决变量与其他多个变量之间的线性相关性问题,因此需要与其他方法配合使用。

条件指数(Condition Index)
条件指数是通过特征值分析进一步量化共线性程度的指标。它基于特征值分解,将特征值的平方根用于衡量:

\[\text{Condition Index} = \frac{\sqrt{\lambda_{\text{max}}}}{\sqrt{\lambda_i}} \]

其中,$ \lambda_i $ 是第 $ i $ 个特征值。条件指数反映了设计矩阵的不同维度之间的线性相关性。

判断依据

  • 当条件指数大于 10 时,表明存在一定程度的共线性;
  • 当条件指数大于 30 时,认为存在严重的多重共线性。

条件指数与条件数相似,但它通过对各个特征值进行逐步分解,能够进一步帮助识别出具体的变量之间共线性的程度。

逐步回归法(Stepwise Regression)
在多重共线性情况下,某些自变量的显著性检验可能会受到影响。因此,逐步回归法(Stepwise Regression)可以通过引入或剔除变量来检测和处理共线性问题。逐步回归法逐步将自变量引入回归模型或将不显著的变量移除,并在每一步都检测回归模型的显著性和自变量之间的相关性。虽然逐步回归法不能直接量化共线性,但它可以通过调整模型结构来减弱共线性的影响。此外,逐步回归法可以与 VIF 等方法结合使用,以提高模型的稳定性。

消除多重共线性的方法
删除高度相关的自变量:如果某些自变量之间高度相关,可以通过删除其中一个变量来减少共线性。
主成分回归(PCR):通过主成分分析将自变量降维,将相关性较强的变量转化为若干无相关的主成分,再用这些主成分进行回归分析。
岭回归(Ridge Regression):在回归模型中引入正则化惩罚项,通过控制回归系数的大小,减少共线性对模型的影响。

二、案例分析1

考虑一家快递公司送货:X1: 运输里程 X2: 运输次数 Y:总运输时间,数据见下表所示。

ID X1: Miles Traveled X2: Number of Deliveries Y: Travel Time (Hours)
1 100 4 9.3
2 50 3 4.8
3 100 4 8.9
4 100 2 6.5
5 50 2 4.2
6 80 2 6.2
7 75 3 7.4
8 65 4 6.0
9 90 3 7.6
10 90 2 6.1
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.outliers_influence import variance_inflation_factor

# 数据输入
data = {
    'X1': [100, 50, 100, 100, 50, 80, 75, 65, 90, 90],  # Miles traveled
    'X2': [4, 3, 4, 2, 2, 2, 3, 4, 3, 2],  # Number of deliveries
    'Y': [9.3, 4.8, 8.9, 6.5, 4.2, 6.2, 7.4, 6.0, 7.6, 6.1]  # Travel time
}

# 创建DataFrame
df = pd.DataFrame(data)

# 自变量(X1 和 X2)
X = df[['X1', 'X2']]
X = sm.add_constant(X)  # 添加常数项

# 因变量(Y)
Y = df['Y']

# OLS回归
model = sm.OLS(Y, X).fit()

# 输出回归结果
summary_table = model.summary2()
print(summary_table)

# Breusch-Pagan 异方差检验
bp_test = het_breuschpagan(model.resid, X)
bp_test_results = {
    'Lagrange multiplier statistic': bp_test[0],
    'p-value': bp_test[1],
    'f-value': bp_test[2],
    'f p-value': bp_test[3]
}
print("\nBreusch-Pagan Test Results:")
for key, value in bp_test_results.items():
    print(f"{key}: {value:.4f}")

# 方差膨胀因子(VIF) 共线性检验
vif = pd.DataFrame()
vif["Variable"] = X.columns
vif["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]

print("\nVariance Inflation Factor (VIF):")
print(vif)

OLS 回归结果:

Dependent Variable: Y
Model: OLS
R-squared: 0.904
Adj. R-squared: 0.876
No. Observations: 10
F-statistic: 32.88
Prob (F-statistic): 0.000276
Log-Likelihood: -6.8398
AIC: 19.6796
BIC: 20.5873

Variable Coefficient Std. Err. t-statistic P-value 95% Confidence Interval
const -0.8687 0.9515 -0.9129 0.3916 [-3.1188, 1.3814]
X1 0.0611 0.0099 6.1824 0.0005 [0.0378, 0.0845]
X2 0.9234 0.2211 4.1763 0.0042 [0.4006, 1.4463]

Durbin-Watson: 2.515
Condition No.: 435

Breusch-Pagan 异方差检验结果:

  • Lagrange Multiplier Statistic: 0.369
  • p-value: 0.8315 (无异方差问题)

方差膨胀因子(VIF):

Variable VIF
const 27.5636
X1 1.027
X2 1.027

该模型表明没有显著的异方差和多重共线性问题。

三、案例分析2

这里采用波士顿房价数据集进行多元回归分析,各个特征变量(自变量)含义如下表所示。

Attribute Description Attribute Description
CRIM 该区域内犯罪率(每人均犯罪次数) DIS 该区域内房屋到市中心的加权距离
ZN 该区域内住宅地的占地面积比例 RAD 该区域内房屋到径向高速公路的距离指数
INDUS 该区域内非零售商业用地比例 TAX 该区域内每一万美元的财产税率
CHAS 是否靠河,1为靠河,0为不靠河 PTRATIO 该区域内教师与学生数之比
NOX 一氧化氮浓度(每千万分之一) B 该区域内黑人所占比例
RM 该区域内房屋的平均房间数量 LSTAT 该区域内人口中有多少比例属于低收入阶层
AGE 该区域内房屋的平均房龄 PRICE 自主房屋房价中位数(被解释变量-因变量)

3.1 描述性分析

# 导入必要的库
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

# 读取波士顿房价数据
data_url = "http://lib.stat.cmu.edu/datasets/boston"
raw_df = pd.read_csv(data_url, sep="\s+", skiprows=22, header=None)

# 数据处理
data = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]])  # 特征
target = raw_df.values[1::2, 2]  # 房价(目标变量)

# 创建DataFrame
columns = ['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']
df = pd.DataFrame(data, columns=columns)
df['PRICE'] = target  # 添加房价列

# 计算相关矩阵
corr_matrix = df.corr()

# 绘制相关矩阵的热力图
plt.figure(figsize=(12, 8))
sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', linewidths=0.5)
plt.title('Correlation Heatmap of Boston Housing Dataset')
plt.show()

可以看出RM,lSTAT,INDUS,PTRATIO都是与PRICE的相关性比较大的特征。

3.2 多元回归分析

波士顿房价数据集是20世纪70年代波士顿郊区房价的中位数,统计了当时教区部分的犯罪率、房产税等共计13个输入变量,统计出房价(1个输出变量)。哪些指标对房价的影响较大?通过对波士顿房价数据集的作图和定性分析,可以得出以下结论:

  • 房屋的平均房间数 (RM):对房价的影响较大,系数为正数,说明房间越多,房价越高。
  • 距离五个波士顿就业中心的加权距离 (DIS):对房价的影响较大,系数为负数,说明距离就业中心越远,房价越低。
  • 教育水平较高的人口所占比例 (LSTAT):对房价的影响较大,系数为负数,说明该比例越高,房价越低。
  • 住宅的年龄 (AGE):对房价的影响较大,系数为负数,说明住宅越老,房价越低。

其他自变量对房价的影响较小,甚至为负数,说明这些因素与房价呈反比例关系。为简化计算仅对上面四个变量做回归,需要注意的是,这样的结论仅代表波士顿房价数据集的情况,不一定适用于其他地区。

# 导入必要的库
import pandas as pd
import numpy as np
import statsmodels.api as sm
import seaborn as sns
import matplotlib.pyplot as plt

# 读取波士顿房价数据
data_url = "http://lib.stat.cmu.edu/datasets/boston"
raw_df = pd.read_csv(data_url, sep="\s+", skiprows=22, header=None)

# 数据处理
data = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]])  # 特征
target = raw_df.values[1::2, 2]  # 房价(目标变量)

# 创建DataFrame
columns = ['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']
df = pd.DataFrame(data, columns=columns)
df['PRICE'] = target  # 添加房价列

# 选择变量:RM, DIS, LSTAT, AGE 和 PRICE
X = df[['RM', 'DIS', 'LSTAT', 'AGE']]  # 自变量
y = df['PRICE']  # 因变量

# 添加常数项(截距项)
X = sm.add_constant(X)

# 使用OLS进行回归分析
model = sm.OLS(y, X).fit()

# 输出回归结果
print(model.summary())

# 绘制残差图
plt.figure(figsize=(10,6))
sns.residplot(x=model.fittedvalues, y=model.resid, lowess=True, line_kws={'color': 'red'})
plt.xlabel('Fitted values')
plt.ylabel('Residuals')
plt.title('Residuals vs Fitted values')
plt.show()

# 绘制预测值与实际值对比图
plt.figure(figsize=(10,6))
plt.scatter(model.fittedvalues, y, alpha=0.5, color='b')
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'r--', lw=2)
plt.xlabel('Predicted')
plt.ylabel('Actual')
plt.title('Predicted vs Actual')
plt.show()
  OLS Regression Results                            
==============================================================================
Dep. Variable:                  PRICE   R-squared:                       0.649
Model:                            OLS   Adj. R-squared:                  0.646
Method:                 Least Squares   F-statistic:                     231.6
Date:                Mon, 14 Oct 2024   Prob (F-statistic):          2.09e-112
Time:                        14:41:23   Log-Likelihood:                -1575.4
No. Observations:                 506   AIC:                             3161.
Df Residuals:                     501   BIC:                             3182.
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.0799      3.437      1.187      0.236      -2.672      10.832
RM             4.9904      0.449     11.126      0.000       4.109       5.872
DIS           -0.6588      0.175     -3.767      0.000      -1.002      -0.315
LSTAT         -0.6848      0.054    -12.721      0.000      -0.791      -0.579
AGE           -0.0255      0.014     -1.771      0.077      -0.054       0.003
==============================================================================
Omnibus:                      124.062   Durbin-Watson:                   0.869
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              355.890
Skew:                           1.171   Prob(JB):                     5.24e-78
Kurtosis:                       6.376   Cond. No.                     1.08e+03
==============================================================================

结果解释

  • R-squared (R²):
    • 值为 0.649,表示模型解释了约 64.9% 的房价变动。这是一个中等的拟合度,表明还有 35.1% 的变动未被模型解释。
  • Adj. R-squared (调整后的 R²):
    • 值为 0.646,相较于 R²,这个指标考虑了自变量的数量。调整后的 R² 用于评估模型的拟合度,当增加自变量时,调整后的 R² 可能会降低。
  • F-statistic:
    • 值为 231.6,表明模型整体是显著的(P值为 2.09e-112)。这意味着至少有一个自变量与因变量 PRICE 之间存在显著关系。
  • 各个自变量的系数 (coef):
    • RM (房间数): 系数为 4.9904,表明房间数每增加一个单位,房价将增加约 4.99 千美元,且 p 值为 0.000,说明显著。
    • DIS (与就业中心的距离): 系数为 -0.6588,说明距离每增加一个单位,房价将降低约 0.66 千美元,且 p 值也显著。
    • LSTAT (低收入人群比例): 系数为 -0.6848,表明比例每增加一个单位,房价将降低约 0.68 千美元,且 p 值显著。
    • AGE (房屋年龄): 系数为 -0.0255,虽然 p 值为 0.077(接近 0.05),但未达到显著水平,表明房屋年龄对房价的影响不明显。
  • Omnibus 和 Jarque-Bera 检验:
    • 这些检验用于检查残差的正态性。Omnibus 的 p 值为 0.000,表明残差不服从正态分布,Jarque-Bera 检验的 p 值也非常小,这进一步支持了这一结论。
  • Durbin-Watson:
    • 值为 0.869,接近于 0,表明可能存在自相关。理想情况下,Durbin-Watson 值应在 1.5 至 2.5 之间。
  • Condition Number (条件数):
    • 值为 1.08e+03,条件数大于 30 时,表明模型中可能存在多重共线性问题。

3.3 检查结果

import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.diagnostic import het_breuschpagan, acorr_breusch_godfrey
from sklearn.datasets import fetch_california_housing

# 加载数据集:使用加州房价数据作为替代
housing = fetch_california_housing(as_frame=True)
df = housing.frame

# 选择自变量和因变量
X = df[['MedInc', 'HouseAge', 'AveRooms', 'AveOccup']]
y = df['MedHouseVal']

# 添加常数项
X = sm.add_constant(X)

# 拟合 OLS 回归模型
model = sm.OLS(y, X).fit()

# 打印回归结果
print(model.summary())

# 1. 共线性检测(VIF)
def calculate_vif(X):
    vif = pd.DataFrame()
    vif["Variables"] = X.columns
    vif["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
    return vif

vif = calculate_vif(X)
print("\n共线性检测结果 (VIF):")
print(vif)

# 2. 异方差性检测 (Breusch-Pagan 检验)
bp_test = het_breuschpagan(model.resid, model.model.exog)
bp_test_results = pd.Series(bp_test, index=['LM Stat', 'LM p-value', 'F-stat', 'F p-value'])
print("\n异方差性检测结果 (Breusch-Pagan Test):")
print(bp_test_results)

# 3. 自相关性检测 (Breusch-Godfrey 检验)
bg_test = acorr_breusch_godfrey(model, nlags=1)
bg_test_results = pd.Series(bg_test, index=['LM Stat', 'LM p-value', 'F-stat', 'F p-value'])
print("\n自相关性检测结果 (Breusch-Godfrey Test):")
print(bg_test_results)

检查共线性:VIF 值超过 10 通常表明存在严重的多重共线性。
检查异方差性:Breusch-Pagan 检验的 p 值如果小于 0.05,表明存在异方差性。
检查自相关:Breusch-Godfrey 检验的 p 值如果小于 0.05,说明存在自相关。

有待于结合使用其他模型展开分析。

四、道格拉斯生产函数的多元回归

道格拉斯生产函数通常表示为:$$Y = A \cdot L^{\alpha} \cdot K^{\beta}$$
其中:

  • \(Y\)表示总产值,
  • \(L\)表示劳动投入(工资总额),
  • \(K\)表示资本投入(投资总额),
  • \(A\)是技术系数,
  • \(\alpha\)\(\beta\)分别是劳动和资本的弹性系数,即贡献率

弹性系数的解释

  • α (Alpha):劳动要素的产出弹性,表示劳动对产出的贡献率。它表明当劳动投入增加 1% 时,总产出将增加\(\alpha\)%。
  • β (Beta):资本要素的产出弹性,表示资本对产出的贡献率。它表明当资本投入增加 1% 时,总产出将增加\(\beta\)%。

贡献率的含义

  • α + β = 1:如果这两个系数之和为 1,表明生产过程存在规模报酬不变,即劳动和资本投入比例增加时,产出也会成比例增加。
  • α + β > 1:则表示存在规模报酬递增,即投入的增长会带来更大比例的产出增长。
  • α + β < 1:则表示存在规模报酬递减,即投入的增长带来的产出增长比例相对较小。

为了将其转换为可以使用普通最小二乘法(OLS)进行回归分析的形式,我们对等式两边取对数,得到:

\[\ln Y = \ln A + \alpha \ln L + \beta \ln K \]

因此,这是一个线性回归模型,其中:\(\ln Y\)是因变量,\(\ln L\)\(\ln K\)是自变量。

下表为1990-2011年上海市工业企业劳动工资和资本支出以及总产值的数据,可以分别回归出劳动要素和资本要素的贡献率。

年份 工资总额 (L) 投资总额 (K) 总产值 (Y) 年份 工资总额 (L) 投资总额 (K) 总产值 (Y)
1990 146.78 227.08 781.66 2001 678.29 1994.73 5210.12
1991 172.84 258.3 893.77 2002 733.31 2187.06 5741.03
1992 217.21 357.38 1114.32 2003 803.84 2452.11 6694.23
1993 279.33 653.91 1519.23 2004 837.39 3084.66 8072.83
1994 357.89 1123.29 1990.86 2005 1146.97 3542.55 9247.66
1995 440.75 1601.79 2499.43 2006 1475.93 3925.09 10572.24
1996 492.7 1952.05 2957.55 2007 1802.17 4458.61 12494.01
1997 510.1 1977.59 3438.79 2008 2184.2 4829.45 14069.87
1998 510.35 1964.83 3801.09 2009 2594.2 5273.33 15046.45
1999 583.54 1856.72 4188.73 2010 3018.55 5317.67 17165.98
2000 614.53 1869.67 4771.17 2011 4505.57 5067.09 19195.69
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

# 创建数据
data = {
    'Year': np.arange(1990, 2012),
    'L': [146.78, 172.84, 217.21, 279.33, 357.89, 440.75, 492.7, 510.1, 510.35, 583.54, 614.53, 678.29, 733.31, 803.84, 837.39, 1146.97, 1475.93, 1802.17, 2184.2, 2594.2, 3018.55, 4505.57],
    'K': [227.08, 258.3, 357.38, 653.91, 1123.29, 1601.79, 1952.05, 1977.59, 1964.83, 1856.72, 1869.67, 1994.73, 2187.06, 2452.11, 3084.66, 3542.55, 3925.09, 4458.61, 4829.45, 5273.33, 5317.67, 5067.09],
    'Y': [781.66, 893.77, 1114.32, 1519.23, 1990.86, 2499.43, 2957.55, 3438.79, 3801.09, 4188.73, 4771.17, 5210.12, 5741.03, 6694.23, 8072.83, 9247.66, 10572.24, 12494.01, 14069.87, 15046.45, 17165.98, 19195.69]
}

# 转换为DataFrame
df = pd.DataFrame(data)

# 取对数
df['log_Y'] = np.log(df['Y'])
df['log_L'] = np.log(df['L'])
df['log_K'] = np.log(df['K'])

# 设置自变量和因变量
X = df[['log_L', 'log_K']]
X = sm.add_constant(X)  # 添加常数项
y = df['log_Y']

# 进行回归分析
model = sm.OLS(y, X).fit()

# 输出回归结果
print(model.summary())

# 绘制残差图
plt.figure(figsize=(8, 6))
plt.plot(df['Year'], model.resid, marker='o')
plt.axhline(0, color='r', linestyle='--')
plt.xlabel('Year')
plt.ylabel('Residuals')
plt.title('Residual Plot')
plt.grid(True)
plt.show()
 OLS Regression Results                            
==============================================================================
Dep. Variable:                  log_Y   R-squared:                       0.981
Model:                            OLS   Adj. R-squared:                  0.978
Method:                 Least Squares   F-statistic:                     478.0
Date:                Mon, 14 Oct 2024   Prob (F-statistic):           5.66e-17
Time:                        17:31:59   Log-Likelihood:                 13.382
No. Observations:                  22   AIC:                            -20.76
Df Residuals:                      19   BIC:                            -17.49
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.1834      0.250      4.726      0.000       0.659       1.708
log_L          0.6285      0.090      7.021      0.000       0.441       0.816
log_K          0.4159      0.088      4.729      0.000       0.232       0.600
==============================================================================
Omnibus:                        0.870   Durbin-Watson:                   0.297
Prob(Omnibus):                  0.647   Jarque-Bera (JB):                0.870
Skew:                          -0.354   Prob(JB):                        0.647
Kurtosis:                       2.329   Cond. No.                         85.3
==============================================================================

结果解释
异方差性:需要通过 Breusch-Pagan 检验确认,但结果并未直接显示异方差问题。
共线性:通过 VIF 检测共线性,当前系数和标准误差情况看,不存在严重共线性。
自相关性:从 Durbin-Watson 统计量来看,存在显著的自相关性问题,需要进一步检验和修正。
自相关性可以考虑使用以下方法改进模型:使用自回归模型(AR模型),或者使用 Newey-West HAC 标准误来修正标准误差,使得结果更加稳健。这些方法可以帮助消除或降低自相关的影响,使模型更为稳健(Robust)。

总结

多元线性回归的核心在于通过多个自变量的线性组合来预测因变量的变化,但模型的有效性依赖于假设的满足。当假设未被满足时,估计的回归系数可能不再准确,并且显著性检验结果会失效。异方差性和多重共线性是多元线性回归中常见的两个问题,通过适当的检验和修正方法,可以提高模型的解释能力和预测效果。此外,研究者应结合实际问题,选择适合的模型和方法来处理回归分析中的复杂性。多元线性回归虽然是一种经典的统计方法,但其应用必须谨慎。通过合理的假设检验和对潜在问题的处理,可以确保回归结果的可靠性,并为决策和政策分析提供有力的依据。

参考资料

  1. 7.3 多元回归分析(multiple regression)
  2. 关于波士顿房价数据集的回归预测
posted @ 2024-10-14 15:21  郝hai  阅读(291)  评论(0编辑  收藏  举报