之前上学时计量经济学的模型实现总是用Eviews等软件实现。但是对于点击鼠标得到结果的方式,总是让自己感觉没有参与模型建立的过程。所以准备利用python写代码从底层,进行编写进行计量经济分析,让自己更了解预算过程。暂时准备写以下几篇,后面再慢慢补充;

  1. 多元线性回归和显著性检验(参数估计、T检验、F检验、拟合优度)
  2. 多重共线性(导致结果、检验——方差膨胀因子、补救措施——岭回归)
  3. 异方差(导致结果、检验——White、补救措施——广义线性回归)
  4. 自相关(导致结果、检验——D-W、补救措施——广义线性回归)

写的比较仓促,代码中如有错误欢迎指正!

一、计量经济分析

个人理解计量经济学可以认为是统计学的一个应用或分支。我们知道在统计层面上,因素的关系大概可以分两种:因果关系相关关系

现在十分流行的机器学习、深度学习中关注的事物间的相关关系(比如身高和体重),事物间不一定需要是一个导致另一个关系,只要需要再生产环境相关关系是有效的就可以拿来用。

而计量经济分析则主要侧重的是因果关系,因为我们的主要目的是分析然后进行预测。既然是分析,我们就需要知道因素之间的具体关系,一个因素的变化是怎样导致另一个因素的变化。这也是计量与机器学习的最明显差异。

既然计量经济分析侧重的是因果关系,因此我们对模型的检验就会更加深入。具体的,比如在机器学习中,我们得到一个模型,模型再训练集、测试集、时间外样本集等AUC等指标均表现良好,并且线上可以稳定运行。我们就可以拿来使用。而在计量经济学中,我们需要对模型进行多种检验,确保模型中的参数意义明确和有效。

计量经济分析的方法步骤可以总结如下:1.模型设定(设定一个合理模型)——>2.参数估计(最小二乘等)——>3.模型检验(参数显著性、是否有多重共线、异方差,自相关和参数符号是否符合经济意义等)——>4.模型应用(利用模型进行因素间内在结构分析或预测)

本系列只对原理和形式进行简单介绍,然后利用python进行实现。如果想要入门,可以看看国内本科的一些计量经济学教材。想要深入一些,可以看看伍德里奇的《计量经济学导论》。若想要深究,可以拜读格林的《计量经济分析》。

二、多元线性回归

2.1.多元线性回归形式
一个因素的变化通常是由多个因素共同决定的,比如经典的生产函数认为,产出是由资本和劳动力决定的。如果我们研究多个因素对另一个变量的作用,我们设置如下多元线性回归模型:

3 python 计量经济学 python做计量经济学_数学建模
其矩阵形式为:
3 python 计量经济学 python做计量经济学_数学建模_02

3 python 计量经济学 python做计量经济学_python_03被称为应变量或被解释变量,3 python 计量经济学 python做计量经济学_多元线性回归_04是解释变量,3 python 计量经济学 python做计量经济学_3 python 计量经济学_05扰动项,3 python 计量经济学 python做计量经济学_数学建模_063 python 计量经济学 python做计量经济学_多元线性回归_04前的参数,表示3 python 计量经济学 python做计量经济学_多元线性回归_04一单位变动引起的3 python 计量经济学 python做计量经济学_3 python 计量经济学_09变动程度,如果3 python 计量经济学 python做计量经济学_多元线性回归_043 python 计量经济学 python做计量经济学_3 python 计量经济学_09都是对数形式:3 python 计量经济学 python做计量经济学_多元线性回归_123 python 计量经济学 python做计量经济学_3 python 计量经济学_13,则3 python 计量经济学 python做计量经济学_数学建模_06表示的经济学中弹性。

我们估计模型使用的最小二乘法,但是若保证估计是有效的,模型必须满足线性回归的六大假定(不细说了,每违反一种假定就会犯响应错误,就需要相应补救):

1.模型符合线性模式
2.X满秩(无多重共线)
3.零均值价值:3 python 计量经济学 python做计量经济学_统计学_15(自变量外生)
4.同方差:3 python 计量经济学 python做计量经济学_统计学_16
5.无自相关:3 python 计量经济学 python做计量经济学_统计学_17
6.球形扰动:ε_i是正态分布

2.2.多元线性回归参数估计

参数估计使用最小二乘法,得到参数估计为:
3 python 计量经济学 python做计量经济学_统计学_18
3 python 计量经济学 python做计量经济学_统计学_19 (后面t检验会用到,3 python 计量经济学 python做计量经济学_3 python 计量经济学_20可用残差3 python 计量经济学 python做计量经济学_3 python 计量经济学_21估计得到)

2.3.多元线性回归检验

对于多元线性回归的检验,和统计的假设检验一样,就是构造符合几大分布的统计量(正态、T分布、F分布、卡方分布)的统计量,然后进行假设检验,其中零假设为参数等于0(或方程不显著),我们的目标是拒绝零假设,这样参数才是有意义,才可以进行下一步的分析

2.3.1.单个参数检验:
构造统计量:T=3 python 计量经济学 python做计量经济学_统计学_22
T符合3 python 计量经济学 python做计量经济学_多元线性回归_233 python 计量经济学 python做计量经济学_数学建模_24为样本量,3 python 计量经济学 python做计量经济学_数学建模_25为参数个数,查表可知参数是否显著等于0

2.3.2.整体检验:

1.拟合优度3 python 计量经济学 python做计量经济学_多元线性回归_26
其中3 python 计量经济学 python做计量经济学_3 python 计量经济学_27是回归平方和,3 python 计量经济学 python做计量经济学_统计学_28是总平方和,3 python 计量经济学 python做计量经济学_python_29,即总平方和=回归平方和+误差平方和。3 python 计量经济学 python做计量经济学_多元线性回归_30是一个0~1之间的数字,3 python 计量经济学 python做计量经济学_多元线性回归_30越大,代表方程拟合越好(实际上,随着变量增多,3 python 计量经济学 python做计量经济学_多元线性回归_30一定会增大,仅此还有加入了变量个数惩罚的调整拟合优度)。

2.F检验:3 python 计量经济学 python做计量经济学_python_33
因为3 python 计量经济学 python做计量经济学_多元线性回归_30并没有自由度的概念,因此我们无法查表得知3 python 计量经济学 python做计量经济学_多元线性回归_30多大算大,因此我们构造3 python 计量经济学 python做计量经济学_3 python 计量经济学_36统计量,3 python 计量经济学 python做计量经济学_3 python 计量经济学_36统计量服从3 python 计量经济学 python做计量经济学_3 python 计量经济学_38,于是我们可以查表得到拒绝域,从而判断方程整体是否显著。
多说一句,F分布是两个卡方分布除以各自自由度的比构成的。F检验的思想非常常用,比如方差分析等都是利用F检验的思想,将方差分为不同的类别,查看不同类别的方差(例如组内,组间)是否存在不同。

三、多元线性回归的python实现

我们使用python的sklearn中的房价数据集进行多元线性回归。

from sklearn.datasets import fetch_california_housing as fch

