I need to estimate a Seemingly Unrelated Regression (SUR) model that accounts for the temporal dimension of the data. Specifically, I am having trouble with the equation formula and incorporating the spatial dimension of the data. For a bit more context, the SUR model allows for different coefficients for each time period and accounts for inter-temporal dependence through the covariance matrix of the system of regression equations.

I have data for 5 years (in the “year” column), 32 regions (in the “CVE_ENT” column), and 4 effects (“std_region_effect”, “std_industrymix_effect”, “std_interaction_effect”, “std_allocation_effect”). The dependent variable is state productivity (in the “estate_labprod” column).

I am having trouble specifying and estimating such a model, despite having read the documentation for the spsur package and consulting various GPTs for hours. Any help on this would be greatly appreciated.

Here is a sample of my data (gdf) and the spatial weights matrix (2 years and 2 regions):

gdf_sample2 <- st_sf(
  year = c(1998, 1998, 2003, 2003, 1998, 2003),
  CVE_ENT = c(1, 1, 1, 1, 2, 2),
  estate_labprod = c(136, 136, 234, 234, 200, 210),
  std_region_effect = c(234.9, 234.9, 125.1, 125.1, 105.7, 110.2),
  std_industrymix_effect = c(8.97, 8.97, 4.01, 4.01, 10.2, 9.5),
  std_interaction_effect = c(305.3, 305.3, 141.8, 141.8, 508.7, 495.6),
  std_allocation_effect = c(-336, -336, -272, -272, -169, -175),
  geometry = st_sfc(
    st_multipolygon(list(st_polygon(list(rbind(c(-102, 22), c(-101, 22), c(-101, 23), c(-102, 23), c(-102, 22)))))),
    st_multipolygon(list(st_polygon(list(rbind(c(-102, 22), c(-101, 22), c(-101, 23), c(-102, 23), c(-102, 22)))))),
    st_multipolygon(list(st_polygon(list(rbind(c(-102, 22), c(-101, 22), c(-101, 23), c(-102, 23), c(-102, 22)))))),
    st_multipolygon(list(st_polygon(list(rbind(c(-102, 22), c(-101, 22), c(-101, 23), c(-102, 23), c(-102, 22)))))),
    st_multipolygon(list(st_polygon(list(rbind(c(-103, 21), c(-102, 21), c(-102, 22), c(-103, 22), c(-103, 21)))))),
    st_multipolygon(list(st_polygon(list(rbind(c(-103, 21), c(-102, 21), c(-102, 22), c(-103, 22), c(-103, 21))))))
  )
)

# Queen contiguity matrix (standarized)
neighbors_sample <- poly2nb(gdf_sample2, queen = TRUE)
neighbors_stand2 <- nb2listw(neighbors2, style = "W", zero.policy = TRUE)

# this is a kind of formula that i want but won't work
Formula(estate_labprod ~ std_region_effect |
          estate_labprod ~ std_industrymix_effect |
          estate_labprod ~ std_interaction_effect |
          estate_labprod ~ std_allocation_effect)

spsurml_result <- spsurml(formula = formula, #do not know how to do it
                          data = gdf,
                          listw = neighbors_stand2, #Queen contiguity standarized matrix
                          type = "sim", #without spatial effects
                          method = "eigen",
                          G = 4,  # Equations
                          N = 2,  # spatial unites (sample)
                          Tm = 2  # time periods (sample)
)

In summary, my question is based on:

• I don’t know how to incorporate the temporal dimension.
• I don’t know how to specify the formula.

The specification of the functional form would be something like this: