这篇博客，伤到我的脑细胞了——线性回归分析案例

最新推荐文章于 2023-05-07 19:38:49 发布
最新推荐文章于 2023-05-07 19:38:49 发布
阅读量1.1k 收藏 7
点赞数 2
分类专栏： 分析报告库
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_42521211/article/details/105725543
版权

目录

一、背景介绍

 在宏观经济学的理论中，投资是构成GDP的一部分，我国经济增长一般都用我国的GDP总量来衡量。从这一点来说，投资与经济增长之间应该存在着相关关系。由于投资与经济增长存在着相互作用，因此，从理论上讲固定资产投资对经济增长存在着相互作用、相互影响的关系。固定资产投资的增加可以推动经济总量的增长，同时经济增长的水平在一定程度上也决定着固定资产总量的水平。为了探讨其具体的定量关系，本文将根据北京市统计年鉴的2005年-2018年各行业固定资产投资额和GDP数据，进行回归分析，探讨行业固定资产投资额对行业GDP的贡献情况，并分析全市GDP数据与各行业投资额之间的量化关系。
在这里插入图片描述
数据来源：北京市统计局1

下述代码 阅读/运行 须知：
为了第二、第三及第四部分的内容展示上稍微精简一些，将部分可能重复用到的代码进行了打包工作，统一在第五部分展示，即 分析函数化Codes 部分，因此如需运行第二、第三及第四部分的代码，请先运行第五部分的代码。

二、数据探索

# 读取数据
sheet_names = ["北京市投资与gdp","各行业gdp","各行业固定资产投资"]
names = ["年份","全市","建筑业","批发和零售业","交通运输、仓储和邮政业","住宿和餐饮业" ,
         "信息传输、软件和信息技术服务业","金融业","房地产业","租赁和商务服务业",
         "科学研究和技术服务业","水利、环境和公共设施管理业","居民服务、修理和其他服务业",
         "教育","卫生和社会工作","文化、体育和娱乐业","公共管理、社会保障和社会组织"]

# data = pd.read_excel('20200306 各行业投资与经济增长回归(慕).xls',sheet_name=sheet_names[1],
#                      skiprows=3,names =names,index_col=0,dtypes={names[1]:np.float})
industry_gdp_data = pd.read_excel('20200306 各行业投资与经济增长回归(慕).xls',sheet_name=sheet_names[1],
                     header=0,index_col=0)
industry_fac_data = pd.read_excel('20200306 各行业投资与经济增长回归(慕).xls',sheet_name=sheet_names[2],
                     header=0,index_col=0)

可视化了解各行业GDP与固定资产投资的关系

# 画散点图
fig = plt.figure(figsize=(180,70))
for i in range(1,len(names)):
    ax = fig.add_subplot(2,8,i)
    ax.scatter(industry_fac_data[names[i]],industry_gdp_data[names[i]],s=500)
    plt.title(names[i],fontsize=56)
    plt.xlabel(names[i]+'投资额',fontsize=48)
    plt.ylabel(names[i]+'GDP',fontsize=48)
    plt.xticks(fontsize=48)
    plt.yticks(fontsize=48)
plt.show()
# fig.tight_layout() #自动调整格式

在这里插入图片描述

# 全市各行业的GDP与对应行业投资额的相关性
corr_gdp_industry_list = []
for i in range(1,len(names)):
    corr_gdp_industry_list.append(industry_fac_data[names[i]].corr(industry_gdp_data[names[i]]))

figure = plt.figure(figsize=(20,10))
corr_bar = plt.bar(names[1:len(names)],corr_gdp_industry_list)
plt.xticks(rotation=60,fontsize=12)
for bar in corr_bar:
    height = bar.get_height()
    plt.text(bar.get_x() + bar.get_width() / 2, height, str(('%.2f')%height), ha='center', va='bottom',fontsize=14)
plt.axhline(0.9,color='r')
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.title('全市各行业的GDP与对应行业投资额相关性条形图',fontsize=16)
plt.show()

在这里插入图片描述
小结：
通过上图各行业GDP与固定资产投资的关系散点图、相关性结果可知

  • 交通运输仓储和邮政业、信息传输软件和信息技术服务业、房地产业、水利环境和公共设施管理业 4个行业的投资额与其各自的GDP，具有较好的相关性，超过0.9，其他行业的gdp和投资额的相关性还需要进一步来判断。
  • 另外，全市的投资额和GDP，也具有很强的相关性。

箱型图，了解各行业GDP数据分布情况

title_name = '各行业GDP数据箱线图'
sub_ylabel = "GDP"
industry_boxplot_mu(industry_gdp_data,title_name,sub_ylabel)

在这里插入图片描述


plt.figure(figsize=(15,10))
industry_gdp_data.ix[:,1:].boxplot(vert=False,
                                   showmeans=True, # 以点的形式显示均值
                                  meanprops = {'marker':'D','markerfacecolor':'indianred','linewidth':1} # 设置均值点的属性，点的形状、填充色
                                  )

plt.xlabel('GDP',fontsize=14)
plt.ylabel('行业',fontsize=14)
plt.xticks(rotation=0,fontsize=12)
plt.yticks(fontsize=12)
plt.title('各行业GDP数据箱线图',fontsize=18)

在这里插入图片描述
小结：
通过各行业GDP的箱线图可知，各行业的GDP数据均无异常，但部分行业的GDP存在左偏或者右偏现象。

  • 批发和零售业、住宿和餐饮业、左偏现象比较明显；
  • 信息传输软件和信息服务业、金融业、科学研究和技术服务业、水利环境额 公共设施管理业、教育、卫生和社会工作，右偏现象较明显

箱型图，了解各行业投资额数据分布情况

title_name = '各行业 投资额数据箱线图'
sub_ylabel = "投资额"
industry_boxplot_mu(industry_fac_data,title_name,sub_ylabel)

在这里插入图片描述


plt.figure(figsize=(15,10))
industry_fac_data.ix[:,1:].boxplot(vert=False,
                                   showmeans=True, # 以点的形式显示均值
                                  meanprops = {'marker':'D','markerfacecolor':'indianred','linewidth':1} # 设置均值点的属性，点的形状、填充色
                                  )
