python里很多模块都有OLS的实现,之前总结过一次,详见《从零开始学量化(五):用Python做回归》。今天这个是自己用numpy实现OLS,WLS的一些内容。
自己写的好处是比内置的要快一些,倒也不是说内置的代码写的不好,而是内置的函数一般比较完善,会算出来很多东西,但有的时候就只需要个回归系数,需要个预测值,结果内置函数给你整一堆t值,r方,mse,summary拖慢你的速度,也没啥意义。
01
OLS的beta
先从OLS开始,OLS比较常用的是statsmodels里的OLS模块,除此之外,numpy里本身也有一个函数能实现。这里从定义出发直接算一个,另外做一个简单测试对比numpy和statsmodels里的速度差异。
OLS的beta定义:
公式推导就省略了,随便找概率书都有,直接代码。
def getOLSBeta(x,y):
x = np.hstack([x.reshape(-1,1),np.ones([x.shape[0],1])])
beta = np.linalg.inv(x.T.dot(x)).dot(x.T).dot(y)
return beta
接下来分别用这个函数和statsmodel的OLS模块算beta,做一个简单测试。
首先随机生成两个(10000,1)的向量作为x,y:
np.random.seed(1)
x = np.random.uniform(0,1, (10000,))
y = np.random.uniform(0,1, (10000,))
然后分别用这两个函数算beta,重复100次,计算总的时间,结果对比如下
OLS需要86735microseconds,numpy函数的结果:
自定义的函数需要14958microseconds,自定义的要快6-7倍。
02
WLS的beta
同样的道理,定义WLS的beta函数,这个就不做测试了,不用想都知道肯定是比statsmodel里的WLS更快一些。WLS的beta表达式:
python代码如下
def getWLSBeta(x,y,w):
x = np.hstack([x.reshape(-1,1),np.ones([x.shape[0],1])])
x1 = np.sqrt(w.reshape(-1,1)) * x
y1 = np.sqrt(w) * y
beta = np.linalg.inv(x1.T.dot(x1)).dot(x1.T).dot(y1)
return beta
03
OLS的预测值
OLS的预测值,有两种,一般大家只看点预测,也就是拟合出来的值,这个很简单,不管是新来的点还是回归数据里的点,直接和beta乘一下就出来了,没啥难度。
还有一个是区间预测,给定一个置信度alpha,预测值是会有一个范围的,点预测实际上是预测区间的中点,预测区间有多宽,取决于方程本身的拟合程度,拟合的越好,区间越窄,说明预测的结果准确性越高。所以这里用numpy实现OLS的区间预测,推导也不写了,直接写公式
numpy计算的函数如下
def getOLSfitted(x,y,x_new = None,alpha = 0.05):
x = np.hstack([x.reshape(-1,1),np.ones([x.shape[0],1])])
y_pred_se = np.linalg.inv(x.T.dot(x))
beta = y_pred_se.dot(x.T).dot(y)
y_pred = x.dot(beta)
sigma2_est = np.sum((y - y_pred)**2) / (x.shape[0] - x.shape[1])
if x_new == None:
x_new = x
y_pred_new = y_pred
else:
x_new = np.hstack([x_new.reshape(-1,1),np.ones([x_new.shape[0],1])])
y_pred_new = x_new.dot(beta)
y_pred_se = sigma2_est + (x_new * np.dot(y_pred_se * sigma2_est,x_new.T).T).sum(1)
y_pred_se = np.sqrt(y_pred_se)
y_pred_lower = y_pred_new - stats.t.ppf(q = 1- alpha/2 , df = x.shape[0] - x.shape[1])*y_pred_se
y_pred_upper = y_pred_new + stats.t.ppf(q = 1- alpha/2 , df = x.shape[0] - x.shape[1])*y_pred_se
return y_pred_se,y_pred,y_pred_lower,y_pred_upper
解释一下,函数输入x,y为拟合方程的,x_new为用来预测的点,如果不输入的话返回的是x的拟合值,alpha是置信度。返回值有四个,se,y的点预测,y的区间预测下限,y的区间预测上限。写的比较粗糙,自己用的话可以再根据需求改一改。
另外statsmodel里也可以直接求OLS的预测区间,需要用到wls_prediction_std函数,所以还是之前的那个例子,做一个测试。
wls_prediction_std的结果
需要130651microseconds。
自定义函数的结果:
需要81749microseconds,要快一些,但还有提升空间。
03
WLS的预测值
WLS同样的道理,可以直接用wls_prediction_std计算,另外也自己定义了getWLSfitted函数,对比结果如下,还是之前的例子,再随机生成一个权重w。
np.random.seed(1)
x = np.random.uniform(0,1, (10000,))
y = np.random.uniform(0,1, (10000,))
w = np.random.uniform(0,1, (10000,))
wls_prediction_std函数的结果
需要166555microseconds。
自定义函数的结果
需要98703microseconds,所以还是自定义的要快一些,但还有提升空间。
本文所有代码包括getWLSfitted的函数定义请在后台回复“WLS“获取。
当然除了以上这些, 还可以自己定义算OLS 的t值,OLS的控制等等,这里就略过啦。