data=fch() #导入数据
house_data=pd.DataFrame(data.data) #将自变量转换成dataframe格式,便于查看
house_data.columns=data.feature_names  #命名自变量
house_data['value'] = data.target
house_data = house_data[0:100]
house_data = house_data[['value', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population']]
print(house_data.head())

我们想用,房屋年龄,房间数量,洗漱间数量和人口来预测一个房子的房价。

3 python 计量经济学 python做计量经济学_统计学_39


接下来我们利用python,进行多元线性回归的参数估计和检验进行底层编写,我们将多元线性回归整体封装为一个类,将多元线性回归模型设计成一个对象:

# -*- coding: utf-8 -*-
"""
Created on Fri Apr  3 15:40:38 2020

@author: nbszg
"""

import numpy as np
import pandas as pd
import scipy.stats as st

class mul_linear_model(object):
    #intercept定义是否有截距项,在计量经济分析中,基本都要有截距项,否则计算出的拟合优度是没有意义的
    def __init__(self, y, X, intercept=True):
        data_x = X.copy()
        data_y = y.copy()
        self.data_x = data_x
        self.y = np.array(data_y)
        self.intercept = intercept
        
        if intercept == False:
            self.columns = data_x.columns
            self.X = np.mat(data_x)
        else:
            #插入截距项
            data_x.insert(0,'intercept',np.ones(len(data_x)))   
            #每个变量名称
            self.columns = data_x.columns
            #转换为矩阵,方便后面运算
            self.X = np.mat(data_x)
        
        #X的转置,方便后面运算
        self.XT=np.mat(self.X).T 
        self.N = self.X.shape[0]   #样本量
        self.K = self.X.shape[1]  #变量个数
        
    
    #拟合模型
    def fit(self, output=True):
        #使用最小二乘法(X'X)的逆成X'y
        self.b = np.array(np.dot(np.dot(np.linalg.inv(np.dot(self.XT,self.X)),self.XT),np.mat(self.y).T))
        if output:
            for i in range(self.K):
                print("variable {0}'s cofe estimate is : {1}".format(self.columns[i], self.b[i][0]))
        
        self.S = np.array(np.linalg.inv(np.dot(self.XT,self.X)))
        return self.b
    
    #预测函数
    def predict(self, X_p):
        #直接使用条件期望作y=Xb为预测值
        y_hat =np.array(np.dot(X_p,self.b))
        return y_hat
    
    #参数T检验
    def T_test(self, alpha):
        #计算预测值
        y_hat = np.array(np.dot(self.X,self.b).T)[0]
        #计算e'e,即残差平方和
        error_square = np.sum(np.square(y_hat-self.y))
        #计算扰动项的方差估计,即残差平方和
        sigma_hat = error_square/(self.N-self.K)
        #计算(X'X)的逆,用于估计参数b的方差
        #S = self.S
        for i in range(self.K):
            sk = self.S[i][i]
            #b的方差估计为sigma_hat * sk,构造t统计量,服从自由度为n-k的T分布
            Tk = self.b[i][0]/np.sqrt(sigma_hat * sk)
            #查表进行检验
            if abs(Tk) > abs(st.t.ppf(alpha/2, df = self.N-self.K)):
                print("T of variable {0} is: {1}, refuse H0".format(self.columns[i], Tk))
            else:
                print("T of variable {0} is: {1}, can't refuse H0".format(self.columns[i], Tk))
    #拟合优度计算
    def R_square(self):
        y_ba = np.mean(self.y)
        y_hat = np.array(np.dot(self.X,self.b).T)[0]
        #计算总平方和
        self.SST = np.sum(np.square(self.y-y_ba))
        #计算误差平方和
        self.SSE = np.sum(np.square(self.y-y_hat))
        #计算离差平方和
        self.SSR = self.SST - self.SSE
        #利用拟合优度公式计算拟合优度
        self.R_2 = self.SSR/self.SST
        
        return self.R_2
    
    #整体F检验
    def F_test(self, alpha):
        #先求出拟合优度
        self.R_square()
        #利用公式求F统计量,附送自由优度为(K-1, N-K)的F分布
        F = (self.R_2/(self.K-1))/((1-self.R_2)/(self.N-self.K))
        #查表进行检验
        if abs(F) > abs(st.f.ppf(1-alpha, dfn = self.K-1, dfd = self.N-self.K)):
            print("F of function is: {0}, refuse H0".format(F))
        else:
            print("F of function is: {0}, can't refuse H0".format(F))

编写好了多元线性回归,我们使用刚才的房屋数据进行试验:

linear_model = mul_linear_model(house_data['value'],house_data[['HouseAge', 'AveRooms', 'AveBedrms', 'Population']])

linear_model.fit()
linear_model.T_test(0.05)
print(linear_model.R_square())
linear_model.F_test(0.05)

得到的结果如下:

3 python 计量经济学 python做计量经济学_多元线性回归_40


想知道编写的对不对,我们直接通statsmodels.formula.apiols模块进行验证,两行代码就可以解决刚才编写的所有问题:

import statsmodels.api as sm #最小二乘
from statsmodels.formula.api import ols #加载ols模型

lm=ols('value~ HouseAge + AveRooms + AveBedrms + Population',data=house_data).fit()
print(lm.summary())

3 python 计量经济学 python做计量经济学_python_41


检验结果如上图,其中红框的部分是我们的估计参数和检验结果,可以看到使用statsmodels包与我们从底层编写的计算结果是一致的。证明了我们代码的正确性。