Pooled model SUR

I am trying to estimate a pooled SUR in R.
I am not able to replicate the results given in text.

Link for dataset https://drive.google.com/file/d/0B11woa7YyVb1VjJUcXlZeHI0M1E/view?usp=sharing

I am also attaching my code and the example in book.

![image2|350x500]

rm(list=ls())

costdata= read.csv("TableF10-2.csv")
install.packages("micEcon")
library(micEcon)
library(tidyverse)
costdata= costdata %>%
  mutate(lnpkpm= log(costdata$Pk)/log(costdata$Pm)) %>%
  mutate(lnplpm= log(costdata$Pl)/log(costdata$Pm)) %>%
  mutate(lnpepm= log(costdata$Pe)/log(costdata$Pm)) %>%
  mutate(sk= (costdata$K*costdata$Pk)/costdata$Cost) %>%
  mutate(sl= (costdata$L*costdata$Pl)/costdata$Cost) %>%
  mutate(se= (costdata$E*costdata$Pe)/costdata$Cost)

View(costdata)
library(systemfit)

eqn_cpaital= sk~ lnpkpm+ lnplpm+ lnpepm
eqn_labour= sl~ lnpkpm+ lnplpm+ lnpepm
eqn_energy= se~ lnpkpm+ lnplpm+ lnpepm
system1=list(capital= eqn_cpaital, labour= eqn_labour, energy= eqn_energy)

#SUR
costdata_SUR= systemfit(system1, method = "SUR", data = costdata,  methodResidCov= "noDfCor", pooled = TRUE, residCovWeighted = TRUE )
summary(costdata_SUR)

Any help is appreciated.
Thanks!

image4

What output are you getting when you run summary(costdata_SUR)?
Is the script running all the way through?
It's a bit hard to read the text from the images β€” but there is a nice little SUR FAQ section from UCLA's IDRE
https://stats.idre.ucla.edu/r/faq/how-can-i-perform-seemingly-unrelated-regression-in-r/

1 Like

It's running through.

systemfit results 
method: SUR 

        N DF SSR detRCov   OLS-R2 McElroy-R2
system 72 60   0       0 0.803285   0.893478

         N DF SSR MSE    RMSE       R2   Adj R2
capital 24 20   0   0 3.0e-05 0.629199 0.573579
labour  24 20   0   0 6.7e-05 0.794454 0.763623
energy  24 20   0   0 2.1e-05 0.920042 0.908048

The covariance matrix of the residuals used for estimation
            capital      labour      energy
capital 7.69646e-10 1.19362e-09 4.26354e-10
labour  1.19362e-09 3.78213e-09 9.95445e-10
energy  4.26354e-10 9.95445e-10 3.58730e-10

The covariance matrix of the residuals
            capital      labour      energy
capital 7.69646e-10 1.19362e-09 4.26354e-10
labour  1.19362e-09 3.78213e-09 9.95445e-10
energy  4.26354e-10 9.95445e-10 3.58730e-10

The correlations of the residuals
         capital   labour   energy
capital 1.000000 0.699601 0.811410
labour  0.699601 1.000000 0.854604
energy  0.811410 0.854604 1.000000


SUR estimates for 'capital' (equation 1)
Model Formula: sk ~ lnpkpm + lnplpm + lnpepm

                Estimate   Std. Error  t value   Pr(>|t|)    
(Intercept)  2.74321e-04  4.93844e-05  5.55481 1.9469e-05 ***
lnpkpm       1.31729e-05  5.60644e-06  2.34960  0.0291753 *  
lnplpm      -8.49374e-05  2.60368e-05 -3.26220  0.0039001 ** 
lnpepm       5.62066e-05  8.91539e-06  6.30445 3.7346e-06 ***
---
Signif. codes:  0 β€˜***’ 0.001 β€˜**’ 0.01 β€˜*’ 0.05 β€˜.’ 0.1 β€˜ ’ 1

Residual standard error: 3e-05 on 20 degrees of freedom
Number of observations: 24 Degrees of Freedom: 20 
SSR: 0 MSE: 0 Root MSE: 3e-05 
Multiple R-Squared: 0.629199 Adjusted R-Squared: 0.573579 


SUR estimates for 'labour' (equation 2)
Model Formula: sl ~ lnpkpm + lnplpm + lnpepm

                Estimate   Std. Error  t value   Pr(>|t|)    
