Grupo_9_py.ipynb

Kernel: Python 3 (ipykernel)

import pandas as pd
import numpy as np
import scipy.stats as stats
from scipy.stats import t
import pyreadr

cps2012_env = pyreadr.read_r("../data/cps2012.RData")
type(cps2012_env)  #diccionario

collections.OrderedDict

cps2012 = cps2012_env['data'] #se puede leer en dataframe ahora
cps2012.drop(['year'], axis=1, inplace = True )  #dropear columna de año 2012 para todos los valores

cps2012['married'] = cps2012['married'].astype('int')   #convertir bool en enteros
cps2012['ne'] = cps2012['ne'].astype('int')
cps2012['sc'] = cps2012['sc'].astype('int')

cps2012=cps2012.head(2000)
cps2012

Valores a evaluar:

#Se pide mostrar DataFrame de variables explicativas
X = cps2012.iloc[ : , 1: ]  
#X.shape

#Se pide lista de valores para seleccionar 'x' de la Dataframe de variables explicativas
lista = [ 1, 5, 17,20]   #lista con valores de 0-20
lista

#Se pide el vector de la variable Y
y = cps2012.iloc[:, 0]

 class RegClass( object ):  # == RegClass():
    
    def __init__( self, X : pd.DataFrame , y : pd.Series , lista: list, robust = True ):
    
        if not isinstance( X, pd.DataFrame ):
            raise TypeError( "X must be a pd.DataFrame." )

        if not isinstance( y , pd.Series ):
            raise TypeError( "y must be a pd.Series." )
            
        # asignando los cuatro atributos de la clase:
        self.X = X
        self.y = y 
        self.lista = lista
        self.robust = robust
        
        # agregar columna de unos en self.X:
        self.X[ 'Intercept' ] = 1
        cols = self.X.columns.tolist()
        new_cols_orders = [cols[ -1 ]] + cols[ 0:-1 ] 
        self.X = self.X.loc[ : , new_cols_orders ] 
        
        
        # agregar variables de x seleccionadas, y y nombres de x
        self.X_np=self.X.iloc[:,lista].values
        self.y_np=y.values.reshape(-1,1)
        self.columns=self.X.iloc[:,lista].columns.tolist()
        
    def reg_beta_OLS( self ):
        
        # reasignar nombres a atributos    
        X_np=self.X_np
        y_np=self.y_np
            
        # beta_ols
        beta_ols = np.linalg.inv( X_np.T @ X_np ) @ ( X_np.T @ y_np )
        self.beta_ols=beta_ols
            
        # columnas de X
        index_names = self.columns
            
        # Output
        beta_OLS_output = pd.DataFrame( beta_ols , index = index_names , columns = [ 'Coef.' ] )
                 
        # agregar Dataframe de coeficientes como atributo 
        self.beta_OLS_output = beta_OLS_output
        
        return beta_OLS_output
           
    def var_standar(self): 
     
        #llamar método anterior
        self.reg_beta_OLS()  
                
        #reasignar nombres a atributos del método anterior
        X_np = self.X_np
        y_np = self.y_np
        index_names = self.columns
        beta_OLS = self.beta_OLS_output.values.reshape(- 1, 1) # convierte dataframe en vector

        #### Operaciones ####
        e = y_np - ( X_np @ beta_OLS )             # vector de errores: error_i = Y_i - Y_estimado_i
        self.e=e

        N = X.shape[ 0 ]                           # numero de filas
        k = X.shape[ 1 ]                           # numero de columnas
        ee=(e.T @ e)[ 0 ]                          # sumatoria e^2
        self.ee=ee
        error_var = (ee)/( N - k )                 # s^2= e^2/(n-k)

        ### 1. matriz de varianza y cov estándar  ###
        var_OLS = error_var*np.linalg.inv( X_np.T @ X_np ) # var_OLS(betas) = s^2 * (X'X)-1

        ### 2. error estandar de cada coeficiente ###
        sd = np.sqrt( np.diag(var_OLS) )           # desv(betas)=diagonal_var(betas)^(1/2)
        self.sd=sd
                
        ### 3. intervalos de confianza ###
        superior= beta_OLS.ravel() + 1.96 * sd
        self.superior=superior
        inferior= beta_OLS.ravel() - 1.96 * sd
        self.inferior=inferior

        
    def var_robust(self):
        
        #llamar a método anterior
        self.reg_beta_OLS()  
                
        #y reasignar nombres a atributos del método anterior
        X_np = self.X_np
        y_np = self.y_np
        index_names = self.columns
        beta_OLS = self.beta_OLS_output.values.reshape(- 1, 1) # convertir dataframe a vector

        #### Operaciones ####
        e = y_np - ( X_np @ beta_OLS )             # vector de errores: error_i = Y_i - Y_estimado_i
        self.e=e

        N = X.shape[ 0 ]                           # numero de filas
        k = X.shape[ 1 ]                           # numero de columnas
        ee=(e.T @ e)[ 0 ]                          # sumatoria e^2
        self.ee=ee

        self.w=np.eye(N)*ee
        
        ### 1. matriz de varianza y cov robusta  ###
        var_robust=np.linalg.inv(X_np.T @ X_np) @ X_np.T @(np.eye(N)*ee) @ X_np @ np.linalg.inv(X_np.T@X_np) 
        self.vr=var_robust                         # var_OLS(betas) = (X'X)-1* X' * White * X * (X'X)-1
                
        ### 2. error estandar de cada coeficiente ###
        sd = np.sqrt( np.diag(var_robust) )        # desv(betas)=var(betas)^(1/2)
        self.sd=sd
                
        ### 3. intervalos de confianza ###
        superior= beta_OLS.ravel() + 1.96 * sd
        self.superior=superior
        inferior= beta_OLS.ravel() - 1.96 * sd
        self.inferior=inferior

    def reg_estad(self):
        
        #llamar a método anterior
        self.reg_beta_OLS()
        
        N = self.X_np.shape[ 0 ]                    # Numero de filas
        
        y_est=self.X_np @ self.beta_ols             # fila de Y_estimados
        
        self.SCR = np.sum(np.square(self.y_np - y_est ))           # Sumatoria Cuadrado de Residuos
        
        self.SCT = np.sum(np.square(self.y_np - np.mean(y_est)))   # Sumatoria Cuadrado Totales      
           
        #### Root-MSE ####  
        root_mse = np.sqrt(self.SCT/ N)
        self.root_mse=root_mse
    
        #### R cuadrado ####  *****
        R2 = 1- self.SCR/self.SCT
        self.R2=R2
        
        fit = {"Root_MSE":root_mse, "R2": R2}     
        
        return fit
           
    def reg_OLS(self):
        
        # a. Se corren las funciones y ordenan valores para colocarlos en una tabla:
        self.reg_estad()
        self.reg_beta_OLS() 
        
             # a.1. coeficiente estimados
        beta_ols = self.beta_ols
        
        # b. Varianza de acuerdo a errores estandar o robustos
        
            # b.1. errores estandar y límites
        self.var_standar()
        
        sd_standar=self.sd.reshape( -1, 1 )
        
        superior_standar=self.superior.reshape( -1, 1 )
        inferior_standar=self.inferior.reshape( -1, 1 )
        
            # b.2. errores robustos y límites
        self.var_robust()
             
        sd_robust=self.sd.reshape( -1, 1 )
        
        superior_robust=self.superior.reshape( -1, 1 )
        inferior_robust=self.inferior.reshape( -1, 1 )
        
        # c. Colocar valores en una tabla de acuerdo a errores estandar o robustos
        
        if self.robust:      
            table_data ={  "Coef."    : beta_ols.ravel() ,
                       "Std.Err." : sd_standar.ravel(),
                       "Interv. sup."  : superior_standar.ravel(),
                       "Interv. inf."  : inferior_standar.ravel() }
        
        # d. crear dataframe 
            reg_OLS = pd.DataFrame( table_data , index = self.columns )
             
        # e. crear diccionario
            dic = {"OLS": reg_OLS, "root MSE": self.root_mse, "R2": self.R2}
            
        else:
            table_data ={  "Coef."    : beta_ols.ravel() ,
                       "Std.Err." : sd_robust.ravel(),
                       "Interv. sup."  : superior_robust.ravel(),
                       "Interv. inf."  : inferior_robust.ravel() }
        
        # d. crear dataframe 
            reg_OLS = pd.DataFrame( table_data , index = self.columns )
             
        # e. crear diccionario
            dic = {"OLS": reg_OLS, "root MSE": self.root_mse, "R2": self.R2}
            
        return dic      
   