plt.title('各行业投资额数据箱线图',fontsize=16)
plt.xlabel('GDP',fontsize=14)
plt.ylabel('行业',fontsize=14)
plt.xticks(rotation=0,fontsize=12)
plt.yticks(fontsize=12)

plt.show()

在这里插入图片描述
小结：
通过各行业投资额的箱线图可知，行业的投资额历年数据，也存在一定的偏态，并且部分行业的投资额变化较大。

  • 诸如建筑业、批发和零售业、交通运输仓储和邮政业、信息传输软件和信息技术服务业、租赁和商务服务业、文化体育和娱乐业，根据箱线图判断，存在异常数据；
  • 除过上述含异常值的行业，其余行业中，金融业、房地产业、科学研究和技术服务业、教育、卫生和社会工作，投资额存在左偏现象；
  • 而水利环境和公共设施管理业、居民服务修理和其他服务业、公共管理社会保障和社会组织，投资额存在右偏现象。
    可视化了解各行业GDP之间的关系
    （1）画矩阵散点图
'''
kdeplot(核密度估计图):核密度估计(kernel density estimation)是在概率论中用来估计未知的密度函数，
              属于非参数检验方法之一。通过核密度估计图可以比较直观的看出数据样本本身的分布特征。
distplot :displot()集合了matplotlib的hist()与核函数估计kdeplot的功能，
              增加了rugplot分布观测条显示与利用scipy库fit拟合参数分布的新颖用途。
'''
g = sns.PairGrid(industry_gdp_data)
g.map_diag(sns.distplot) #设置对角线图表
g.map_upper(plt.scatter, color = 'r')  #设置对角线上端图表
g.map_lower(sns.kdeplot, cmap='Blues_d') #设置对角线下端图表

在这里插入图片描述
（2） 画相关系数矩阵图

gdp_corr = industry_gdp_data.corr()  # gdp相关系数矩阵

vmax=gdp_corr.max().max()
vmin=gdp_corr.min().min()
print("最小相关系数：",vmin)

plt.figure(figsize=(15,15))

'''
annot=True，意思是显式热力图上的数值大小
square=True，意思是将图变成一个正方形，默认是一个矩形
cmap="Blues"是一种模式，就是图颜色配置方案
vmax是显示最大值
'''
sns.heatmap(gdp_corr,annot=True,vmax=vmax,vmin=vmin,square=True,cmap='Blues')

在这里插入图片描述
小结：
通过散点图和相关系数矩阵图，均可以发现，各行业的GDP之间具有较强的正相关，普遍相关系数，均在0.8以上。

可视化了解各行业固定资产投资的关系
（1） 画矩阵散点图

'''
kdeplot(核密度估计图):核密度估计(kernel density estimation)是在概率论中用来估计未知的密度函数，
              属于非参数检验方法之一。通过核密度估计图可以比较直观的看出数据样本本身的分布特征。
distplot :displot()集合了matplotlib的hist()与核函数估计kdeplot的功能，
              增加了rugplot分布观测条显示与利用scipy库fit拟合参数分布的新颖用途。
'''
g = sns.PairGrid(industry_fac_data)
g.map_diag(sns.distplot) #设置对角线图表
g.map_upper(plt.scatter, color = 'r')  #设置对角线上端图表
g.map_lower(sns.kdeplot, cmap='Blues_d') #设置对角线下端图表

在这里插入图片描述
（2）画相关系数矩阵图

fac_corr = industry_fac_data.corr()  # gdp相关系数矩阵

vmax=fac_corr.max().max()
vmin=fac_corr.min().min()
print("最小相关系数：",vmin)
# print("最大相关系数：",vmax)


plt.figure(figsize=(15,15))
sns.heatmap(fac_corr,annot=True,vmax=vmax,vmin=vmin,square=True,cmap='Blues')
plt.title("各行业固定资产投资的关系相关系数矩阵图",fontsize=18)

在这里插入图片描述
小结：
对于各行业的投资额，通过散点图矩阵和相关系数矩阵可以发现，

  • 行业间的投资额的相关性，相对行业GDP的相关性要弱一些，并且部分行业间投资额具有负相关性。
  • 建筑业、文化体育和娱乐业、公共管理社会保障和社会组织，相对来说，普遍与其他行业的相关性较弱。

通过上述分析，信息传输软件和信息技术服务业、房地产业、水利环境和公共设施管理业 3个行业的投资额与其各自的GDP，具有较好的线性关系；大部分行业的历年GDP和投资额，存在左偏或者右偏现象；行业GDP之间相关性普遍较强，而行业投资额间相关性则分化比较大。

三、信息传输、软件和信息技术服务业回归分析

1. 投资额和GDP可视化

industry_name = "信息传输、软件和信息技术服务业" 
ytrain = industry_gdp_data[industry_name]
xtrain = industry_fac_data[industry_name]

2. 回归分析

（1）模型一：原始数据

model,y_pre = simple_regression_mu(ytrain,xtrain)
model.summary()

在这里插入图片描述
在这里插入图片描述
小结
模型通过检验，模型的系数通过检验，但是模型的常数项没有通过检验。另外，D.W.的值为2.152，根据DW检验表，DW取值在（dU, 4- dU）之间，接受原假设H0 ,认为ut 非自相关，即在信息传输、软件和信息技术服务业，残差非自相关，但也有可能残差具有较高的自相关性。

# 拟合结果可视化
visual_regression_result(range(len(ytrain)),ytrain,y_pre)
# visual_regression_result(model,sm.add_constant(xtrain))

在这里插入图片描述
** 残差分析**

线性模型的残差通常具有正太分布。我们可以绘制残差密度来检查正态性。


fig = plt.figure(figsize=(15,5))
ax = fig.add_subplot(121)
# model.resid.plot.density()  #残差密度图
model.resid.hist(bins=5,alpha = 0.5)
model.resid.plot(kind = 'kde', secondary_y=True)
ax.set_title('密度直方图')
# plt.grid()

