Total factor productivity (TFP) estimation in R

Question

I'm trying to calculate TFP with the estprod package in R because this package allows me to calculate with Gross-Output (like with levinsohn petrin) but I have trouble getting omega

While the prodest package can calculate omega but I can't add Gross-Output.

Has anyone calculated the TFP in R?

[updated] Reproducible example

library(prodest)
library(estprod)

data = structure(list(Y = c(12.8268952320028, 12.8521337216858, 12.8297798993317, 12.6357098405348, 
                            12.4004643557221, 12.3600068590732, 12.826153741752, 12.8990697648315, 
                            12.8393196511736, 12.7619091581355), 
                      Labor = c(18.2268433780233, 18.1280743937883, 
                            16.2004953940248, 18.4010948157752, 16.1374389570964, 18.3902253199657, 
                            16.2359220634517, 16.2468898407345, 16.2720069307031, 16.2673373011924), 
                      Capital = c(12.1642906069818, 12.088450259469, 
                           12.0237634086603, 12.0750443461238, 12.0098651639831, 12.0601371534517, 
                           11.9679864443247, 11.9299321784366, 12.0016971178669, 12.1209493577051), 
                      Materials = c(18.1472736836291, 17.7735308597061, 
                          17.1908332586684, 18.3914858811768, 18.2747180360049, 16.9304018090212, 
                          18.4414110290762, 18.0188590654862, 18.1427756804037, 15.805209131618), 
                      id = c(1070101, 1070101, 1070101, 1070101, 
                              1070101, 1070101, 1070102, 1070102, 1070102, 1070102), 
                      year = c(2013,2014, 2015, 2016, 2017, 2018, 2013, 2014, 2015, 2016)), 
                 class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, -10L))

# Y, Labor, Capital, Materials are in logarithms
mod1 = estprod::levinsohn_petrin(data,
                                 formula = Y ~ Labor | Capital | Materials, 
                                 id = "id", 
                                 time = "year", 
                                 reps = 20,
                                 gross = T) #estprod with Gross-Output = T

#how can i calculate omega with estprod?

mod2 = prodest::prodestLP(data$Y, # prodest does not have the parameter Gross-Output
             fX = data$Labor,
                 sX = data$Capital, 
             pX = data$Materials,
                 idvar = data$id,
             timevar = data$year,
                 opt='optim',
                 exit = F,
                 tol = 1e-100)  

omega = prodest::omega(mod2) #with 'prodest' I can calculate omega

Answer 1

Following the prodtest manual they define omega as residuals from the model, that is

$$\omega_{it} = y_{it} − (\alpha + w_{it} \beta + k_{it}\gamma)=y_{it}-\hat{y}_{it}$$

estprod does not have any comprehensible documentation but one way how we can get the above would be to extract the coefficients from the model and then calculate $y-\hat{y}$ manually.

We can do that with the broom package:

install.packages("broom")
require(broom) 
coefs<-tidy(mod1)
# A tibble: 3 x 3
  statistic   bias std.error
      <dbl>  <dbl>     <dbl>
1     -3.71   4.77      2.35
2     10.5  -11.8       4.26
3      1.73  -2.14      1.64

the names are off but the first column gives the same coefficients as when I inspect the model with summary so they are correct.

Next we can just can just calculate $y-\hat{y}$ as

omega_1 = (data$Y - data$Labor*coefs$statistic[1] - data$Capital*coefs$statistic[2] - data$Materials*coefs$statistic[3])  

the code is inelegant but it should do the job.

As a side note it might be worth while considering to estimate the TFP without these packages. Normally packages are the R's biggest strength but in this case estprod does not have any good accompanying documentation or manual and the prodest does not have the possibility to add the gross output yet as in the link the devs mentioned. In cases like this it might be better to just program your own function.