%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import numpy as np
import pandas as pd
import statsmodels.api as sm

# Simulo la decisión de trabajo a partir de una log normal.
# Parámetros de la distribución lognormal bivariante
mean = [0, 0]  # Media para la distribución normal subyacente (ambos iguales)
var1 = 1  # Varianza para la primera variable
var2 = 4  # Varianza para la segunda variable
correlation = 0.75  # Correlación entre las dos variables

tamano_pobl = 10000

# Matriz de covarianza
cov = [[var1, correlation * np.sqrt(var1 * var2)],
       [correlation * np.sqrt(var1 * var2), var2]]

# Generar muestras de la normal bivariante
np.random.seed(42)  # Para reproducibilidad
normal_samples = np.random.multivariate_normal(mean, cov, size=tamano_pobl)

# Transformar a lognormal bivariante
lognormal_samples = np.exp(normal_samples)

# Separar las dos variables para análisis
e_r = np.log(lognormal_samples[:, 0]) # r: salario de reserva
e_m = np.log(lognormal_samples[:, 1]) # m: salario de mercado

# Recta de regresión
cov_x = [[.1, 0],
       [0,.1]]

x = np.random.multivariate_normal(mean, cov_x, size=tamano_pobl)

x_r = np.log(6) + x[:, 0]
x_m = np.log(6) + x[:, 1]

w_r = 1 + 2.5*x_r + e_r
w_m = 1 + 3.5*x_m + e_m

# Graficar los datos
plt.figure(figsize=(8, 6))
plt.scatter(w_r, w_m, alpha=0.3, s=5)
plt.title("Distribución conjunta  del salario de reserva y salario de mercado (logs)")
plt.xlabel(" Salario de reserva ($w_r$)")
plt.ylabel(" Salario de mercado ($w_m$)")
plt.grid(alpha=0.5)
plt.show()

# Grafico del salario de mercado (potencial) con respecto a la educacion  
plt.figure(figsize=(8, 6))
plt.scatter(x_m, w_m, alpha=0.3, s=5)
plt.title("Relación salarios potencial-educación (en logaritmos)")
plt.ylabel(" Salario de mercado (potencial) ($w_m$)")
plt.xlabel(" Educación $X$ ")
plt.grid(alpha=0.5)
plt.show()

# Recta de regresión verdadera del salario de mercado
wm_pred_verdadera = 1 + 3.5*x_m 

# Ordenar x_m y wm_pred según x_m
x_m_s, wm_pred_v_s = zip(*sorted(zip(x_m, wm_pred_verdadera)))

# Graficar los puntos de datos y la recta de regresión
plt.scatter(x_m, w_m, label="Datos originales", color="blue", s = 0.1)  # Puntos de datos
plt.plot(x_m_s, wm_pred_v_s, label="Recta de regresión verdadera", color="red")  # Recta de regresión

# Añadir la fórmula al gráfico
#Extraer coeficientes del modelo
intercepto = 1  # Coeficiente de la constante
pendiente = 3.5   # Coeficiente de la variable independiente
formula = f"$w_m = {intercepto:.2f} + {pendiente:.2f}X$"
plt.text(max(x_m) , np.mean(w_m) + 5, formula, fontsize=10, color="red", bbox=dict(facecolor='white', alpha=0.5))


plt.xlabel("Variable independiente (Educación $X_m$)")
plt.ylabel("Variable dependiente (Salarios de mercado (potencial) $w_m$)")
plt.title("Salarios de mercados (potencial) y educación")

plt.legend()
plt.grid()
plt.show()

#  Ajusto una regresión a la poblacion
X = x_m
X = sm.add_constant(X)  # Agregar una constante para el intercepto
y = w_m

modelo = sm.OLS(y, X).fit()

# Mostrar el resumen de los resultados
print(modelo.summary())

#Grafico
#Predecir los valores usando el modelo
wm_pred = modelo.predict(X)

# Ordenar x_0 y w0_pred según x_0
x_m_s, wm_pred_s = zip(*sorted(zip(x_m, wm_pred)))

# Graficar los puntos de datos y la recta de regresión
plt.scatter(x_m, w_m, label="Datos", color="blue", s = 0.1)  # Puntos de datos
plt.plot(x_m_s, wm_pred_s, label="Recta de regresión ajustada", color="black")  # Recta de regresión
plt.plot(x_m_s, wm_pred_v_s, label="Recta de regresión verdadera", color="red")  # Recta de regresión
# Añadir la fórmula al gráfico
#Extraer coeficientes del modelo
intercepto = modelo.params[0]  # Coeficiente de la constante
pendiente = modelo.params[1]   # Coeficiente de la variable independiente
formula = f"$w_m = {intercepto:.2f} + {pendiente:.2f}X$"
plt.text(min(x_m) , np.mean(w_m) - 6, formula, fontsize=10, color="black", bbox=dict(facecolor='white', alpha=0.5))