Demostración:

A = RegClass( X, y ,lista, robust=True)   # correr datos

• Usando funcion general

A.X_np

array([[1.00000e+00, 0.00000e+00, 2.34256e+01, 1.00000e+00],
       [1.00000e+00, 0.00000e+00, 8.10000e+01, 1.00000e+00],
       [0.00000e+00, 0.00000e+00, 1.30321e+01, 1.00000e+00],
       ...,
       [1.00000e+00, 0.00000e+00, 1.30321e+01, 1.00000e+00],
       [1.00000e+00, 1.00000e+00, 1.00000e-04, 1.00000e+00],
       [0.00000e+00, 1.00000e+00, 2.56000e-02, 1.00000e+00]])

• Usando la funcion reg_beta_OLS para hallar los coeficientes;

$$ \begin{aligned}

\widehat{\beta} = (\widehat{\beta_1}, \widehat{\beta_2}, \widehat{\beta_3}, ..., \widehat{\beta_k}) \end{aligned} $$

A.reg_beta_OLS()

• Usando la funcion var_standar

A.var_standar()
A.superior     # límite superior 

array([-1.44421809e-01, -1.63013264e-01, -2.36178048e-04,  2.93652786e+00])

• Usando la funcion var_robust

A.var_robust()
A.superior     # límite superior 

