[Kapitel 8] Ende des Kapitels Problem der metrischen Ökonomie (Yukikaku), Antwort von Python

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

Grundsätzlich mache ich Antworten, während ich mir die offiziellen Dokumente von Statistikmodellen und linearen Modellen ansehe, Github. Ich weiß überhaupt nichts über Python. Wenn Sie also Fehler haben, weisen Sie bitte darauf hin und fragen Sie.

Das Problem am Ende von Kapitel 8 ist die Reproduktion der Tabellen 8-1 und 8-5.

Beispiel für eine empirische Analyse, Reproduktion von Tabelle 8-1

Die Reproduktion von Tabelle 8-1 erfordert die Behandlung, wenn die erläuterte Variable eine binäre Variable ist. Wir werden das lineare Wahrscheinlichkeitsmodell, das Probit-Modell und das Logit-Modell als Mittel vergleichen.

Lineares Wahrscheinlichkeitsmodell


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


#Daten lesen und bestätigen
df = pd.read_stata('piaac.dta')
df.head()

In der Analyse in Tabelle 8-1 ist [lfs] auf Frauen beschränkt und kehrt zu [Bildung], [Alter], [Paar], [Kind] zurück. Holen Sie sich nur die Spalten und Datensätze, die Sie benötigen.


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

Das Problem hierbei ist, dass die erklärte Variable [lfs] eine qualitative Variable ist. Lassen Sie uns vorerst überprüfen, welchen Wert die erklärte Variable annimmt.


#['lfs']Überprüfen Sie den Typ von
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]

#Überprüfen Sie auch das Alter
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 scheint vier Arten von Werten anzunehmen: Angestellt, Arbeitslos, Nicht erwerbstätig, Nicht bekannt. Wie dies zu interpretieren ist, ist von Natur aus eine sehr heikle Angelegenheit, aber Tabelle 8-1 fällt nicht aus. Nicht bekannt und scheint als nicht arbeitende Person zu gelten. Was zu tun ist, wenn die erläuterte Variable nicht beachtet wird (Stichprobenauswahl wird zum Problem), wird bei der Implementierung des Hekit-Modells berücksichtigt. Auch die Arbeitslosen sind buchstäblich keine Arbeitskräfte, aber da alle Daten nur für Personen über 16 Jahre gelten, können sie von Natur aus arbeiten. Ich entscheide mich jedoch aus irgendeinem Grund, nicht zu arbeiten, zum Beispiel wegen niedriger angebotener Löhne. Hier betrachten wir den Fall, in dem die erläuterte Variable eine binäre Variable verwendet. Erstellen Sie daher eine Dummy-Variable mit Employed = 1 und other = 0 zur Analyse.


#Dummy variable Spalte['emp_dummy']Schaffung
df_1.loc[df_1['lfs']=='Employed','emp_dummy']=1
df_1.loc[df_1['lfs']!='Employed','emp_dummy']=0

# Überprüfen Sie den Datentyp für alle Fälle
df_1.dtypes
#lfs          category
educ         category
age          category
couple        float32
child         float32
emp_dummy     float64
dtype: object

Aus irgendeinem Grund sind Bildung und Alter kategorisch ... In eine numerische Variable konvertieren.

#Datentyp ändern
df_1[['educ','age']] = df_1[['educ','age']].apply(pd.to_numeric)

Ignorieren Sie zunächst die fehlenden Werte und führen Sie die Analyse durch.


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

Geben Sie bei Verwendung des weißen Standardfehlers für die Freiheitskorrektur HC1 für cov_type an. Geben Sie HC0 an, um den normalen weißen Standardfehler zu verwenden.


#Überprüfen Sie das Ergebnis
result_1_1.summary()

	coef	std err	t	P>t	[0.025	[0.025
Intercept	0.3759	0.105	3.566	0.000	0.169	0.583
educ	0.0195	0.006	3.312	0.001	0.008	0.031
age	0.0020	0.001	1.726	0.085	-0.000	0.004
couple	-0.1178	0.031	-3.743	0.000	-0.180	-0.056
child	0.0137	0.013	1.067	0.286	-0.012	0.039

Das Ergebnis kann als genau das gleiche bezeichnet werden. Tabelle 8-1 im Lehrbuch scheint die fehlenden Werte zu ignorieren und die Analyse durchzuführen.

Bonus: Berücksichtigen Sie fehlende Werte

Ich werde es auch analysieren, um mit fehlenden Werten umzugehen. Dies ist ein sehr heikles Thema, aber lassen Sie uns zunächst den Durchschnittswert eingeben.


#Füllen Sie fehlende Werte mit Durchschnittswerten aus
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())

#OLS-Schätzung
result_1_nonan = ols(formula = 'emp_dummy ~ educ + age + couple + child',
                     data = df_1).fit(cov_type='HC1',use_t=True)
#Überprüfen Sie das Ergebnis
result_1_nonan_summary()

	chef	std err	t	P>t	[0.025	0.975]
Intercept	0.0189	0.062	0.302	0.763	-0.104	0.141
educ	0.0450	0.004	10.951	0.000	0.037	0.053
age	0.0032	0.001	3.895	0.000	0.002	0.005
couple	-0.1035	0.022	-4.735	0.000	-0.146	-0.061
child	-0.0006	0.009	-0.059	0.953	-0.019	0.018

Der Wert hat sich etwas geändert. Schließlich scheinen bei der Analyse von Lehrbüchern alle Datensätze einschließlich fehlender Werte gelöscht und analysiert zu werden.

Probit-Modell

Folgen Sie von nun an dem Lehrbuch und ignorieren Sie alle fehlenden Werte.


#Datenformung
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)

