[Chapitre 8] Problème de fin de chapitre de l'économie métrique (Yukikaku), réponse de python

[Measurement Economics](https://www.amazon.co.jp/ Quantitative Economics-New-Liberal-Arts-Selection / dp / 4641053855 / ref = asc_df_4641053855_nodl /? Tag = jpgo-22 & linkCode = df0 & hvadid = 342654745883 & hvpos = & hvn = g & hvrand = 3840609228244811524 & hvpone = & hvptwo = & hvqmt = & hvdev = m & hvdvcmdl = & hvlocint = & hvlocphy = 19009363 & hvtargid = pla-796815219068 & psc = 1 & th = 1) & psc = 1)

En gros, je fais des réponses en regardant les documents officiels des modèles de statistiques et des modèles linéaires, github. Je ne connais pas du tout python, donc si vous avez des erreurs, veuillez le signaler et demander.

Le problème à la fin du chapitre 8 est la reproduction des tableaux 8-1 et 8-5.

Exemple d'analyse empirique, reproduction du tableau 8-1

La reproduction du Tableau 8-1 nécessite le traitement lorsque la variable expliquée est une variable binaire. Nous comparerons le modèle de probabilité linéaire, le modèle probit et le modèle logit comme moyennes.

Modèle de probabilité linéaire

import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
sns.set()

import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.formula.api import logit, probit

#Lire et confirmer les données
df = pd.read_stata('piaac.dta')
df.head()

Dans l'analyse du tableau 8-1, [lfs] est limité aux femmes et est retourné à [éduc], [âge], [couple], [enfant]. Obtenez uniquement les colonnes et les enregistrements dont vous avez besoin.

df_1 = df[df['gender']== 'Female']
df_1 = df_1[['lfs','educ','age','couple','child']]

Le problème ici est que la variable expliquée [lfs] est une variable qualitative. Pour le moment, vérifions quelle valeur prend la variable expliquée.

#['lfs']Vérifiez le type de
df_1['lfs'].unique()
#[Employed, Out of the labour force, Unemployed, Not known]
Categories (4, object): [Employed < Unemployed < Out of the labour force < Not known]

#Vérifiez également l'âge
df_1['age'].unique()
#[45, 29, 35, 48, 60, ..., 41, 37, 19, 58, 28]
Length: 50
Categories (50, int64): [16 < 17 < 18 < 19 ... 62 < 63 < 64 < 65]

Lfs semble prendre quatre types de valeurs: Employé, Chômeur, Hors de la main-d'œuvre, Inconnu. Comment interpréter cela est une question très délicate par nature, mais le tableau 8-1 ne baisse pas. Inconnu et semble être considéré comme des personnes non actives. Que faire si la variable expliquée n'est pas observée (la sélection de l'échantillon devient un problème) sera prise en compte lors de la mise en œuvre du modèle Hekit. De plus, la population hors du marché du travail est littéralement une population non active, mais comme toutes les données ne concernent que les personnes de plus de 16 ans, elles peuvent travailler par nature. Cependant, je choisis de ne pas travailler pour une raison quelconque, comme les bas salaires offerts. Ici, nous allons considérer le cas où la variable expliquée prend une variable binaire, donc créez une variable fictive avec Employed = 1 et autre = 0 pour l'analyse.

#Colonne de variable factice['emp_dummy']Création
df_1.loc[df_1['lfs']=='Employed','emp_dummy']=1
df_1.loc[df_1['lfs']!='Employed','emp_dummy']=0

# Vérifiez le type de données au cas où
df_1.dtypes
#lfs          category
educ         category
age          category
couple        float32
child         float32
emp_dummy     float64
dtype: object

Pour une raison quelconque, l'éducation et l'âge sont catégoriques ... Convertissez en une variable numérique.

#Changer le type de données
df_1[['educ','age']] = df_1[['educ','age']].apply(pd.to_numeric)

Tout d'abord, ignorez les valeurs manquantes et effectuez l'analyse.

#Estimation OLS
result_1_1 = ols(formula = 'emp_dummy ~ educ + age + couple + child',
                data = df_1).fit(cov_type='HC1',use_t=True)

Lorsque vous utilisez l'erreur standard de blanc de correction de liberté, spécifiez HC1 pour cov_type. Spécifiez HC0 pour utiliser l'erreur standard du blanc normal.

#Vérifiez le résultat
result_1_1.summary()

On peut dire que le résultat est exactement le même. Il semble que le tableau 8-1 du manuel ignore les valeurs manquantes pour l'analyse.

Bonus: considérez les valeurs manquantes

Je vais également l'analyser pour traiter les valeurs manquantes. C'est une question très délicate, mais pour l'instant, mettons la valeur moyenne.

#Remplissez les valeurs manquantes avec des valeurs moyennes
df_1['educ'] = df_1['educ'].fillna(df_1['educ'].mean())
df_1['age'] = df_1['age'].fillna(df_1['age'].mean())
df_1['couple'] = df_1['couple'].fillna(df_1['couple'].mean())
df_1['child'] = df_1['child'].fillna(df_1['couple'].mean())

#Estimation OLS
result_1_nonan = ols(formula = 'emp_dummy ~ educ + age + couple + child',
                     data = df_1).fit(cov_type='HC1',use_t=True)
#Vérifiez le résultat
result_1_nonan_summary()

La valeur a un peu changé. Après tout, dans l'analyse des manuels, il semble que tous les enregistrements, y compris les valeurs manquantes, soient supprimés et analysés.

Modèle Probit

À partir de là, suivez le manuel et ignorez toutes les valeurs manquantes.

#Mise en forme des données
df = pd.read_stata('piaac.dta')
df_2 = df[df['gender']== 'Female']
df_2 = df_2[['lfs','educ','age','couple','child']]
df_2.loc[df_2['lfs']=='Employed','emp_dummy']=1
df_2.loc[df_2['lfs']!='Employed','emp_dummy']=0
df_2[['educ','age']] = df_2[['educ','age']].apply(pd.to_numeric)

Si vous utilisez un modèle probit, vous pouvez simplement utiliser "probit" au lieu de "ols", si vous profitez des statsmodels.

###modèle probit
result_1_2 = probit(formula = 'emp_dummy ~ educ + age + couple + child',
                    data = df_2).fit()
#Vérifiez le résultat
result_1_2.summary().tables[1]

Cette fois, le résultat est presque le même que celui du tableau 8-1.

Vient ensuite le calcul de l'effet marginal. Si vous utilisez statsmodels, l'effet marginal du modèle Probit est terminé avec la méthode "get_margeff".

#Effet marginal
result_1_2.get_margeff().summary().tables[1]

C'est le même résultat que le manuel.

Modèle Logit

Tout comme avec Probit, vous utilisez simplement "logit". L'effet marginal peut être calculé de la même manière avec le module "get_margeff".

#Modèle Logit
result_1_4 = logit(formula = 'emp_dummy ~ educ + age + couple + child',
                    data = df_2).fit()

#Vérifiez le résultat
result_1_4.summary().tables[1]

#Effet marginal
result_1_4.get_margeff().summary().tables[1]

Bonus: Créez un modèle Probit en utilisant GenericLikelihoodModel.

statsmodels a une classe utile "GenericLikelihoodModel", qui dit officiellement 　The GenericLikelihoodModel class eases the process by providing tools such as automatic numeric differentiation and a unified interface to scipy optimization functions. Using statsmodels, users can fit new MLE models simply by “plugging-in” a log-likelihood function.

Il semble que ce soit pratique lorsque vous souhaitez estimer le maximum de vraisemblance à l'aide d'une fonction de vraisemblance maison.

Étant donné que l'exemple du document officiel n'est mis en œuvre qu'en modifiant l'ensemble de données, veuillez également vous reporter au document officiel.

#importer
from scipy import stats
from statsmodels.base.model import GenericLikelihoodModel

#Supprimer les enregistrements contenant des valeurs manquantes
df_4 = df_2.dropna()

#Créez une matrice de variables explicatives.
exog = df_4.loc[:,['educ','age','couple','child']]
exog = sm.add_constant(exog,prepend=True)
#Créez un vecteur de variables expliquées.
endog = df_4.loc[:,'emp_dummy']

Si la variable latente est déterminée par la formule de droite, $ Y_ {i} ^ {*} = X_ {i} '\ beta + e_ {i} $ Modèle de sélection binaire La fonction de vraisemblance générale de type log est

logL(\beta)=\sum_{i=1}^{n}[Y_{i}logF(X_{i}'\beta)+(1-Y_{i})log(1-F(X'_{i}\beta)]

Donc, tout ce que vous avez à faire est de changer F () en une fonction de distribution normalement distribuée et de vous y plonger.

class MyProbit(GenericLikelihoodModel):
    def loglike(self,params):
        exog = self.exog
        endog = self.endog
        q= 2*endog-1
        return stats.norm.logcdf(q*np.dot(exog,params)).sum()#Log de la probabilité du modèle Probit

Il semble que vous puissiez estimer avec juste cela.

#Adapté au modèle Probit
result_mine = MyProbit(endog, exog).fit(maxiter=1000)

#Vérifiez le résultat
result_mine.summary().tables[1]

Le résultat est le même que lors de l'utilisation du module probit de statsmodels. Il n'y a pas de module officiellement fusionné pour les modèles Tobit, mais si vous utilisez GenericLikelihoodModels, il semble que vous puissiez le créer simplement en augmentant la fonction de vraisemblance.

Exemple d'analyse empirique, reproduction du tableau 8-5

Bien qu'il s'agisse d'un modèle Heckit, statsmodels et linearmodels n'avaient pas de module capable d'effectuer l'estimation en deux étapes et l'estimation la plus probable en l'état, pour autant que j'ai lu le document. Je vous serais reconnaissant si vous pouviez me le faire savoir. Cependant, bien que pas officiellement fusionné, il y avait un module appelé [Heckman] dans la branche.

Utilisez ceci.

import Heckman 

#Lecture et formatage des données

read_stata('piaac.dta')
df_3 = df[df['gender']== 'Female']
df_3 = df_3[['educ','couple','child','age']]
df_3 = df_3.dropna()
df_3['wage'] = df['wage']
df_3[['educ','age']] = df_3[['educ','age']].apply(pd.to_numeric)
df_3['experience'] = df_3['age']-df_3['educ']-6
df_3['experience_2'] = (df_3['experience']**2)/100
df_3['logwage'] = np.log(df_3['wage'])

#Confirmation des données
df_3

Certaines personnes n'ont pas respecté leur salaire pour une raison quelconque. Il est possible que vous ne travailliez pas en premier lieu parce que le salaire offert est inférieur au salaire réservé. Si les salaires sont renvoyés aux variables explicatives telles quelles, la sélection de l'échantillon risque de devenir un problème.

Estimation OLS

Quoi qu'il en soit, afin de reproduire le tableau, je vais exécuter OLS tel quel.

#Estimation OLS
result_5_1 = ols(formula = 'logwage ~ educ + experience + experience_2',
                 data = df_3).fit(cov_type='HC1',use_t=True)

#Vérifiez le résultat
result_5_1.summary()

Il y avait 1747 échantillons au total, mais le nombre d'échantillons utilisés est de 984.

Hekit en deux étapes

#Créez un vecteur de variables expliquées.
endog = df_3['logwage']

#Créez une matrice de variables explicatives pour la formule de sélection.
exog_select = df_3[['educ','experience','experience_2','couple','child']]
exog_select = sm.add_constant(exog_select,prepend=True)

#Créez une matrice de variables explicatives pour l'équation de deuxième étape.
exog = df_3[['educ','experience','experience_2']]
exog = sm.add_constant(exog,prepend=True)

#Estimation en deux étapes de Heckman
result_5_2 = Heckman.Heckman(endog,exog,exog_select).fit(method='twostep')

#Vérifiez le résultat
result_5_2.summary()

Méthode la plus probable

#Estimation la plus probable
result_5_4 = Heckman.Heckman(endog,exog,exog_select).fit(method='mle',
                                                         method_mle='nm',
                                                         maxiter_mle=2000)

#Vérifiez le résultat
result_5_4.summary()

Les références

En préparant cette réponse, je me suis référé aux matériaux suivants.

Yoshihiko Nishiyama, Mototsugu Shintani, Daiji Kawaguchi, Ryo Okui "Measurement Economics", Yukaku, 2019

Burce E.Hansen(2020)"ECONOMETRICS"

Mackinnon and White(1985)"Some Heteroskedasticity Consistent Covariance Matrix Estimators with Improved Finite Sample Properties",Journal of Econometrics, 1985, vol. 29, issue 3, 305-325

statmodels

linearmodels

QuantEcon