ax2 = fig.add_subplot(122)
# qq图
stats.probplot(model.resid, dist="norm", plot=plt)
plt.show()

在这里插入图片描述

# 残差正态性检验

from scipy.stats import kstest    # K-S检验：特点是比较严格，基于的原理是CDF，理论上可以检验任何分布
ks_result = kstest(model.resid,'norm')
print('kstest检验：',ks_result)
 
from scipy.stats import shapiro   # Shapiro检验：专门用来检验正态分布  shapiro 不适合做样本数＞5000的正态性检验，检验结果的P值可能不准确
shap_result = shapiro(model.resid)
print('shapiro检验：',shap_result)

from scipy.stats import normaltest  # Normal检验：原理是基于数据的skewness和kurtosis
norm_result = normaltest(model.resid,axis=None)
print('normaltest检验：',norm_result)

kstest检验2： KstestResult(statistic=0.5714285714285714, pvalue=7.456196123263218e-05)
shapiro检验： (0.901134729385376, 0.11714380979537964)
normaltest检验： NormaltestResult(statistic=6.129670593453145, pvalue=0.04666152653777272)

# 定量判断 残差同方差检验
# White's Test
white_test = sm.stats.diagnostic.het_white(model.resid, exog = model.model.exog)

print('模型方差齐性检验统计量：%.2f，P值：%.4f，Fvalue：%.2f，F_pvalue：%.2f：' %(white_test))
# model.model.exog

# 定性判断  残差同方差检验   及残差的自相关性判断
# 判断标准：如果残差在参考线两侧均匀分布，则意味着异方差性较弱；而如果呈现出明显的不均匀分布，则意味着存在明显的异方差性。
residual_scatter_mu(model)

模型方差齐性检验统计量：1.95，P值：0.3779，Fvalue：0.89，F_pvalue：0.44
在这里插入图片描述
小结：
根据怀特检验，残差通过方差齐性检验，但是通过残差序列图可以发现，残差还是存在一定的异方差情况；另外，根据残差序列图和残差的散点图，可以发现残差具有一定的自相关性

模型评价

print('评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：%.3f'
      %(mean_absolute_per_err(y_pre,ytrain)))

长度： 14
评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：0.134

（2） 模型二：消除异方差

#  通过对变量取对数消除异方差
log_ytrain = np.log(ytrain)
log_xtrain = np.log(xtrain)

gdp_fac_graph_mu(log_ytrain,log_xtrain)

在这里插入图片描述

model2 ,log_y_pre= simple_regression_mu(log_ytrain,log_xtrain)
model2.summary()

在这里插入图片描述

# 拟合结果可视化
visual_regression_result(range(len(log_ytrain)),ytrain,np.exp(log_y_pre))

在这里插入图片描述
残差分析

fig = plt.figure(figsize=(15,5))
ax = fig.add_subplot(121)
# model.resid.plot.density()  #残差密度图
model2.resid.hist(bins=5,alpha = 0.5)
model2.resid.plot(kind = 'kde', secondary_y=True)
ax.set_title('密度直方图')
# plt.grid()

ax2 = fig.add_subplot(122)
# qq图
stats.probplot(model2.resid, dist="norm", plot=plt)
plt.show()

在这里插入图片描述

# 残差正态性检验

from scipy.stats import kstest    # K-S检验：特点是比较严格，基于的原理是CDF，理论上可以检验任何分布
ks_result = kstest(model2.resid,'norm')
print('kstest检验：',ks_result)
 
from scipy.stats import shapiro   # Shapiro检验：专门用来检验正态分布  shapiro 不适合做样本数＞5000的正态性检验，检验结果的P值可能不准确
shap_result = shapiro(model2.resid)
print('shapiro检验：',shap_result)

from scipy.stats import normaltest  # Normal检验：原理是基于数据的skewness和kurtosis
norm_result = normaltest(model2.resid,axis=None)
print('normaltest检验：',norm_result)

kstest检验： KstestResult(statistic=0.403305160372756, pvalue=0.014364366464748599)
shapiro检验： (0.931269645690918, 0.3178831934928894)
normaltest检验： NormaltestResult(statistic=1.9191948456802586, pvalue=0.3830470609372596)

# 定量判断 残差同方差检验(残差和自变量之间)
# White's Test
white_test = sm.stats.diagnostic.het_white(model2.resid, exog = model2.model.exog)

print('模型方差齐性检验统计量：%.2f，P值：%.4f，Fvalue：%.2f，F_pvalue：%.2f：' %(white_test))
# model.model.exog

# 定性判断  残差同方差检验   及残差的自相关性判断
# 判断标准：如果残差在参考线两侧均匀分布，则意味着异方差性较弱；而如果呈现出明显的不均匀分布，则意味着存在明显的异方差性。
residual_scatter_mu(model2)

模型方差齐性检验统计量：1.46，P值：0.4829，Fvalue：0.64，F_pvalue：0.55
在这里插入图片描述
模型评价

print('评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：%.3f'
      %(mean_absolute_per_err(log_y_pre,log_ytrain)))

re_log_y_pre = np.exp(log_y_pre)
print('还原数据后，评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：%.3f'
      %(mean_absolute_per_err(re_log_y_pre,ytrain)))

长度： 14
评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：0.018
长度： 14
还原数据后，评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：0.130

（3）模型三：模型变换 加入一阶滞后项，消除自相关


test = industry_fac_data[industry_name][1:]
test2 = pd.Series(industry_fac_data[industry_name][:-1])
test2.index = test.index

add_xtrain = pd.concat([test,test2],axis=1)
add_xtrain.columns = [industry_name,'滞后一年']

add_ytrain = industry_gdp_data[industry_name][1:]


# industry_FAC = sm.add_constant(industry_FAC) # 线性回归增加常数项 y=kx+b

# model = sm.OLS(industry_GDP, industry_FAC) # 普通最小二乘模型，ordinary least square model
# model = model.fit()    #res.model.endog
# # 模型检验及系数检验
# print(industry_name+" 的回归模型检验及系数检验结果：\n")
# model.summary()

model3 ,add_y_pre= simple_regression_mu(add_ytrain,add_xtrain)
model3.summary()