Wenn Sie ein Probit-Modell verwenden, können Sie einfach "probit" anstelle von "ols" verwenden, wenn Sie Statistikmodelle nutzen.

###Probit-Modell
result_1_2 = probit(formula = 'emp_dummy ~ educ + age + couple + child',
                    data = df_2).fit()
#Überprüfen Sie das Ergebnis
result_1_2.summary().tables[1]

	coef	std err	z	P>z	[0.025	0.975]
Intercept	-0.3443	0.287	-1.198	0.231	-0.908	0.219
educ	0.0537	0.016	3.340	0.001	0.022	0.085
age	0.0053	0.003	1.789	0.074	-0.001	0.011
couple	-0.3391	0.097	-3.489	0.000	-0.530	-0.149
child	0.0397	0.032	1.256	0.209	-0.022	0.102

Dieses Mal ist das Ergebnis fast das gleiche wie in Tabelle 8-1.

Als nächstes folgt die Berechnung des Randeffekts. Wenn Sie Statistikmodelle verwenden, ist der Randeffekt des Probit-Modells mit der Methode "get_margeff" beendet.


#Randeffekt
result_1_2.get_margeff().summary().tables[1]

	dy/dx	std err	z	P>z	[0.025	0.975]
educ	0.0197	0.006	3.372	0.001	0.008	0.031
age	0.0019	0.001	1.794	0.073	-0.000	0.004
couple	-0.1243	0.035	-3.524	0.000	-0.193	-0.055
child	0.0145	0.012	1.258	0.209	-0.008	0.037

Dies ist das gleiche Ergebnis wie das Lehrbuch.

Logit-Modell

Genau wie bei Probit verwenden Sie einfach "logit". Der Randeffekt kann auf die gleiche Weise mit dem Modul "get_margeff" berechnet werden.


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

#Überprüfen Sie das Ergebnis
result_1_4.summary().tables[1]

#Randeffekt
result_1_4.get_margeff().summary().tables[1]

	coef	std err	z	P>z	[0.025	0.975]
Intercept	-0.5781	0.471	-1.227	0.220	-1.502	0.345
educ	0.0869	0.026	3.307	0.001	0.035	0.138
age	0.0088	0.005	1.806	0.071	-0.001	0.018
couple	-0.5587	0.163	-3.424	0.001	-0.879	-0.239
child	0.0716	0.058	1.242	0.214	-0.041	0.185

	dy/dx	std err	z	P>z	[0.025	0.975]
educ	0.0195	0.006	3.345	0.001	0.008	0.031
age	0.0020	0.001	1.812	0.070	-0.000	0.004
couple	-0.1255	0.036	-3.463	0.001	-0.197	-0.054
child	0.0161	0.013	1.244	0.213	-0.009	0.041

Bonus: Erstellen Sie ein Probit-Modell mit GenericLikelihoodModel.

statsmodels hat eine nützliche Klasse "GenericLikelihoodModel", die offiziell sagt 　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.

Es scheint praktisch zu sein, wenn Sie die maximale Wahrscheinlichkeit mithilfe einer hausgemachten Wahrscheinlichkeitsfunktion schätzen möchten.

Da das Beispiel des offiziellen Dokuments nur durch Ändern des Datensatzes implementiert wird, lesen Sie bitte auch das offizielle Dokument.


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

#Löschen Sie Datensätze, die fehlende Werte enthalten
df_4 = df_2.dropna()


#Erstellen Sie eine Matrix mit erklärenden Variablen.
exog = df_4.loc[:,['educ','age','couple','child']]
exog = sm.add_constant(exog,prepend=True)
#Erstellen Sie einen Vektor mit erklärten Variablen.
endog = df_4.loc[:,'emp_dummy']

