<jupyter_start><jupyter_text>Creating an OLS function in Python Linear Algebra in _numpy_We reviewed chapter 7 to understand how to generate an inverse matrix.To transform matrix $X$ into $X^{-1}$. To do so, we must transform $[XI]$ into $[IX^{-1}]$. This is accomplished in chapter 7 of the text book. We are going to check the result using _numpy_.<jupyter_code>#invertMatrix.py
import numpy as np
# create array to be transformed into vectors
x1 = np.array([1,2,1])
x2 = np.array([4,1,5])
x3 = np.array([6,8,6])
print("Array 1:", x1, sep="\n")
print("Array 2:", x2, sep="\n")
print("Array 3:", x3, sep="\n")
x1 = np.matrix(x1)
x2 = np.matrix(x2)
x3 = np.matrix(x3)
print("Vector 1", x1, sep = "\n")
print("Vector 2", x2, sep = "\n")
print("Vector 3", x3, sep = "\n")<jupyter_output>Vector 1
[[1 2 1]]
Vector 2
[[4 1 5]]
Vector 3
[[6 8 6]]<jupyter_text>Join x1, x2, and x3 using _np.concatenate()_. We define _axis = 0_ so that the rows are stacked vertically, not  horizontally.<jupyter_code>X = np.concatenate((x1, x2, x3), axis = 0)
print("X:", X, sep = "\n")<jupyter_output>X:
[[1 2 1]
 [4 1 5]
 [6 8 6]]<jupyter_text>Now that we have defined our $X$ matrix, we can generate $X^{-1}$ by using _X.getI()_<jupyter_code>X_inverse = X.getI()
print("X Inverse:", X_inverse, sep = "\n")<jupyter_output>X Inverse:
[[-8.5000000e+00 -1.0000000e+00  2.2500000e+00]
 [ 1.5000000e+00 -7.6861594e-17 -2.5000000e-01]
 [ 6.5000000e+00  1.0000000e+00 -1.7500000e+00]]<jupyter_text>Next, use _np.round()_ to clearly present variables that include digits to the nth place (17).<jupyter_code>X_inverse = np.round(X.getI(), 2)
print("X Inverse:", X_inverse, sep = "\n")<jupyter_output>X Inverse:
[[-8.5  -1.    2.25]
 [ 1.5  -0.   -0.25]
 [ 6.5   1.   -1.75]]<jupyter_text>Finally, we should also see how to generate the transpose of a matrix. The transpose switches the index values $i$ and $j$. In otherwords, switch the row indexes with the column indexes.<jupyter_code>print("X:", X, sep = "\n")
X_transpose = X.getT()
print("X Transpose:", X_transpose, sep = "\n")<jupyter_output>X:
[[1 2 1]
 [4 1 5]
 [6 8 6]]
X Transpose:
[[1 4 6]
 [2 1 8]
 [1 5 6]]<jupyter_text>Diversion on Classes in PythonCreate a sample class that contains math functions.<jupyter_code>class Arithmetic:
    def __init__(self, val):
        self.name = "arithmetic"
        self.val = val
    def add(self, *args):
        total = 0
        try:
            for arg in args:
                total += arg
                
            return total
        except:
            ("pass list or tuple of ints or floats to add()")

    def multiply(self, *args):
        product = 1
        try:
            for arg in args: 
                product *= arg
            return product
        except:
            print("pass only int or float to multiply")
        