array([2.30371053, 3.11630944, 0.03327363, 5.03995024])

• Usando la funcion reg_OLS para ver el diccionario que incluye:

$$ \begin{aligned}

\widehat{\beta} = (\widehat{\beta_1}, \widehat{\beta_2}, \widehat{\beta_3}, ..., \widehat{\beta_k}) \end{aligned} $$

$$Root_{MSE} = \sqrt{ \frac{ \sum(Y_i - \hat{Y_i})^2 }{n} }$$

$$R_{Cuadrado} = {1- \frac{SCR }{SCT} }$$

A = RegClass( X, y ,lista, robust=True)
A.reg_OLS()    # diccionario cuando hay evaluación con errores robustos

{'OLS':                  Coef.  Std.Err.  Interv. sup.  Interv. inf.
 female       -0.200733  0.028730     -0.144422     -0.257045
 nevermarried -0.238444  0.038485     -0.163013     -0.313874
 exp4         -0.001007  0.000393     -0.000236     -0.001778
 ne            2.888145  0.024685      2.936528      2.839763,
 'root MSE': 0.6478257987408599,
 'R2': 0.0428913127155941}

A.reg_OLS()['OLS']   #ver primer elemento del diccionario que es un dataframe

A = RegClass( X, y ,lista,robust=False)
A.reg_OLS()    # diccionario cuando hay evaluación estandar

{'OLS':                  Coef.  Std.Err.  Interv. sup.  Interv. inf.
 female       -0.200733  1.277778      2.303711     -2.705177
 nevermarried -0.238444  1.711609      3.116309     -3.593197
 exp4         -0.001007  0.017490      0.033274     -0.035288
 ne            2.888145  1.097860      5.039950      0.736340,
 'root MSE': 0.6478257987408599,
 'R2': 0.0428913127155941}