Wenn die latente Variable durch die Formel auf der rechten Seite bestimmt wird, ist $ Y_ {i} ^ {*} = X_ {i} '\ beta + e_ {i} $ Die Log-Likelihood-Funktion des Binomialauswahlmodells ist im Allgemeinen

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

Alles, was Sie tun müssen, ist, F () in eine normalverteilte Verteilungsfunktion umzuwandeln und in diese einzutauchen.

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()#Protokollwahrscheinlichkeit des Probit-Modells

Es scheint, dass dies allein zur Schätzung herangezogen werden kann.


#Passend zum Probit-Modell
result_mine = MyProbit(endog, exog).fit(maxiter=1000)

#Überprüfen Sie das Ergebnis
result_mine.summary().tables[1]

	chef	stderr	z	P>z	[0.025	0.975]
const	-0.3444	0.287	-1.198	0.231	-0.908	0.219
educ	0.0537	0.016	3.340	0.001	0.022	0.085
age	0.0053	0.003	1.790	0.073	-0.001	0.011
couple	-0.3391	0.097	-3.489	0.000	-0.530	-0.149
child	0.0397	0.032	1.256	0.209	-0.022	0.102

Das Ergebnis ist das gleiche wie bei Verwendung des Probit-Moduls von Statistikmodellen. Es gibt kein offiziell zusammengeführtes Modul für Tobit-Modelle. Wenn Sie jedoch GenericLikelihoodModels verwenden, können Sie es anscheinend nur durch Erhöhen der Likelihood-Funktion erstellen.

Beispiel für eine empirische Analyse, Reproduktion von Tabelle 8-5

Obwohl es sich um ein Heckit-Modell handelt, hatten Statistikmodelle und lineare Modelle kein Modul, das eine zweistufige Schätzung und die wahrscheinlichste Schätzung durchführen kann, soweit ich das Dokument gelesen habe. Ich wäre Ihnen dankbar, wenn Sie es mich wissen lassen könnten. Obwohl nicht offiziell zusammengeführt, gab es in der Niederlassung ein Modul namens [Heckman].

Benutze das.


import Heckman


#Lesen und Formatieren von Daten

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'])

#Bestätigung der Daten
df_3

Einige Leute haben ihren Lohn aus irgendeinem Grund nicht eingehalten. Es ist möglich, dass Sie überhaupt nicht arbeiten, weil der angebotene Lohn geringer ist als der reservierte Lohn. Wenn die Löhne so wie sie sind an die erklärenden Variablen zurückgegeben werden, wird die Stichprobenauswahl wahrscheinlich zu einem Problem.

OLS-Schätzung

Wie auch immer, um die Tabelle zu reproduzieren, werde ich OLS so ausführen, wie es ist.


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

#Überprüfen Sie das Ergebnis
result_5_1.summary()

	coef	stderr	t	P>t	[0.025	0.975]
Intercept	5.8310	0.236	24.726	0.000	5.368	6.294
educ	0.0862	0.014	6.074	0.000	0.058	0.114
experience	0.0069	0.009	0.762	0.446	-0.011	0.025
experience_2	-0.0126	0.015	-0.820	0.413	-0.043	0.018

No. Observations:	984

Insgesamt gab es 1747 Proben, aber es wurden 984 Proben verwendet.

Zweistufiges Hekit


#Erstellen Sie einen Vektor mit erklärten Variablen.
endog = df_3['logwage']

#Erstellen Sie eine Matrix mit erklärenden Variablen für die Auswahlformel.
exog_select = df_3[['educ','experience','experience_2','couple','child']]
exog_select = sm.add_constant(exog_select,prepend=True)

#Erstellen Sie eine Matrix mit erklärenden Variablen für die Gleichung der zweiten Stufe.
exog = df_3[['educ','experience','experience_2']]
exog = sm.add_constant(exog,prepend=True)


#Heckmans zweistufige Schätzung
result_5_2 = Heckman.Heckman(endog,exog,exog_select).fit(method='twostep')

#Überprüfen Sie das Ergebnis
result_5_2.summary()

Höchstwahrscheinlich Methode


#Höchstwahrscheinlich Schätzung
result_5_4 = Heckman.Heckman(endog,exog,exog_select).fit(method='mle',
                                                         method_mle='nm',
                                                         maxiter_mle=2000)

#Überprüfen Sie das Ergebnis
result_5_4.summary()

Verweise

Bei der Vorbereitung dieser Antwort habe ich auf die folgenden Materialien Bezug genommen.

Yoshihiko Nishiyama, Mototsugu Shintani, Daiji Kawaguchi, Ryo Okui "Messökonomie", 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