(Intercept)  1.46879e-03  1.09474e-04 13.41675 1.8472e-11 ***
lnpkpm      -1.76772e-05  1.24282e-05 -1.42234   0.170337    
lnplpm      -1.59515e-04  5.77179e-05 -2.76370   0.011979 *  
lnpepm       1.45536e-04  1.97634e-05  7.36390 4.0910e-07 ***
---
Signif. codes:  0 β€˜***’ 0.001 β€˜**’ 0.01 β€˜*’ 0.05 β€˜.’ 0.1 β€˜ ’ 1

Residual standard error: 6.7e-05 on 20 degrees of freedom
Number of observations: 24 Degrees of Freedom: 20 
SSR: 0 MSE: 0 Root MSE: 6.7e-05 
Multiple R-Squared: 0.794454 Adjusted R-Squared: 0.763623 


SUR estimates for 'energy' (equation 3)
Model Formula: se ~ lnpkpm + lnplpm + lnpepm

                Estimate   Std. Error  t value   Pr(>|t|)    
(Intercept)  3.52891e-04  3.37154e-05 10.46675 1.4612e-09 ***
lnpkpm      -1.37385e-05  3.82759e-06 -3.58932  0.0018333 ** 
lnplpm      -1.29858e-04  1.77757e-05 -7.30536 4.6046e-07 ***
lnpepm       7.39568e-05  6.08666e-06 12.15065 1.0902e-10 ***
---
Signif. codes:  0 β€˜***’ 0.001 β€˜**’ 0.01 β€˜*’ 0.05 β€˜.’ 0.1 β€˜ ’ 1

Residual standard error: 2.1e-05 on 20 degrees of freedom
Number of observations: 24 Degrees of Freedom: 20 
SSR: 0 MSE: 0 Root MSE: 2.1e-05 
Multiple R-Squared: 0.920042 Adjusted R-Squared: 0.908048 

So it's the actual values that are coming out wrong? Sorryβ€” I honestly don't know much about SUR, so I couldn't tell if you were having a code issue or a calculation one!

Yes the actual values are wrong.

It's been a long time since I read Greene, and only ever ran one SUR (in Stata, not R), so I could be completely off.

First, I think the log quantities in your code should be lnpkpm = log(Pk / Pm) instead of lnpkpm = log(Pk) / log(Pm) given the equations. Next, I think K, L, and E are already factor shares, according to the descripton of those variables on the dataset page. Finally, I think the SUR estimates are in Table 10.3; Table 10.4 involves using those estimates and fitted cost shares to get the relevant elasticities.

That said, I tried my corrections and got estimates similar in magnitude and direction to those in Table 10.3, but not exactly those numbers. I'm not sure why...it might mean I'm totally off, so I would take what I say here with a pinch of salt, and defer to those with more experience with the translog method.

2 Likes

Apologies for not posting table 10.3.

The estimates are still not matching, even after emploing the suggestions.

Any further help is appreciated.
Thanks.

Nobody has any suggestion??

You're outside my domain of expertise!

I'm not sure if the economics Stack Exchange site is "mature" yet, but you might give this question a go either there, on Cross-Validated:

There's also a Principles of Econometrics in R bookdown by Constantin Colonescu, that has an SUR model in this example.

If you can break it down to a specific piece of code that isn't evaluating as expected, you might have a better chance getting help here!

Hopefully these suggestions can meta-help you find the right help!

3 Likes

Thanks,
Happy Birthday!!

1 Like

Howdy. Other than the setup, this isn't really an R or RStudio question. But happy to help. It's been 8 years for me spending much time with Greene

The model setup by eq 10-40 is, I think:

costdata = read.csv(".../TableF10-2.csv")

# setup dependent vars for /delta_...
costdata = costdata %>%
  mutate(
    lnpkpm = log(Pk / Pm),
    lnplpm = log(Pl / Pm),
    lnpepm = log(Pe / Pm)
    )

setup eq 10-40:

df_pooledreg <-
  bind_rows(
    tibble(
      y = costdata$K,
      k = 1,
      l = 0,
      e = 0,
      kk = costdata$lnpkpm,
      kl = costdata$lnplpm,
      ke = costdata$lnpepm,
      ll = 0,
      le = 0,
      ee = 0,
      year = costdata$Year
    ),
    tibble(
      y = costdata$L,
      k = 0,
      l = 1,
      e = 0,
      kk = 0,
      kl = costdata$lnpkpm,
      ke = 0, 
      ll = costdata$lnplpm,
      le = costdata$lnpepm,
      ee = 0,
      year = costdata$Year
    ),
    tibble(
      y = costdata$E,
      k = 0,
      l = 0,
      e = 1,
      kk = 0,
      kl = 0,
      ke = costdata$lnpkpm, 
      ll = 0,
      le = costdata$lnplpm,
      ee = costdata$lnpepm,
      year = costdata$Year
    )
  )

