The following code will attempt to replicate the results of the
numpy.linalg.lstsq() function in Numpy. For this exercise, we will be using a cross sectional data set provided by me in .csv format called “cdd.ny.csv”, that has monthly cooling degree data for New York state. The data is available here (File –> Download).
The OLS regression equation:
where a white noise error term. For this example the population-weighted Cooling Degree Days (CDD) (CDD.pop.weighted), and CDD measured at La Guardia airport (CDD.LGA). Note: this is a meaningless regression used solely for illustrative purposes.
Recall that the following matrix equation is used to calculate the vector of estimated coefficients of an OLS regression:
where the matrix of regressor data (the first column is all 1′s for the intercept), and the vector of the dependent variable data.
Matrix operators in Numpy
matrix()coerces an object into the matrix class.
.Ttransposes a matrix.
dot(X,Y)is the operator for matrix multiplication (when matrices are 2-dimensional; see here).
.Itakes the inverse of a matrix. Note: the matrix must be invertible.
Back to OLS
The following code calculates the 2 x 1 matrix of coefficients, :
## load in required modules: import numpy as np import csv ## read data into a Numpy array df1 = csv.reader(open('/your/file/path/cdd.ny.csv', 'rb'),delimiter=',') b1 = np.array(list(df1))[1:,3:5].astype('float') nrow = b1.shape intercept = np.ones( (nrow,1) ) b2 = b1[:,0].reshape(-1, 1) X = np.concatenate((intercept, b2), axis=1) Y = b1[:,1].T ## X and Y arrays must have the same number of columns for the matrix multiplication to work: print(X.shape) print(Y.shape) ## Use the equation above (X'X)^(-1)X'Y to calculate OLS coefficient estimates: bh = np.dot(np.linalg.inv(np.dot(X.T,X)),np.dot(X.T,Y)) print bh ## check your work with Numpy's built in OLS function: z,resid,rank,sigma = np.linalg.lstsq(X,Y) print(z)
Calculating Standard Errors
To calculate the standard errors, you must first calculate the variance-covariance (VCV) matrix, as follows:
The VCV matrix will be a square k x k matrix. Standard errors for the estimated coefficients are found by taking the square root of the diagonal elements of the VCV matrix.
## Calculate vector of residuals res = as.matrix(women$weight-bh-bh*women$height) res = Y-(bh+X[:,1]*bh) ## Define n and k parameters n = nrow k = X.shape ## Calculate Variance-Covariance Matrix VCV = np.true_divide(1,n-k)*np.dot(np.dot(res.T,res),np.linalg.inv(np.dot(X.T,X))) ## Standard errors of the estimated coefficients stderr = np.sqrt(np.diagonal(VCV))
Now you can check the results above using the
lm() function in R:
df1 = read.csv('/your/file/path/cdd.ny.csv',header=T) coef(lm(CDD.pop.weighted ~ CDD.LGA,data=df1)) ## (Intercept) CDD.LGA ## -7.6210191 0.5937734 summary(lm(CDD.pop.weighted ~ CDD.LGA,data=df1))