GLM multidimensional Python

Question

I have four variables with 50 observations each and I am trying to fit a GLM using statsmodels. There are two endogenous variables, say y1 and y2, and two exogenous variables, say x1 and x2. My first idea was

 Create example data with 5 observations
y1 = np.array([2, 3, 4, 5, 6])
y2 = np.array([1, 2, 2, 3, 4])
x1 = np.array([0.1, 0.2, 0.3, 0.4, 0.5])
x2 = np.array([0.5, 0.4, 0.3, 0.2, 0.1])

# Fit a Poisson GLM for y1
poisson_model_y1 = sm.GLM(y1, sm.add_constant(x1), family=sm.families.Poisson())
poisson_results_y1 = poisson_model_y1.fit()
print("Poisson GLM for y1:")
print(poisson_results_y1.summary())

# Fit a Poisson GLM for y2
poisson_model_y2 = sm.GLM(y2, sm.add_constant(x2), family=sm.families.Poisson())
poisson_results_y2 = poisson_model_y2.fit()
print("\nPoisson GLM for y2:")
print(poisson_results_y2.summary())

Is it possible jointly estimate; that is with a vector of [y1 y2] and a vector [x1 x2].