run "The Pooled Model" eq 10-19 in Greene chapter 10.
Note, eq 10-19 doesn't have an intercept.

pooledreg = glm(
  y ~ 0 + k + l + e + kk + kl + ke + ll +le + ee,
  data = df_pooledreg
)
summary(pooledreg)

I get the following coefficients:

Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
k   0.056223   0.001888  29.773  < 2e-16 ***
l   0.253409   0.001852 136.834  < 2e-16 ***
e   0.041861   0.002291  18.272  < 2e-16 ***
kk  0.030426   0.008067   3.772 0.000349 ***
kl  0.001755   0.004556   0.385 0.701342    
ke -0.003562   0.007481  -0.476 0.635562    
ll  0.075127   0.005280  14.228  < 2e-16 ***
le  0.003252   0.005756   0.565 0.574043    
ee  0.046777   0.017415   2.686 0.009136 ** 
---

Now, Table 10.4. The first line, Fitted shares is just:

fitted_k = predict(
  pooledreg,
  df_pooledreg %>% filter(year == 1959 & k)
)
fitted_l = predict(
  pooledreg,
  df_pooledreg %>% filter(year == 1959 & l)
)
fitted_e = predict(
  pooledreg,
  df_pooledreg %>% filter(year == 1959 & e)
)
# Note that assuming a constant returns to scale production function, 
# `m` is just `1-(k+l+e) `
fitted_m = 1 - (fitted_k + fitted_l + fitted_e) 

fitted_k: 0.05655583
fitted_l: 0.2747627
fitted_e: 0.04487152
fitted_m: 0.6238099


How do we estimate the derivatives involving Materials? i.e. those $/delta_{M...}$? For that, turn to equation 10-38 (well, it's eq 10-38 in my Greene, 6th edition of econometric analysis anyway. they seem to be moving things around). That equation just says,
$\Sigma_{i=1}^M \delta_{ij} = 0$ and $\Sigma_{j=1}^M \delta_{ij} = 0$.

Really simple, you can see this in table 10.4. It just says that $\delta_{KK} + $\delta_{KL} + $\delta_{KE} + $\delta_{KM} = 0.
You'll notice in table 10.4 any derivative with an M has an asterisks next to it, referring to eq 10-38.


Equation 10-39 gives the own and cross elasticities of substitution.
So with $\theta_{kk}$ for example, (/delta_{kk} + s_{i}(s_{i} - 1)) / s_{i}^2, where i is capital share's output in that year.
or with my estimates and your data:

s_k1959 = ((costdata %>% filter(Year == 1959)))$K
delta_kk = pooledreg$coefficients['kk'] %>% as.vector
(delta_kk + s_k1959 * (s_k1959 - 1)) / (s_k1959^2)
[1] -7.214395 # that's under "Implied Elasticities of Substitution, 1959" for Capital/Capital

And for cross price elasticity of substitution, eq 10-39 says
$\theta_{ij}$ is , (/delta_{ij} + s_{i} * s{j}) / (s_{i} * s_{j}),
for example for $\theta_{KL}$:

s_k1959 = ((costdata %>% filter(Year == 1959)))$K
s_l1959 = ((costdata %>% filter(Year == 1959)))$L
delta_kl = pooledreg$coefficients['kl'] %>% as.vector
(delta_kl + (s_k1959 * s_l1959)) / (s_k1959 * s_l1959)
[1] 1.103912 # this is our estimate for "Implied Elasticities of Substitution, 1959" for Labor/Capital 

There are slight deviations. They might come from using an updated dataset...?

3 Likes

Thanks Curtis the values are indeed slightly different with my working also. I just want to know why the pooled SUR from systemfit fails here??

Does the systemfit package throw an error/the code itself fail somehow? If so, you might consider reporting it as a bug or feature request on systemfit's R-Forge page. Here's the direct link to its trackers.

No the output produced is different and the code runs through. I am already intouch with the developers there, Thanks.