Your problem seem almost the same as the issue you linked, just that you will do the indexing on the betas, not the intercept. So you will just make one vector of betas for all the times, and one vector of betas for the locations, then use fancy indexing to expand them to df.shape[0] add them together, like y = (beta_location[loc_idx] + beta_year[year_idx]) * X. Have you tried something like that?