在这里插入图片描述

# 拟合结果可视化
visual_regression_result(range(len(add_ytrain)),add_ytrain,add_y_pre)

在这里插入图片描述
残差分析

# plt.figure()
# model.resid.plot.density()  #残差密度图
# model.resid.hist(bins=5,alpha = 0.5)
# model.resid.plot(kind = 'kde', secondary_y=True)

fig = plt.figure(figsize=(15,5))
ax = fig.add_subplot(121)
# model.resid.plot.density()  #残差密度图
model3.resid.hist(bins=5,alpha = 0.5)
model3.resid.plot(kind = 'kde', secondary_y=True)
ax.set_title('密度直方图')
# plt.grid()

ax2 = fig.add_subplot(122)
# qq图
stats.probplot(model3.resid, dist="norm", plot=plt)
plt.show()

在这里插入图片描述

# 残差正态性检验

from scipy.stats import kstest    # K-S检验：特点是比较严格，基于的原理是CDF，理论上可以检验任何分布
ks_result = kstest(model3.resid,'norm')
print('kstest检验：',ks_result)
 
from scipy.stats import shapiro   # Shapiro检验：专门用来检验正态分布
shap_result = shapiro(model3.resid)
print('shapiro检验：',shap_result)

from scipy.stats import normaltest  # Normal检验：原理是基于数据的skewness和kurtosis
norm_result = normaltest(model3.resid)
print('normaltest检验：',norm_result)

kstest检验： KstestResult(statistic=0.5384615384615384, pvalue=0.0004708705798459567)
shapiro检验： (0.9762569069862366, 0.9561669230461121)
normaltest检验： NormaltestResult(statistic=0.21037569966411002, pvalue=0.9001554126600356)

# 定性判断  残差同方差检验   及残差的自相关性判断
# 判断标准：如果残差在参考线两侧均匀分布，则意味着异方差性较弱；而如果呈现出明显的不均匀分布，则意味着存在明显的异方差性。
residual_scatter_mu(model3)

# 定量判断 残差的同方差检验
# White's Test
white_test = sm.stats.diagnostic.het_white(model3.resid, exog = model3.model.exog)

print('模型方差齐性检验统计量：%.2f，P值：%.4f，Fvalue：%.2f，F_pvalue：%.2f：' %(white_test))
# model.model.exog

模型方差齐性检验统计量：5.58，P值：0.3498，Fvalue：1.05，F_pvalue：0.46
在这里插入图片描述
模型评价

print('评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：%.3f'
      %(mean_absolute_per_err(add_y_pre,add_ytrain)))

长度： 13
评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：0.093

（4）模型四 加入一阶滞后项和数据取对数

  • 模型变换 加入一阶滞后项，消除自相关；
  • 取对数，消除异方差
log_add_xtrain = np.log(add_xtrain)
log_add_ytrain = np.log(add_ytrain)

model4 ,log_add_y_pre= simple_regression_mu(log_add_ytrain,log_add_xtrain)
model4.summary()

在这里插入图片描述


# 拟合结果可视化
visual_regression_result(range(len(log_add_ytrain)),log_add_ytrain,log_add_y_pre)

在这里插入图片描述
残差分析

fig = plt.figure(figsize=(15,5))
ax = fig.add_subplot(121)
# model.resid.plot.density()  #残差密度图
model4.resid.hist(bins=5,alpha = 0.5)
model4.resid.plot(kind = 'kde', secondary_y=True)
ax.set_title('密度直方图')
# plt.grid()

ax2 = fig.add_subplot(122)
# qq图
stats.probplot(model4.resid, dist="norm", plot=plt)
plt.show()

在这里插入图片描述

# 残差正态性检验

from scipy.stats import kstest    # K-S检验：特点是比较严格，基于的原理是CDF，理论上可以检验任何分布
ks_result = kstest(model4.resid,'norm')
print('kstest检验：',ks_result)
 
from scipy.stats import shapiro   # Shapiro检验：专门用来检验正态分布
shap_result = shapiro(model4.resid)
print('shapiro检验：',shap_result)

from scipy.stats import normaltest  # Normal检验：原理是基于数据的skewness和kurtosis
norm_result = normaltest(model4.resid)
print('normaltest检验：',norm_result)

kstest检验： KstestResult(statistic=0.4434958990037281, pvalue=0.007560196567781911)
shapiro检验： (0.9548967480659485, 0.6739445924758911)
normaltest检验： NormaltestResult(statistic=0.5897760812235786, pvalue=0.744614949433436)

# 定性判断  残差同方差检验   及残差的自相关性判断
# 判断标准：如果残差在参考线两侧均匀分布，则意味着异方差性较弱；而如果呈现出明显的不均匀分布，则意味着存在明显的异方差性。
residual_scatter_mu(model4)

# 定量判断 残差的同方差检验
# White's Test
white_test = sm.stats.diagnostic.het_white(model4.resid, exog = model4.model.exog)

print('模型方差齐性检验统计量：%.2f，P值：%.4f，Fvalue：%.2f，F_pvalue：%.2f：' %(white_test))
# model.model.exog

模型方差齐性检验统计量：3.26，P值：0.6606，Fvalue：0.47，F_pvalue：0.79
在这里插入图片描述

print('评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：%.3f'
      %(mean_absolute_per_err(log_add_y_pre,log_add_ytrain)))

re_log_add_y_pre = np.exp(log_add_y_pre)
print('还原数据后，评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：%.3f'
      %(mean_absolute_per_err(re_log_add_y_pre,add_ytrain)))

长度： 13
评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：0.012
长度： 13
还原数据后，评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：0.090

（5）单位根、协整、格兰杰因果关系

# 单位根检验
# 单位根检验是指检验序列中是否存在单位根，因为存在单位根就是非平稳时间序列了。
# 单位根就是指单位根过程，可以证明，序列中存在单位根过程就不平稳，会使回归分析中存在伪回归

x = diff_adf_test(log_xtrain,zhixindu='5%')
y = diff_adf_test(log_ytrain,zhixindu='5%')

