##### Replication of the naïve model
In what follows, we replicate the naïve model in [[Statistical modeling of IV]] using Python's recent *[[PyFixest]]* library.
```python
import pandas as pd
import pyfixest as pf
df = pd.read_csv('d:/causalbook/data/processed-tobacco-tax-as-iv_v101.csv')
pf.feols("pack_sales_per_capita ~ price_per_pack | state + year", data = df).summary()
```
The results are as follows:
```python
Estimation: OLS
Dep. var.: pack_sales_per_capita, Fixed effects: state+year
Inference: CRV1
Observations: 2550
| Coefficient | Estimate | Std. Error | t value | Pr(>|t|) | 2.5% | 97.5% |
|:---------------|-----------:|-------------:|----------:|-----------:|--------:|--------:|
| price_per_pack | -6.625 | 2.002 | -3.310 | 0.002 | -10.646 | -2.605 |
---
RMSE: 13.537 R2: 0.894 R2 Within: 0.063
```
The only difference from the results in R is the standard error. This is due to the difference in the small sample correction as discussed in [[Statistical modeling of IV]]. Previously, we shared the parameters to add to make the R results the same as Stata and `R/plm`. To make the R results exactly the same as the results from `PyFixest`, only one of the small sample corrections (cluster adjustment) needs to be disabled: `ssc = fixest::ssc(cluster.adj = FALSE)` Add this to the `fixest` function in R and you will see that the standard error to be reported is `2.002`.
##### Replication of the IV model
The following is the replication of the IV model in [[Statistical modeling of IV]] using Python's `PyFixest` library:
```python
pf.feols("pack_sales_per_capita ~ 1 | state + year | price_per_pack ~ state_tax_per_pack", data = df).summary()
```
The results are as follows:
```python
Estimation: IV
Dep. var.: pack_sales_per_capita, Fixed effects: state+year
Inference: CRV1
Observations: 2550
| Coefficient | Estimate | Std. Error | t value | Pr(>|t|) | 2.5% | 97.5% |
|:---------------|-----------:|-------------:|----------:|-----------:|--------:|--------:|
| price_per_pack | -7.815 | 1.999 | -3.910 | 0.000 | -11.829 | -3.800 |
```
Note that the results are still identical. To get the same standard error in R, the same parameter-value pair must be added to the R function: `ssc = fixest::ssc(cluster.adj = FALSE)`.
##### Replication using Python's linearmodels
Before `R/fixest` was ported to Python as `PyFixest`, we replicated both the naïve model and the IV model in Python using `linearmodels`. For a panel model with fixed effects, `linearmodels` provides a cumbersome implementation with some hiccups in the standard errors.[^1] `fixest` and `PyFixest` report the correct standard errors adjusted for within-group (state) correlation. As discussed in [[Clustered standard errors]], while we know exactly how the standard errors are computed in `R/fixest` and `Python/PyFixest` and we know it is the correct way, the `Python/linearmodels` documentation is more mysterious. After a deep dive into the behavior of `linearmodels`, we realized that we can get similar standard errors from it if we also include time as a clustering variable. But is that what we want to do here?
> [!note]
> Originally, we wanted to cluster by state. Could we have also clustered by time? Yes, we could. The outcome variable is pack sales, so there may be shocks that are arbitrarily correlated across (not *within*) states, but in a given year. If the outcome in all states is hit by the shocks in the same way, then the time fixed effects we included in the model should solve the problem. But what if there is a dependency between some states (say, neighboring states) that makes them differently vulnerable to the shocks than the other states? Consider New England states, for example. Could the pack sales w.r.t. price in the New England states be affected differently than in the rest of the states in a given year due to some shocks that are not accounted for in the model? This is a far-fetched scenario in most cases, but yes, it may be plausible to apply time clustering.
The real problem here is that the clustering of the standard errors in `Python/linearmodels` is not consistent with `R/fixest`, nor with `Python/PyFixest`, even when time clustering is added to the model. The mystery behind the standard errors is due to a different reason. See [[Oh my! Different standard errors everywhere]] to learn the reason behind such inconsistencies and how to get identical standard errors across languages and packages.
[^1]: Since then, we have also changed our data, but if you want to see the linearmodels code, you can still see it here: [[Replication using linearmodels]].
> [!info]- Last updated: January 18, 2025