formula = f"$w_m = 1 + 3.5 X_m$"
plt.text(max(x_m) , np.mean(w_m) + 3, formula, fontsize=10, color="red", bbox=dict(facecolor='white', alpha=0.5))


plt.xlabel("Variable independiente (Educación $X_m$)")
plt.ylabel("Variable dependiente (Salarios potencial $w_m$)")
plt.title("Salarios de mercado (potencial) y educación: ajuste lineal en la poblacion")

plt.legend()
plt.grid()
plt.show()

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.226
Model:                            OLS   Adj. R-squared:                  0.226
Method:                 Least Squares   F-statistic:                     2921.
Date:                Fri, 13 Dec 2024   Prob (F-statistic):               0.00
Time:                        10:30:06   Log-Likelihood:                -21233.
No. Observations:               10000   AIC:                         4.247e+04
Df Residuals:                    9998   BIC:                         4.248e+04
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.0213      0.117      8.730      0.000       0.792       1.251
x1             3.4851      0.064     54.050      0.000       3.359       3.612
==============================================================================
Omnibus:                        2.205   Durbin-Watson:                   2.019
Prob(Omnibus):                  0.332   Jarque-Bera (JB):                2.184
Skew:                           0.024   Prob(JB):                        0.336
Kurtosis:                       3.054   Cond. No.                         13.6
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

# Vamos a extraer una muestra de la población 

muestra = np.random.uniform(0, 1, size = tamano_pobl) < 0.5

#  Definir las variables independientes y dependiente
X = x_m_m = x_m[muestra]
X = sm.add_constant(X)  # Agregar una constante para el intercepto
y = w_m_m = w_m[muestra]

# Ajustar el modelo de regresión lineal
modelo = sm.OLS(y, X).fit()

# Mostrar el resumen de los resultados
print(modelo.summary())
#Predecir los valores usando el modelo
wm_pred = modelo.predict(X)

# Ordenar x_0 y w0_pred según x_0
x_m_m_s, wm_pred_s = zip(*sorted(zip(x_m_m, wm_pred)))

wm_pred_v_s = 1 + 3.5 * np.array(x_m_m_s)

# Graficar los puntos de datos y la recta de regresión
plt.scatter(x_m, w_m, label="Danos de la muestra", color="blue", s = 0.1)  # Puntos de datos
plt.plot(x_m_m_s, wm_pred_s, label="Recta de regresión: muestra", color="black")  # Recta de regresión

plt.plot(x_m_m_s, wm_pred_v_s, label="Recta de regresión: verdadera", color="red")  # Recta de regresión
# Añadir la fórmula al gráfico
#Extraer coeficientes del modelo
intercepto = modelo.params[0]  # Coeficiente de la constante
pendiente = modelo.params[1]   # Coeficiente de la variable independiente
formula = f"$w_0 = {intercepto:.2f} + {pendiente:.2f}X$"
plt.text(min(x_m) , np.mean(w_m) - 5, formula, fontsize=10, color="black", bbox=dict(facecolor='white', alpha=0.5))


plt.xlabel("Variable independiente (Educación $X$)")
plt.ylabel("Variable dependiente (Salarios de mercado $w_m$)")
plt.title("Salarios de mercado y educación: muestra aleatoria")

plt.legend()
plt.grid()
plt.show()

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.225
Model:                            OLS   Adj. R-squared:                  0.225
Method:                 Least Squares   F-statistic:                     1442.
Date:                Fri, 13 Dec 2024   Prob (F-statistic):          2.22e-277
Time:                        10:30:07   Log-Likelihood:                -10561.
No. Observations:                4977   AIC:                         2.113e+04
Df Residuals:                    4975   BIC:                         2.114e+04
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.1059      0.165      6.709      0.000       0.783       1.429
x1             3.4397      0.091     37.979      0.000       3.262       3.617
==============================================================================
Omnibus:                        3.026   Durbin-Watson:                   2.025
Prob(Omnibus):                  0.220   Jarque-Bera (JB):                3.174
Skew:                           0.012   Prob(JB):                        0.205
Kurtosis:                       3.121   Cond. No.                         13.6
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

# Trabaja: solo si salario de mercado es mayor que el de reserva
# Distribución de los salarios de reserva y de los salarios de mercado (potencial)

# Plot distributions
plt.figure(figsize=(10, 6))