序列1阶平稳，统计量为：-4.2812 ,p值为0.0005, 5%临界值为：-3.3672
序列1阶平稳，统计量为：-4.4592 ,p值为0.0002, 5%临界值为：-3.1894

# 协整检验
print(coint(y,x))
print(coint(x,y))

(-3.499788711567408, 0.03241412320391135, array([-5.04192472, -3.89268694, -3.41677222]))
(-0.141916031919859, 0.9820000262274053, array([-5.04192472, -3.89268694, -3.41677222]))

# 格兰杰因果关系检验
grangercausalitytests(pd.DataFrame([x,y]).T,2)

在这里插入图片描述

# 格兰杰因果关系检验
grangercausalitytests(pd.DataFrame([y,x]).T,2)

在这里插入图片描述

四、多元回归分析

4.1 可视化北京市GDP与各行业投资额的关系

fig = plt.figure(figsize=(60,50))
for i in range(2,len(names)):
    ax = fig.add_subplot(3,5,i-1)
    ax.scatter(industry_fac_data[names[i]],industry_gdp_data[names[1]],s=500)
#     plt.title(names[i],fontsize=56)
    plt.xlabel(names[i]+'投资额',fontsize=48)
    plt.ylabel(names[1]+'GDP',fontsize=48)
    plt.xticks(fontsize=48)
    plt.yticks(fontsize=48)
plt.show()

在这里插入图片描述

# 北京市GDP与各行业投资额的相关系数
corr_gdp_industry_list = []
for i in range(1,len(names)):
    corr_gdp_industry_list.append(industry_gdp_data[names[1]].corr(industry_fac_data[names[i]]))
# corr_series = pd.Series(corr_gdp_industry_list,index=names[2:len(names)])

figure = plt.figure(figsize=(15,5))
corr_bar = plt.bar(names[1:len(names)],corr_gdp_industry_list)
plt.xticks(rotation=60,fontsize=12)
for bar in corr_bar:
    height = bar.get_height()
    plt.text(bar.get_x() + bar.get_width() / 2, height, str(('%.2f')%height), ha='center', va='bottom',fontsize=14)
plt.axhline(0.8,color='r')
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.title('全市GDP与各行业投资额相关性条形图',fontsize=16)
plt.show()

在这里插入图片描述
小结：
通过上述北京市GDP和各行业的投资额散点图可知，交通运输仓储和邮政业、信息传输软件和信息技术服务业、房地产业、水利环境和公共设施管理业、教育、卫生和社会工作等6个行业，具有较强的相关性，达到0.8以上。

4.2 多元回归分析全市GDP和行业投资的关系

筛选上述相关性较强的六个行业与全市GDP关系进行回归分析。

# 
flag = pd.Series(corr_gdp_industry_list)>0.8  
x_name = [names[i] for i in range(2,len(names)) if flag[i-1]]
X = industry_fac_data[x_name]
X = sm.add_constant(X)
Y = industry_gdp_data[names[1]]
model = sm.OLS(Y,X)
model = model.fit()
model.summary()
# industry_gdp_data.index.names

在这里插入图片描述
根据上述结果，剔除交通运输仓储和邮政业、房地产业、卫生和社会工作，重新进行拟合。

x_name = ["信息传输、软件和信息技术服务业","水利、环境和公共设施管理业","教育"]

X = industry_fac_data[x_name]
X = sm.add_constant(X)
Y = industry_gdp_data[names[1]]
model = sm.OLS(Y,X)
model = model.fit()
model.summary()

在这里插入图片描述
对比未剔除交通运输仓储和邮政业、房地产业、卫生和社会工作的回归模型，剔除后，模型效果更好。

Y_p = model.predict(X)
# 拟合结果可视化
visual_regression_result(range(len(Y)),Y,Y_p)
# visual_regression_result(model,sm.add_constant(xtrain))

在这里插入图片描述

fig = plt.figure(figsize=(15,5))
ax = fig.add_subplot(121)
# model.resid.plot.density()  #残差密度图
model.resid.hist(bins=5,alpha = 0.5)
model.resid.plot(kind = 'kde', secondary_y=True)
ax.set_title('密度直方图')
# plt.grid()

ax2 = fig.add_subplot(122)
# qq图
stats.probplot(model.resid, dist="norm", plot=plt)
plt.show()

在这里插入图片描述

# 残差正态性检验

from scipy.stats import kstest    # K-S检验：特点是比较严格，基于的原理是CDF，理论上可以检验任何分布
ks_result = kstest(model.resid,'norm')
print('kstest检验：',ks_result)
 
from scipy.stats import shapiro   # Shapiro检验：专门用来检验正态分布  shapiro 不适合做样本数＞5000的正态性检验，检验结果的P值可能不准确
shap_result = shapiro(model.resid)
print('shapiro检验：',shap_result)

from scipy.stats import normaltest  # Normal检验：原理是基于数据的skewness和kurtosis
norm_result = normaltest(model.resid,axis=None)
print('normaltest检验：',norm_result)

kstest检验： KstestResult(statistic=0.5714285714285714, pvalue=7.456196123263218e-05)
shapiro检验： (0.9852579236030579, 0.9947238564491272)
normaltest检验： NormaltestResult(statistic=0.3540289636912072, pvalue=0.8377676520133555)

# 定量判断 残差同方差检验
# White's Test
white_test = sm.stats.diagnostic.het_white(model.resid, exog = model.model.exog)

print('模型方差齐性检验统计量：%.2f，P值：%.4f，Fvalue：%.2f，F_pvalue：%.2f：' %(white_test))
# model.model.exog

模型方差齐性检验统计量：6.77，P值：0.6608，Fvalue：0.42，F_pvalue：0.87

print('评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：%.3f'
      %(mean_absolute_per_err(Y_p,Y)))

长度： 14
评价模型标准-相对误差绝对值平均MAPE (Mean Absolute Percentage Error)：0.031

五、分析函数化Codes

import pandas as pd
import numpy as np
import scipy
import  matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
import statsmodels.api as sm # 最小二乘
from statsmodels.stats.outliers_influence import summary_table # 获得汇总信息
import statsmodels.formula.api as smf

