Simulate nested logit errors in Python

Question

I am trying to simulate 1,000 draws of the error term 𝜖_𝑗𝑡𝑖 in a two-level nested logit model in Python. For simplicity, I assume that there are 250 individuals, and the choice set consists of 4 products A1, A2, B1, B2 belonging to nests A and B respectively.

I am looking at the specification from the PyBLP documentation, which decomposes the error term into two terms (see link - sorry, I am new and don't have enough credits to post images).

I can draw the second term using, e.g., samples = np.random.gumbel(0, 1, 1000). But I have no idea how to draw the first term such that it does not vary within an individual-nest, yet the two terms (with the second term weighted by 1-rho) still sum up to be Type 1 Extreme Value. I also don't see how this sum can be iid even though the documentation says so, given the existence of the first term.

Perhaps this points to a misunderstanding of how the error term is specified. Grateful if anyone can give tips or correct me.

I am aware that the question has been answered here for R, but do not understand what exactly happened or how that translates to Python.

Thanks a lot!