Add support for two-way clustering about pyfixest HOT 4 CLOSED

from pyfixest.estimation import fepois
from pyfixest.summarize import summary
from pyfixest.utils import get_data

data = get_data(N = 1000, model = "Fepois")
fit_crv1 = fepois("Y ~ X1 + X2 | f1 + f2", data = data, vcov = {"CRV1":"f1+f2"})
fit_crv3 = fepois("Y ~ X1 + X2 | f1 + f2", data = data, vcov = {"CRV3":"f1+f2"})

summary([fit_crv1, fit_crv3])

###

Estimation:  Poisson
Dep. var.: Y, Fixed effects: f1+f2
Inference:  CRV1
Observations:  997

| Coefficient   |   Estimate |   Std. Error |   t value |   Pr(>|t|) |   2.5 % |   97.5 % |
|:--------------|-----------:|-------------:|----------:|-----------:|--------:|---------:|
| X1            |     -0.052 |        0.035 |    -1.502 |      0.133 |  -0.120 |    0.016 |
| X2            |      0.007 |        0.014 |     0.520 |      0.603 |  -0.020 |    0.034 |
---
Deviance: 1066.019
###

Estimation:  Poisson
Dep. var.: Y, Fixed effects: f1+f2
Inference:  CRV3
Observations:  997

| Coefficient   |   Estimate |   Std. Error |   t value |   Pr(>|t|) |   2.5 % |   97.5 % |
|:--------------|-----------:|-------------:|----------:|-----------:|--------:|---------:|
| X1            |     -0.052 |        0.038 |    -1.375 |      0.169 |  -0.126 |    0.022 |
| X2            |      0.007 |        0.015 |     0.480 |      0.631 |  -0.022 |    0.036 |
---
Deviance: 1066.019

fit = fepois("Y ~ X1 + X2", vcov= {'CRV1': "f1+f2"}, data=df)