from scipy import stats
from sklearn.metrics import mean_absolute_error

from statsmodels.tsa.stattools import adfuller   # ADF检验
from statsmodels.tsa.stattools import coint     # 协整检验
from statsmodels.tsa.stattools import grangercausalitytests  # 格兰杰因果关系检验

import warnings
warnings.filterwarnings("ignore")# 忽略静态警告

mpl.style.use('seaborn')

# 为了画图中文可以正常显示
plt.rcParams['font.sans-serif'] = ['simhei'] #指定默认字体
plt.rcParams['axes.unicode_minus'] = False  #解决保存图像时负号'-'显示为方块的问题

###########   函数定义  区域  ###############################################

def industry_boxplot_mu(data,title_name,sub_ylabel):
    '''
    title_name: boxplot的title 
    '''
    fig = plt.figure(figsize=(100,60))
    plt.title(title_name,fontsize=56)
    for i in range(1,len(names)):
        ax = fig.add_subplot(2,8,i)
#         ax.boxplot(industry_gdp_data[names[i]],
        ax.boxplot(data[names[i]],
                   patch_artist=True, # 要求用自定义颜色填充盒形图，默认白色填充
                   showmeans=True, # 以点的形式显示均值
                   boxprops = {'color':'black','facecolor':'#ffffff','linewidth':3}, # 设置箱体属性，填充色和边框色
                   flierprops = {'marker':'o','markerfacecolor':'red','color':'black','linestyle':'-','linewidth':3}, # 设置异常值属性，点的形状、填充色和边框色
                   meanprops = {'marker':'D','markerfacecolor':'indianred','linewidth':15}, # 设置均值点的属性，点的形状、填充色
                   medianprops = {'linestyle':'-','color':'red','linewidth':3}) # 设置中位数线的属性，线的类型和颜色
        # plt.title(names[i],fontsize=56)
        plt.xlabel(names[i],fontsize=48)
        # plt.ylabel(names[i]+'GDP',fontsize=48)
        plt.ylabel(names[i]+sub_ylabel,fontsize=48)
        plt.xticks(fontsize=48)
        plt.yticks(fontsize=48)
    plt.show()

'''
可视化行业的历年GDP、历年投资额，以及投资额和GDP的散点图
'''
def gdp_fac_graph_mu(ytrain,xtrain):
    '''
    ytrain：行业GDP数据
    xtrain：行业投资额数据
    '''
#     industry_GDP = industry_gdp_data[industry_name]
#     industry_FAC = industry_fac_data[industry_name]

    # 画图
    figure = plt.figure(figsize=(20,5))
    ax = figure.add_subplot(121)
    ax.plot(xtrain.index,xtrain,'rp-',label='投资额')
    ax.set_ylabel('投资额')
    ax.set_xlabel('年份')
    ax.legend(loc=0)
    
    ax = ax.twinx()
    ax.plot(ytrain.index,ytrain,'*-',label='GDP')
    ax.set_ylabel('GDP')
    # ax.set_xlabel('年份')
    ax.legend(loc=1)
    ax.set_title('时间序列图',fontsize=14)
    plt.xticks(fontsize=12)
    plt.yticks(fontsize=12)


    ax2 = figure.add_subplot(122)
    ax2.scatter(xtrain,ytrain,s=50)
    plt.title('投资额与GDP关系散点图',fontsize=14)
    plt.xlabel('投资额',fontsize=12)
    plt.ylabel('GDP',fontsize=12)
    plt.xticks(fontsize=12)
    plt.yticks(fontsize=12)
    
    plt.show()
    
'''
对各行业的投资额和GDP，进行一元回归建模
'''
def simple_regression_mu(ytrain,xtrain):
    '''
    params
    ytrain：行业GDP数据
    xtrain：行业投资额数据
    
    return 
    model: 拟合的回归模型
    '''
#     industry_GDP = industry_gdp_data[industry_name]
#     industry_FAC = industry_fac_data[industry_name]
    
    xtrain = sm.add_constant(xtrain) # 线性回归增加常数项 y=kx+b

    model = sm.OLS(ytrain, xtrain) # 普通最小二乘模型，ordinary least square model
    model = model.fit()    #res.model.endog
    # 模型检验及系数检验
#     print(industry_name+" 的回归模型检验及系数检验结果：\n",model.summary())
    y_pre = model.predict(xtrain)
    return model,y_pre

# '''
# 对各行业的投资额和GDP取自然对数后，进行一元回归建模
# '''
# def simple_regression_log_mu(ytrain,xtrain):
#     '''
#     params
#     ytrain：行业GDP数据
#     xtrain：行业投资额数据
    
#     return 
#     model: 拟合的回归模型
#     '''
#     industry_GDP = np.log(industry_gdp_data[industry_name])
#     industry_FAC = np.log(industry_fac_data[industry_name])
    
#     industry_FAC = sm.add_constant(industry_FAC) # 线性回归增加常数项 y=kx+b

#     model = sm.OLS(industry_GDP, industry_FAC) # 普通最小二乘模型，ordinary least square model
#     model = model.fit()    #res.model.endog
#     # 模型检验及系数检验
# #     print(industry_name+" 取自然对数后，回归模型检验及系数检验结果：\n",model.summary())
#     return model,industry_FAC

'''
回归结果可视化
'''
def visual_regression_result(xtrain,ytrain,y_pre):
    '''
    xtrain: 训练数据x
    ytrain:训练数据y
    y_pre:训练预测值
    '''
#     ytest =model.predict(xtrain)

    #绘制拟合结果和误差图
    figure = plt.figure(figsize=(10,5))
    plt.scatter(xtrain,ytrain,marker = 'd',color = 'k',label='原始值')
    plt.scatter(xtrain,y_pre,marker = 'p',color = 'g',label='预测值')

    plt.plot(xtrain,y_pre,color = 'r',label='预测线')
    plt.plot([xtrain,xtrain],[ytrain,y_pre] ,color = 'grey')
    plt.legend(fontsize = 14)
    plt.xticks(fontsize=12)
    plt.yticks(fontsize=12)
    plt.xlabel("序列",fontsize=12)
    plt.ylabel(industry_name+"GDP",fontsize=12)
    plt.title(industry_name+'线性回归拟合',fontsize=16)
    plt.show()

    