# Distribution of observed outcomes for the treated
sns.histplot(w_r, color="red", label="Salario reserva", kde=True, stat="density", alpha=0.35, bins=50)

# Distribution of potential outcomes for the treated
sns.histplot(w_m, color="blue", label="Salario mercado", kde=True, stat="density", alpha=0.25, bins=50)

# Customize the plot
plt.title("Distribuciones marginales de salario reserva y mercado (potencial)", fontsize=16)
plt.xlabel("Salario", fontsize=14)
plt.ylabel("Density", fontsize=14)
plt.legend(fontsize=12)
plt.grid(False)  # Disable the grid
plt.gca().spines['top'].set_visible(False)  # Remove the top border
plt.gca().spines['right'].set_visible(False)  # Remove the right border

# Show the plot
plt.tight_layout()
plt.show()

# Decisión de trabajar 
I = w_m > w_r + 1
np.mean(I)

0.6553

# Plot distributions
plt.figure(figsize=(10, 6))

# Distribution of observed outcomes for the treated
sns.histplot(w_m[I], color="red", label="$Observado$", kde=True, stat="density", alpha=0.35, bins=50)

# Distribution of potential outcomes for the treated
sns.histplot(w_m, color="blue", label="Potential", kde=True, stat="density", alpha=0.25, bins=50)

# Customize the plot
plt.title("Distribución del salario mercado: potencial y observado", fontsize=16)
plt.xlabel("Salarios", fontsize=14)
plt.ylabel("Density", fontsize=14)
plt.legend(fontsize=12)
plt.grid(False)  # Disable the grid
plt.gca().spines['top'].set_visible(False)  # Remove the top border
plt.gca().spines['right'].set_visible(False)  # Remove the right border

# Show the plot
plt.tight_layout()
plt.show()

# Impacto en la regresión

w_ms = w_m[I] # qué estoy haciendo aqui
x_ms = x_m[I]
plt.figure(figsize=(13, 10))
plt.xlabel("schooling")
plt.ylabel("w0")
plt.plot(np.unique(x_m),
         np.poly1d(np.polyfit(x_m, w_m, 1))(np.unique(x_m)),
         color='darkblue', linewidth=1.5, label =" Regresión verdadera: datos de mercado potenciales")
plt.scatter(x_m, w_m, alpha = 0.5)

plt.plot(np.unique(x_ms),
         np.poly1d(np.polyfit(x_ms, w_ms, 1))(np.unique(x_ms)),
         color='red',linewidth=1.5, label =" Regresión datos observados: $w_m > w_r$")
plt.scatter(x_ms, w_ms, facecolor = 'none', s = 300, color = "brown")

plt.title("Ajuste en salarios de los que trabajan y potenciales", fontsize=16)
plt.xlabel("Salarios", fontsize=14)
plt.ylabel("Density", fontsize=14)
plt.legend()

plt.show()

#  Ajusto una regresión a la poblacion
X = x_ms
X = sm.add_constant(X)  # Agregar una constante para el intercepto
y = w_ms

modelo = sm.OLS(y, X).fit()

# Mostrar el resumen de los resultados
print(modelo.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.103
Model:                            OLS   Adj. R-squared:                  0.103
Method:                 Least Squares   F-statistic:                     754.5
Date:                Fri, 13 Dec 2024   Prob (F-statistic):          2.55e-157
Time:                        10:30:12   Log-Likelihood:                -12875.
No. Observations:                6553   AIC:                         2.575e+04
Df Residuals:                    6551   BIC:                         2.577e+04
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.4907      0.141     31.807      0.000       4.214       4.768
x1             2.0344      0.074     27.469      0.000       1.889       2.180
==============================================================================
Omnibus:                       70.566   Durbin-Watson:                   2.010
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               74.597
Skew:                           0.233   Prob(JB):                     6.33e-17
Kurtosis:                       3.236   Cond. No.                         16.0
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

# Probabilidad de que un individuo no trabaje condisional a su salario de reserva

dens_c = sm.nonparametric.KDEMultivariateConditional(endog=[w_m],
                                                    exog=[w_r], dep_type='c', indep_type='c', bw='normal_reference')


c = 1- dens_c.cdf([w_m - 1 ])
plt.scatter(w_m, c , s = 0.5)
plt.plot(np.unique(w_m),
         np.poly1d(np.polyfit(w_m, c, 3))(np.unique(w_m)),
         color='red', linewidth=1.5)
plt.xlabel("Salario de mercado")
plt.ylabel("Pobabilidad")
plt.title("Probabilidad de no trabajar condicional al salario\n de mercado")
plt.xlim(2.5, 12.5)


plt.show()