I want to estimate a system of 4 equations, where the Y-Variables sum to one approximately, and so the error terms are correlated. The data set contains ~180 variables and ~1600 observations. I tried running sureg with i.variable in the specification to add in fixed effects. So there are three issues that I am facing:

1. The co-efficients with systemfit package in R and Stata are different. (They were the same when I ran felm() in R and reghdfe in Stata.) So, I am not able to understand which one is consistent.
2. I want to add one more factor variable, which is correlate with another X-variable.
3. I want to add clustered standard errors to sureg.

Could you please help me with these issues, or point me to some research papers which used SUR with fixed effects?

I am relatively new to data, so any help on this front would be very valuable to me.

I am happy to discuss the data problem in more detail in personal chat.

Thank You