#     '''
# 回归结果可视化
# '''
# def visual_regression_result(xtrain,ytrain,y_pre):
#     '''
#     xtrain: 训练数据x
#     ytrain:训练数据y
#     y_pre:训练预测值
#     '''
# #     ytest =model.predict(xtrain)

#     #绘制拟合结果和误差图
#     figure = plt.figure(figsize=(12,7))
#     plt.scatter(xtrain,ytrain,marker = 'd',color = 'k',label='原始值')
#     plt.scatter(xtrain,y_pre,marker = 'p',color = 'g',label='预测值')

#     plt.plot(xtrain,y_pre,color = 'r',label='预测线')
#     plt.plot([xtrain,xtrain],[ytrain,y_pre],
#              color = 'grey')
#     plt.legend(fontsize = 14)
#     plt.xticks(fontsize=12)
#     plt.yticks(fontsize=12)
#     plt.xlabel(industry_name+"投资额",fontsize=12)
#     plt.ylabel(industry_name+"GDP",fontsize=12)
#     plt.title(industry_name+'线性回归拟合',fontsize=16)
#     plt.show()
'''
可视化回归后的残差散点图
'''
# def residual_scatter_mu(industry_name,model):
def residual_scatter_mu(model):
    '''
    industry_name：str,行业名称
    model: 拟合的模型
    '''
    # 残差均值  0
    err_mean = np.mean(model.resid)
    # 残差标准差
    err_std = np.std(model.resid)

    # 残差图
    figure = plt.figure(figsize=(20,5))
    ax = figure.add_subplot(121)
    plt.scatter(model.resid.index,model.resid,marker='d',color='k')
    plt.plot(model.resid.index,model.resid)
    plt.axhline(err_std)
    plt.axhline(-err_std)
    plt.axhline(0,color='r')
    
    plt.xlabel("年份",fontsize=12)
    plt.ylabel("残差",fontsize=12)
    plt.title(industry_name+"残差序列图",fontsize=14)
    
    ax2 = figure.add_subplot(122)
    plt.scatter(model.resid[:len(model.resid)-1],model.resid[1:len(model.resid)])
    plt.xlabel("滞后残差",fontsize=12)
    plt.ylabel("残差",fontsize=12)
    plt.title("一阶滞后变量的散点图",fontsize=14)
#     ax2 = figure.add_subplot(122)
#     plt.scatter(industry_gdp_data[industry_name],model.resid)
#     plt.axhline(err_std)
#     plt.axhline(-err_std)
#     plt.xlabel(industry_name+"GDP",fontsize=12)
#     plt.ylabel("残差",fontsize=12)
#     plt.title(industry_name+"残差散点图",fontsize=14)
    plt.show()


'''
相对误差绝对值平均
'''
def mean_absolute_per_err(y_pre,y): 
    '''
    y_pre:预测值
    y: 原始数据
    '''
    result = 0
    for i in range(len(y_pre)):
        result = result + np.abs((y_pre[i]-y[i])/y[i])
    result = result/len(y_pre)
    print('长度：',len(y_pre))
    return result

'''adfuller()返回值
df : float
    Test statistic
pvalue : float
    MacKinnon's approximate p-value based on MacKinnon (1994, 2010)
usedlag : int
    Number of lags used
nobs : int
    Number of observations used for the ADF regression and calculation of
    the critical values
critical values : dict
    Critical values for the test statistic at the 1 %, 5 %, and 10 %
    levels. Based on MacKinnon (2010)
'''
def diff_adf_test(data,max_jieshu=3,zhixindu='5%'):
    '''
    data:检验序列
    max_jieshu:差分的最高阶数
    '''
    adf_result = adfuller(data)
    if adf_result[0] < adf_result[4]['10%']:  # 平稳性判断
        print('原序列平稳，统计量为：%.4f ,p值为%.4f, %s临界值为：%.4f'%(adf_result[0],adf_result[1],zhixindu,adf_result[4][zhixindu]))
        return data
    for i in range(1,max_jieshu+1):
        temp = data.diff(i)[i:]
        adf_result = adfuller(temp)
        if adf_result[0] < adf_result[4]['10%']:  # 平稳性判断
            print('序列%d阶平稳，统计量为：%.4f ,p值为%.4f, %s临界值为：%.4f'%(i,adf_result[0],adf_result[1],zhixindu,adf_result[4][zhixindu]))
            return temp
    print('原序列在差分%d阶以内，均不平稳'%max_jieshu)
    
########################################################################

六、题外话

这部分的分析工作用了很长时间，也查看过一些资料，当真正准备把它当一个分析报告想展示的时候，自己突然发现不知道该展示什么，或者说要来说明什么问题。这也是自己第一次尝试分享展示这种分析案例，有很多粗糙的地方和不足，多多练手，多多参考大家的作品，提升自己的实力，并希望能够不断的分享新的有帮助的内容给阅读者。

  1. 北京市统计局 ↩︎

  2. Python做假设检验 ↩︎

积跬步，慕至千里
关注
  • 2
    点赞
  • 踩
  • 7
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 1
    评论
