IV regression in Python with clustered robust SE (maybe)

ivreg_in_python.py

# Native Python packages
import numpy as np
import pandas as pd

# R2pi package
from rpy2.robjects import numpy2ri, pandas2ri
from rpy2.robjects.packages import importr

# R imports (Assuming AER and ivpack are installed in R)
base = importr("base")
aer = importr("AER")
ivpack = importr("ivpack")

# Automatic conversion btw R Vectors and Python objects
numpy2ri.activate()
pandas2ri.activate()

# Some fake data (IV not really needed here)
n = 10000
df = pd.DataFrame(data=np.random.multivariate_normal(
                            mean=[0] * 5,
                            cov=np.ones((5,5)) + np.eye(5),
                            size=n),
                  columns =["Y", "X1", "X2", "Z1", "Z2"])
clusterid = np.tile([1,2,3,4,5], (n//5, ))


# Example regression
model = aer.ivreg(formula="Y ~ X1 + X2 | Z1 + Z2", data=df)

# Example clustering (Needs review! I've never used this package before.)
robust_se = ivpack.cluster_robust_se(model, clusterid)

# Tidying up
robust_se = pd.DataFrame(
        data=np.array(robust_se),
        columns=['Estimate', 'Std. Error', 't value', 'Pr(>|t|)'],
        index=['(Intercept)', 'X1', 'X2']
)
print(robust_se)