arithmetic = Arithmetic(3)
print(arithmetic)
print(arithmetic.add(1,2,3,4))
print(arithmetic.multiply(1,2,3,4))
print(arithmetic.name)
print(arithmetic.val)
print(arithmetic.__dict__)<jupyter_output><__main__.Arithmetic object at 0x00000204321CA400>
10
24
arithmetic
3
{'name': 'arithmetic', 'val': 3}<jupyter_text>Building an OLS Regression Function$\beta = (X'X)^{-1}X'Y$An OLS regression posits a relationship between a dependent variable, $Y$, and independent variables - each $x_i$ - contained in $X$.<jupyter_code>#econFreedomRegression.py
import pandas as pd

data = pd.read_csv("cleanedEconFreedomData.csv", 
    index_col = "Country Name")
data
from regression import Regression
reg = Regression()
print(reg)
reg.stats.mean([1,2,3,4,5,6])<jupyter_output>You created an instance of stats()
<regression.Regression object at 0x0000020431D48DF0><jupyter_text>Our _regression.py_ file so far:<jupyter_code>class Regression:
    def __init__(self):
        self.stats = stats()
    
    # if constant is True, add a column of ones to 
    # estimate a constant
    def OLS(self, reg_name, data, y_name, beta_names,
            min_val = 0, max_val = None,
            constant = True):
    
        self.min_val = min_val
        if max_val != None:
            self.max_val = max_val
        else:
            self.max_val = len(data)
        self.reg_name = reg_name
        self.y_name = y_name
        
        self.beta_names = copy.copy(beta_names)
        self.data = data.copy()
        if constant:
            self.add_constant()

                    
    def add_constant(self):
        self.data["Constant"] = 1
        self.beta_names.append("Constant")<jupyter_output><empty_output><jupyter_text>Next, let's call the regression and view the data that now has a column of ones.<jupyter_code>#econFreedomRegression.py
# . . .
# dependent variable -- determined by x_vars
y_var = ["GDP per Capita (PPP)"]
# independent variables
x_vars = ["Trade Freedom", "Property Rights", 
          "Inflation (%)",
          "Public Debt (% of GDP)"]

reg.OLS(reg_name = "GDP per Capita", data = data, 
        y_name = y_var, beta_names = x_vars)
reg.data<jupyter_output><empty_output><jupyter_text>Next we will estimate the influence of each x variable (exogenous) on the y variable (endogenous):    $\beta = (X'X)^{-1}X'Y$$Y =$ Observations for Endogenous Variable$X =$ Observations for Exogenous Variables$X' =$ $X$-transpose$(X'X)^{-1} =$ Inverse of $X'X$The _regression.py_ file should include the following:<jupyter_code>class Regression:
    def __init__(self):
        self.stats = stats()
    
    # if constant is True, add a column of ones to 
    # estimate a constant
    def OLS(self, reg_name, data, y_name, beta_names,
            min_val = 0, max_val = None,
            constant = True):
    
        self.min_val = min_val
        if max_val != None:
            self.max_val = max_val
        else:
            self.max_val = len(data)
        self.reg_name = reg_name
        self.y_name = y_name
        
        self.beta_names = copy.copy(beta_names)
        self.data = data.copy()
        if constant:
            self.add_constant()
        self.build_matrices()
        self.estimate_betas_and_yhat()

    def add_constant(self):
        self.data["Constant"] = 1
        self.beta_names.append("Constant")
        
    def build_matrices(self):
        # Transform dataframes to matrices
        self.y = np.matrix(self.data[self.y_name][self.min_val:self.max_val])
        # create a k X n nested lest containing vectors from each exog var
        self.X = np.matrix(self.data[self.beta_names])
        self.X_transpose = np.matrix(self.X).getT()
        # (X'X)**-1
        X_transp_X = np.matmul(self.X_transpose, self.X)
        self.X_transp_X_inverse = X_transp_X.getI()
        # X'y
        self.X_transp_y = np.matmul(self.X_transpose, self.y)

    def estimate_betas_and_yhat(self):
        # betas = (X'X)**-1 * X'y
        self.betas = np.matmul(self.X_transp_X_inverse, self.X_transp_y)
        # y-hat = X * betas
        self.y_hat = np.matmul(self.X, self.betas)
        self.data[self.y_name[0] + " estimator"] =\
            [i.item(0) for i in self.y_hat]
        # create a table for the estimates
        self.estimates = pd.DataFrame(self.betas, index = self.beta_names,
                                      columns = ["Coefficient"])
        # identify y vairiable in index
        self.estimates.index.name = "y = " + self.y_name[0]<jupyter_output><empty_output><jupyter_text>Save the file and rerun your code so that you can call these statistics.<jupyter_code>reg.betas<jupyter_output><empty_output><jupyter_text>Next estimate the beta coefficients that indicate the influence of the x (exogenous) variables on the y (endogenous) variable.$y ̂ = X\beta$<jupyter_code>reg.estimates<jupyter_output><empty_output><jupyter_text>Next, we will use the sum squared statistics (ssr, sse, and sst) to calculate different variables that reflect the goodness of fit of the regression and its parts.We define these values as follows:$SSR = \sum_{i=0}^{n} (y ̂ _{i} - y ̅ )^2$$SSE = \sum_{i=0}^{n} (y_{i} - y ̂ _{i})^2$$SST = \sum_{i=0}^{n} (y_{i} - y ̅ _{i})^2$It happens that the sum of the squared distances between the estimated values and mean of observed values and the squared distances between the observed and estimated values add up to the sum of the squared distances between the observed values and the mean of observed values. We indicate this as:$SST = SSR + SSE$The _regression.py_ should look as follows:<jupyter_code>class Regression:
    def __init__(self):
        self.stats = stats()
    

    # if constant is True, add a column of ones to 
    # estimate a constant
    def OLS(self, reg_name, data, y_name, beta_names,
            min_val = 0, max_val = None,
            constant = True):
    
        self.min_val = min_val
        if max_val != None:
            self.max_val = max_val
        else:
            self.max_val = len(data)
        self.reg_name = reg_name
        self.y_name = y_name
        
        self.beta_names = copy.copy(beta_names)
        self.data = data.copy()
        if constant:
            self.add_constant()
        self.build_matrices()
        self.estimate_betas_and_yhat()
        self.calculate_regression_stats()
        
    def add_constant(self):
        self.data["Constant"] = 1
        self.beta_names.append("Constant")
        
    def calculate_regression_stats(self):
        self.sum_square_stats()
        
    def build_matrices(self):
        # Transform dataframes to matrices
        self.y = np.matrix(self.data[self.y_name][self.min_val:self.max_val])
        # create a k X n nested lest containing vectors from each exog var
        self.X = np.matrix(self.data[self.beta_names])
        self.X_transpose = np.matrix(self.X).getT()
        # (X'X)**-1
        X_transp_X = np.matmul(self.X_transpose, self.X)
        self.X_transp_X_inverse = X_transp_X.getI()
        # X'y
        self.X_transp_y = np.matmul(self.X_transpose, self.y)

    def estimate_betas_and_yhat(self):
        # betas = (X'X)**-1 * X'y
        self.betas = np.matmul(self.X_transp_X_inverse, self.X_transp_y)
        # y-hat = X * betas
        self.y_hat = np.matmul(self.X, self.betas)
        self.data[self.y_name[0] + " estimator"] =\
            [i.item(0) for i in self.y_hat]
        # create a table for the estimates
        self.estimates = pd.DataFrame(self.betas, index = self.beta_names,
                                      columns = ["Coefficient"])
        # identify y vairiable in index
        self.estimates.index.name = "y = " + self.y_name[0]

    def sum_square_stats(self):
        ssr_list = []
        sse_list = []
        sst_list = []
        mean_y = self.stats.mean(self.y).item(0)
        for i in range(len(self.y)):
            # ssr is sum of squared distances between the estimates
            # and the avergage of y values
            y_hat_i = self.y_hat[i]
            y_i = self.y[i]
            r = y_hat_i - mean_y
            e = y_i - y_hat_i
            t = y_i - mean_y
            ssr_list.append((r) ** 2)
            sse_list.append((e) ** 2)
            sst_list.append((t) ** 2)
        
        # since the ssr, sse, and sst use values from 
        # matrices, select the value within the resultant
        # matrix using matrix.item(0)
        self.ssr = self.stats.total(ssr_list).item(0)
        self.sse = self.stats.total(sse_list).item(0)
        self.sst = self.stats.total(sst_list).item(0)<jupyter_output><empty_output><jupyter_text>Save the file and rerun your code so that you can call these statistics.<jupyter_code>reg.ssr, reg.sse, reg.sst<jupyter_output><empty_output><jupyter_text>We estimate the variance of the error term of the estimator for the dependent variable. $\sigma^2 = \frac{SSE}{n-k}$$n = $number of observations$k = $number of independent variablesAn increase in the number of exogenous variables tends ot increase the fit of a model. By dividing the $SSE$ by degrees of freedom, $n-k$ , improvements in fit that result from increases in the number of variables are offset in part by a reduction in degrees of freedom.We can calculate the covariance matrix:$\sigma^2 (X'X)^{-1}$The diagonals of the covariance matrix represent the standard errors squared for each variable. The standard what we will use to calculate t-statistics and p-values. The t-statisitic of a coefficient is found by comparing the size of the estimated coefficient to its standard error:$tstat_{\beta_i}=\frac{\beta_i}{SE_{\beta_i}}$The larger the coefficient compared to the error, the more reliable is the statistic, as implied by a large t-stat. We draw the p-value associated with a particular t-stat from a table in light of the degrees of freedom associated with the regression. The p-values provides a rating of the estimate in light of the t-stat together with the number of degrees of freedom.<jupyter_code>class Regression:
    def __init__(self):
        self.stats = stats()
    

    # if constant is True, add a column of ones to 
    # estimate a constant
    def OLS(self, reg_name, data, y_name, beta_names,
            min_val = 0, max_val = None,
            constant = True):
    
        self.min_val = min_val
        if max_val != None:
            self.max_val = max_val
        else:
            self.max_val = len(data)
        self.reg_name = reg_name
        self.y_name = y_name
        
        self.beta_names = copy.copy(beta_names)
        self.data = data.copy()
        if constant:
            self.add_constant()
        self.build_matrices()
        self.estimate_betas_and_yhat()
        self.calculate_regression_stats()
        
    def add_constant(self):
        self.data["Constant"] = 1
        self.beta_names.append("Constant")
        
    def calculate_regression_stats(self):
        self.sum_square_stats()
        self.calculate_degrees_of_freedom()
        self.calculate_estimator_variance()
        self.calculate_covariance_matrix()
        self.calculate_t_p_error_stats()
        
    def build_matrices(self):
        # Transform dataframes to matrices
        self.y = np.matrix(self.data[self.y_name][self.min_val:self.max_val])
        # create a k X n nested lest containing vectors from each exog var
        self.X = np.matrix(self.data[self.beta_names])
        self.X_transpose = np.matrix(self.X).getT()
        # (X'X)**-1
        X_transp_X = np.matmul(self.X_transpose, self.X)
        self.X_transp_X_inverse = X_transp_X.getI()
        # X'y
        self.X_transp_y = np.matmul(self.X_transpose, self.y)

    def estimate_betas_and_yhat(self):
        # betas = (X'X)**-1 * X'y
        self.betas = np.matmul(self.X_transp_X_inverse, self.X_transp_y)
        # y-hat = X * betas
        self.y_hat = np.matmul(self.X, self.betas)
        self.data[self.y_name[0] + " estimator"] =\
            [i.item(0) for i in self.y_hat]
        self.data["residual"] = self.data[self.y_name[0]].sub(self.data[self.y_name[0] + " estimator"])
        # create a table for the estimates
        self.estimates = pd.DataFrame(self.betas, index = self.beta_names,
                                      columns = ["Coefficient"])
        # identify y vairiable in index
        self.estimates.index.name = "y = " + self.y_name[0]

    def sum_square_stats(self):
        ssr_list = []
        sse_list = []
        sst_list = []
        mean_y = self.stats.mean(self.y).item(0)
        for i in range(len(self.y)):
            # ssr is sum of squared distances between the estimates
            # and the avergage of y values
            y_hat_i = self.y_hat[i]
            y_i = self.y[i]
            r = y_hat_i - mean_y
            e = y_i - y_hat_i
            t = y_i - mean_y
            ssr_list.append((r) ** 2)
            sse_list.append((e) ** 2)
            sst_list.append((t) ** 2)
        
        # since the ssr, sse, and sst use values from 
        # matrices, select the value within the resultant
        # matrix using matrix.item(0)
        self.ssr = self.stats.total(ssr_list).item(0)
        self.sse = self.stats.total(sse_list).item(0)
        self.sst = self.stats.total(sst_list).item(0)
        
    def calculate_degrees_of_freedom(self):
        # Degrees of freedom compares the number of observations to the number  
        # of exogenous variables used to form the prediction          
        self.lost_degrees_of_freedom = len(self.estimates)
        self.num_obs = self.max_val + 1 - self.min_val
        self.degrees_of_freedom = self.num_obs - self.lost_degrees_of_freedom

    def calculate_estimator_variance(self):
        # estimator variance is the sse normalized by the degrees of freedom  
        # thus, estimator variance increases as the number of exogenous  
        # variables used in estimation increases(i.e., as degrees of freedom   
        # fall) 
        self.estimator_variance = self.sse / self.degrees_of_freedom
    
    def calculate_covariance_matrix(self):
        # Covariance matrix will be used to estimate standard errors for  
        # each coefficient.  
        # estimator variance * (X'X)**-1  
        self.cov_matrix = float(self.estimator_variance) * self.X_transp_X_inverse
        self.cov_matrix = pd.DataFrame(self.cov_matrix,
                                  columns = self.beta_names,
                                  index = self.beta_names)
        
    def calculate_t_p_error_stats(self):
        results = self.estimates
        stat_sig_names = ["SE", "t-stat", "p-value"]
        # create space in data frame for SE, t, and p
        for stat_name in stat_sig_names:
            results[stat_name] = np.nan
        # generate statistic for each variable
        for var in self.beta_names:
            # SE ** 2 of coefficient is found in the diagonal of the cov_matrix
            results.loc[var]["SE"] = self.cov_matrix[var][var] ** (1/2)     
            # t-stat = Coef / SE
            results.loc[var]["t-stat"] = \
                results["Coefficient"][var] / results["SE"][var]
            # p-values is estimated using a table that transforms t-stat in   
            # light of degrees of freedom  
            # 2 is for 2 tail...
            # 5 is to round to 5 decimal places
            results.loc[var]["p-value"] = np.round(
                t.sf(np.abs(results.loc[var]["t-stat"]), 
                     self.degrees_of_freedom +1) * 2, 5)
        ratings = [.05, .01, .001]
        significance = ["" for name in self.beta_names]
        for i in range(len(self.beta_names)):
            var_name = self.beta_names[i]
            for rating in ratings:
                if results.loc[var_name]["p-value"] < rating:
                    significance[i] = significance[i] + "*"
        results["significance"] = significance<jupyter_output><empty_output><jupyter_text>Save the file and rerun your code so that you can call these statistics.<jupyter_code>reg.cov_matrix.to_csv("regCovMatrix.csv")
reg.cov_matrix<jupyter_output><empty_output><jupyter_text>Next, we have included the code for generating standard errors, t-stats, and p-values.<jupyter_code>reg.estimates.to_csv("regResults.csv")
reg.estimates<jupyter_output><empty_output><jupyter_text>$rootMSE = \sqrt{\sigma^2}$$R^2 = \frac{SSR}{SST}$ $ F =\frac{\frac{SST - SSE}{K - 1}}{\frac{SSE}{N - K}}$The F-statistic tests the likelihood of whether or not the values of our estimated parameters are all zero:$\beta_1 = \beta_2 = . . . = \beta_{n-1} = \beta_n$<jupyter_code>reg.stats_DF.to_csv("regStatistics.csv")
reg.stats_DF
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages

def plot_scatter_with_estimator(title, data, x_vars, y_var, pp = None):
    # set default font size
    plt.rcParams.update({"font.size":19})
    # use a loop to call each exogenous
    y = y_var[0]
    for x in x_vars:
        fig, ax = plt.subplots(figsize=(12,8))
        y_label1 = "Estimate"
        y_label2 = "Observation"
        data.plot.scatter(x=x, y=y + " estimator", ax = ax, c = "r",
                         s = 10, label = y_label1, legend = False)
        
        # erase the y_axis label so that "estimator" is not present
        plt.ylabel("")
        data.plot.scatter(x = x, y = y, ax = ax, s = 10, 
                          label = y_label2, legend = False)
        plt.legend(loc="upper left", labels = [y_label1, y_label2],
                  bbox_to_anchor = (0, 1.17))
        
        plt.title("Observations and\nOLS Estimates")
        # remove lines marking unit on the axis
        ax.xaxis.set_ticks_position("none")
        ax.yaxis.set_ticks_position("none")
        if pp != None: pp.savefig(fig, bbox_inches = "tight")
pp = PdfPages("RegressionPlots.pdf")
plot_scatter_with_estimator("Plot", reg.data, x_vars, y_var, pp)

### plot


pp.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11263259402121505, 'fraction_numerical': 0.016538090646094503, 'mean_word_length': 2.551541095890411, 'pattern_counts': {'":': 1, '<': 69, '<?xml version=': 0, '>': 68, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '239431', 'license': 'MIT', 'n_tokens_mistral': 7147, 'n_tokens_neox': 6817, 'n_words': 1820, 'path': 'In Class Projects/Old Project 5 - OLS (and Linear Algebra).ipynb', 'repo_name': 'hunterluepke/Learn-Python-for-Stats-and-Econ', 'size': '20740'}