专栏目录
回归分析案例详解
09-11
根据下表某猪场25头肥猪的数据资料，瘦肉量y对其它变量进行多元回归分析。
Excel 回归分析实例
09-14
Excel回归分析应用实例，包括回归分析的四种方法以及一个实例，实例运用了三种模型来对比。本例简单明了，希望对想要用Excel做回归分析的朋友有点小小的帮助，本例为菜鸟级的介绍。
回归分析及应用实例
05-02
回归分析原理和实例举证，详细的回归分析过程。
案例：回归分析-R实现
R语言中文社区
05-27 5910
作者 | 结实来源 | DATA江湖1）业务需求：某牙膏制造企业为了更好的拓展产品市场，有效管理库存，董事会要求销售部门根据市场调查，找出公司生产牙膏销售量与销售价格以及...
机器学习案例--回归分析
xiayto的博客
03-31 1649
目录 1、最小二乘法 1.1 原理推导 1.2 代码样例 2、SKlearn 2.1 代码样例 2.2 多项式扩展 2.3 正则化 2.4 逻辑回归和多酚类 1、最小二乘法解 1.1 算法原理 回归算法推导 根据中心极限定理，误差服从正态分布，将误差带入之后用极大似然估计，取对数得到目标函数，目标函数是最小二乘的形式，求导得出的解是： θ=(XTX)−1XTYθ=(XTX)−1...
机器学习——线性回归介绍及案例实战的数据集（客户价值数据表）
04-06
本资源适用于多元线性回归的Python代码实现，是小编的机器学习——线性回归介绍及案例实战的数据集
关于多元线性回归分析——Python&SPSS
12-20
在这个案例中，我们关注的是如何使用Python和SPSS进行多元线性回归分析。数据集`Folds5x2_pp.csv`包含了五个变量：AT（温度），V（压力），AP（湿度），RH（压强）和PE（输出电力）。目标是找到一个线性模型，该模型...
数学建模——线性回归分析PPT教案.pptx
10-05
"数学建模——线性回归分析PPT教案" 本资源是关于数学建模和线性回归分析的PPT教案，共81页，涵盖了电力市场的输电阻塞管理、线性回归分析等多个方面的内容。本教案适合高等教育和研究机构的学生和教师使用。 从...
MATLAB线性回归分析案例
最新发布
06-07
这个资源是一个使用MATLAB进行线性回归分析的教学案例，它提供了一个实际应用场景，通过分析广告支出与产品销售量之间的关系，演示了如何利用MATLAB软件进行数据处理、可视化以及建模预测。以下是对资源的简介： - ...
18.线性回归及多项式回归分析四个案例分享1
08-03
在本文中，我们探讨了四个使用Python进行线性回归分析的案例，涉及到数据预处理、模型训练、预测以及结果可视化。 1. **线性回归预测Pizza价格案例**： - 这个案例中，我们利用一元线性回归模型预测不同直径的...
回归分析原理及例子
12-13
回归分析是一种处理变量的统计相关关系的一种数理统计方法。回归分析的基本思想是: 虽然自变量和因变量之间没有严格的、确定性的函数关系, 但可以设法找出最能代表它们之间关系的数学表达形式
R完整回归案例，以二手车数据为例
07-22
以狗熊会二手车数据为例，利用R语言进行回归分析，涉及原始数据的清洗、描述统计分析、经典回归模型建立、异常点诊断、交叉验证等技术，附有原始数据以及R完整代码、注释。
回归分析案例分析全流程
m0_37228052的博客
05-31 2774
一、案例说明 1.案例数据 在“工资影响因素”的调查问卷中，调查了每个人的起始工资、工作经验、受教育年限、受雇月数、职位等级以及当前工资六个方面。 2.分析目的 目的是建立以当前工资为因变量的回归模型，并得出结论。[案例来源于：SPSS统计分析（第5版）卢纹岱,朱红兵主编，案例有一些变动 具体请看分析。] 二、数据清理 在数据分析之前，首先需要进行数据查看，包括数据中是否有异常值，无效样本等。如果有异常值则需要进行处理，然后再进行分析。另外如果数据中有无效样本也需要进行处理后再进行分析。无效..
数据分析: 线性回归案例
logan_li2010的博客
12-22 1431
import numpy as np import pandas as pd import matplotlib.pyplot as plt %matplotlib inline from sklearn.datasets import load_boston boston = load_boston() data = boston['data'] target = boston['targ...
[机器学习]回归分析案例 --- 网游排名分析
stary_yan的博客
07-28 1272
[机器学习]回归分析案例 — 网游排名分析项目描述案例来自狗熊会。 基础案例 | 玩家必看之网游排行big picture从项目描述以及数据上看，可以很明显地发现，这是一个非常文本化的数据，也就是无法直接使用机器学习来解决问题。那么，这就非常依赖于将文本量化的能力。本文，我会尝试着完整地把整个项目实现一遍，并给出我的一些理解和感悟。Research首先，我们需要对文本数据有一个大体的认识。各字段的
【SPSS】回归分析详细操作教程（附案例实战）
m0_64336780的博客
05-07 3万+
本文主要讲解使用SPSS进行回归分析，包括线性回归、曲线估计、二元Logistic回归、多元Logistic回归分析，使用案例来讲解具体操作过程
gamma分布_负二项分布及其应用
weixin_39853210的博客
11-24 2822
本文是一篇手记，记录了我对负二项分布及其应用的理解。目录如下：1.理解"负"的含义知乎回答有提到：https://www.zhihu.com/question/24253978?sort=created文档中有提到：http://www.johndcook.com/negativebinomial.pdf负二项级数：https://brilliant.org/wiki/negative-binom...
基于Python回归模型的多重共线性分析
DL11007的博客
02-24 2918
基于Python回归模型的多重共线性分析
回归分析——简单线性回归实例讲解（SPSS）
热门推荐
qq_39783601的博客
09-06 7万+
什么是回归分析？ 回归分析是研究自变量与因变量之间数量变化关系的一种分析方法，它主要是通过因变量Y与影响它的自变量Xi（i1,2,3…）之间的回归模型，衡量自变量Xi对因变量Y的影响能力的，进而可以用来预测因变量Y的发展趋势。 ????????????相关分析与回归分析的联系与区别： 联系： 两者均为研究与测度两个或两个以上变量之间关系的方法。 在实际工作中，一般先进行相关分析，计算相关函数，然后建立回归模型，最后用回归模型进行推算或预测。 区别： ✨相关分析研究的都是随机变量，并且不分自变量与因变量；
线性回归分析：从简单问题到多元线性回归
线性回归分析 线性回归分析是机器学习中的一种重要方法，用于处理回归问题。线性回归模型可以预测连续值，广泛应用于各个领域。下面我们将详细介绍线性回归分析的相关知识点： 一、什么是线性回归？ 线性回归是一...
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

积跬步，慕至千里

你的鼓励将是我创作的最大动力！

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值