Единственный способ, которым я могу думать, это сделать вручную.
Образцы данных
import pandas as pd
import numpy as np
import statsmodels.api as sm
from linearmodels.datasets import wage_panel
from linearmodels.panel import PanelOLS
data = wage_panel.load()
Сначала мы пойдем по стопам Статы, создавая макеты для каждого изгод фиксировал эффекты, и мы пропускаем первое значение, лексикографически отсортированное (выполняется с аргументом drop_first=True
).Важно использовать np.unique
, чтобы затем получить метки, поскольку это тоже сортирует.Нет необходимости в statsmodels
для добавления константы, просто сделайте это сами:
data = wage_panel.load()
data = pd.concat([data, pd.get_dummies(data.year, drop_first=True)], axis=1)
exog_vars = ['exper','union','married'] + np.unique(data.year)[1::].tolist()
data = data.set_index(['nr', 'year'])
exog = data[exog_vars].assign(constant=1)
Если мы попытаемся запустить эту модель, она потерпит неудачу из-за идеальной коллинеарности.Поскольку мы делаем в рамках -регрессии, мы не можем просто проверить коллинеарность на экзоге, нам нужно сначала деинвестировать в панели, так как коллинеарность подразумеваемых независимых переменных имеет значение.Я сделаю копию здесь, чтобы не связываться с exog
, который мы будем использовать в окончательной регрессии.
exog2 = exog.copy()
exog2 = exog2 - exog2.groupby(level=0).transform('mean') + exog2.mean()
Теперь мы можем видеть, что существует совершенная коллинеарность;когда мы регрессируем в пределах де-значения переменной exog против каждой другой переменной, мы получаем идеальный R-квадрат для некоторых регрессий:
for col in exog2:
print(col, sm.OLS(exog2[col], exog2.drop(columns=col)).fit().rsquared)
#exper 1.0
#union 0.004326532427527674
#married 0.18901677578724896
#1981 1.0
#1982 1.0
#1983 1.0
#1984 1.0
#1985 1.0
#1986 1.0
#1987 1.0
Теперь, как Stata или другие программные пакеты решают, какую переменнуюуронить для меня загадка.В этом случае вы, вероятно, предпочтете сбросить один из своих манекенов над другими переменными.Давайте просто выберем 1987, чтобы мы могли получить тот же результат, что и Stata в конце.
exog2 = exog2.drop(columns=1987)
for col in exog2:
print(col, sm.OLS(exog2[col], exog2.drop(columns=col)).fit().rsquared)
#exper 0.48631
#union 0.0043265
#married 0.1890
#1981 0.34978
#1982 0.28369
#1983 0.2478680
#1984 0.2469180
#1985 0.2846552
#1986 0.35067075
#constant 0.94641
Итак, мы рассмотрели коллинеарность и можем вернуться к модели.Поскольку мы включили ежегодные фиксированные эффекты вручную, мы можем удалить time_effects
из модели.
mod = PanelOLS(data.lwage, exog.drop(columns=1987), entity_effects=True, time_effects=False)
print(mod.fit())
PanelOLS Estimation Summary
================================================================================
Dep. Variable: lwage R-squared: 0.1689
Estimator: PanelOLS R-squared (Between): -0.0882
No. Observations: 4360 R-squared (Within): 0.1689
Date: Sat, Mar 09 2019 R-squared (Overall): 0.0308
Time: 00:59:14 Log-likelihood -1355.7
Cov. Estimator: Unadjusted
F-statistic: 85.946
Entities: 545 P-value 0.0000
Avg Obs: 8.0000 Distribution: F(9,3806)
Min Obs: 8.0000
Max Obs: 8.0000 F-statistic (robust): 85.946
P-value 0.0000
Time periods: 8 Distribution: F(9,3806)
Avg Obs: 545.00
Min Obs: 545.00
Max Obs: 545.00
Parameter Estimates
==============================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
------------------------------------------------------------------------------
exper 0.0639 0.0033 19.593 0.0000 0.0575 0.0703
union 0.0834 0.0194 4.2887 0.0000 0.0453 0.1215
married 0.0583 0.0184 3.1759 0.0015 0.0223 0.0944
1981 0.0497 0.0201 2.4755 0.0133 0.0103 0.0890
1982 0.0399 0.0191 2.0888 0.0368 0.0025 0.0774
1983 0.0194 0.0187 1.0369 0.2998 -0.0172 0.0559
1984 0.0230 0.0187 1.2309 0.2184 -0.0136 0.0595
1985 0.0081 0.0191 0.4259 0.6702 -0.0294 0.0457
1986 0.0036 0.0201 0.1809 0.8565 -0.0357 0.0430
constant 1.1692 0.0231 50.566 0.0000 1.1239 1.2145
==============================================================================
Что соответствует выходу Stata.На самом деле не было никакой причины отказаться от 1987 года, мы могли бы выбрать что-то другое, но по крайней мере с этим мы можем увидеть результаты, соответствующие xtreg.