Replicate Blundell and Bond (2000) results using R

Question

I want to replicate Blundell and Bond (2000) Table III in R. I'm using the function pgmm from package plm, which (apparently) replicates the approach of Stata's xtabond2 (Roonman, 2009). The estimates I get are qualitatively similar, but not the same, the Arellano-Bond test of serial correlation suggests the moment conditions are invalide....and I'm sad to say I don't yet see how to recover the production function coefficients.

The authors of plm have a useful companion textbook: Millo and Croissant 2018 Panel Data Econometrics in R. I get the data from there (associated packaged pder), and follow section 7.5 for set-up and diagnostics. Here is a self-contained example:

library(plm)
library(pder) # contains the Blundell and Bond dataset
#detach(package:dplyr) # if loaded, dplyr::lag masks plm::lag

data('RDPerfComp')
str(RDPerfComp) #n.b. - variables already in logs (see ?RDPerfComp)

And now to replicate some of the columns from Blundell and Bond 2000 Table 3:

# "OLS levels"
ols<-plm(y~lag(y,1) + lag(n,0:1) + lag(k,0:1),
             data=RDPerfComp, index=c("id","year"), 
             effect="twoways", model = "pooling")
summary(ols)

# "Within groups"
within_groups<-plm(y~lag(y,1) + lag(n,0:1) + lag(k,0:1),
                   data=RDPerfComp, index=c("id","year"), 
                   effect="twoways", model = "within")
summary(within_groups)

# "SYS t-2"
sys_t2<-pgmm(
  y ~ lag(y,1) + lag(n,0:1) + lag(k,0:1) |
    lag(y,2:99) + lag(n,2:99) + lag(k,2:99),
  data=RDPerfComp, index=c("id","year"),
  effect="twoways",model="onestep",
  transformation="ld"
)
summary(sys_t2,robust=TRUE,time.dummies=FALSE) # both options can be provided directly in pgmm()

# "SYS t-3"
sys_t3<-pgmm(
  y ~ lag(y,1) + lag(n,0:1) + lag(k,0:1) |
    lag(y,3:99) + lag(n,3:99) + lag(k,3:99),
  data=RDPerfComp, index=c("id","year"),
  effect="twoways",model="onestep",
  transformation="ld"
)
summary(sys_t3,robust=TRUE)

The results are very similar for the pooled OLS model, almost identical for fixed effects/within model, and only "qualitatively similar" (magnitudes, signs) for the SYS-GMM models, with some exceptions (notable the estimate for the y_t-1 coefficient). More concerning, the Sargan-Hansen test appears to reject the null of valid instruments for both the SYS-GMM models:

sargan(sys_t2) # same statistic as reported by summary()
sargan(sys_t3)

...and the Arellano-Bond moment condition test appears to reject the null of no serial correlation, suggesting that the moment conditions are invalid

mtest(sys_t2) # you can make this robust using vcovc =
mtest(sys_t3)

So my questions:

1. Am I specifying these wrong (probably - my first try with gmm) or are we just looking at the results of different machinery?
2. What could explain the diagnostic test results, which appear to say the SYS-GMM models aren't valid?
3. How do we get Beta_n, Beta_k, and rho (the production function parameters, as given in the middle panel of Table III) from these results? (EDIT: OK, as stated in the paper, via testing and imposing the coefficient restrictions using "a minimum distance estimator" see eq'n 2.3 and text below in original paper)

(yes these are 'advanced undergrad' type questions, no this is not coursework)

Answer 1

3. How do we recover parameters from production function estimates (INCOMPLETE ANSWER - will be updated with how to do this in R once I have time to figure it out, or if somebody else knows...)

Blundell and Bond aren't estimating the parameters of a Cobb-Douglas production function. They're estimates the parameters of a "dynamic (common factor) representation" of a production function, specifically:

$(1) \ \ y_{it} = \beta_n n_{i,t} - \rho \beta_n n_{i,t-1} + \beta_k k_{i,t} - \rho \beta_k k_{i,t} + \rho y_{i,t-1} + (\gamma_t - \rho \gamma_{t-1}) + (\eta_i(1-\rho) + e_{it} + m_{it} - \rho m_{i,t-1})$

...which they get by quasi-differencing their production function and doing something weird with the $v_{it}$ term (not shown above). Changing the symbols, they re-write this as:

$(2) \ \ y_{it} = \pi_1n_{it} + \pi_2 n_{i,t-1} + \pi_3 k_{it} + \pi_4 k_{i,t-1} + \pi_5 y_{i,t-1} + \gamma^*_t + (\eta^*_i +w_{it})$

Where the * presumably indicates steady-state values, but again, I wish I had access to a teacher to explain precisely why they write that. Anyway, what their SYS-GMM estimator returns is a vector of parameter values $\boldsymbol{\pi} = (\pi_1,...,\pi_5)'$. We now have to use this to get the vector of parameters $\boldsymbol{\theta} = (\beta_n,\beta_k,\rho)$.

The only place I can find an explanation of how one does that is this paper applying SYS-GMM to Cuba. They refer to Woolridge 2002 (I don't have it), and assert that one can consider the problem as defined by the relationship $\boldsymbol{\pi} = h(\boldsymbol{\theta})$. We have $\boldsymbol{\pi}$, we need $\boldsymbol{\theta}$, and we will have to employ some function h() to make the link. It sounds like h() is taken to be an instance of a general class of minimum distance estimators, and we can (according to Woolridge 2002, apparently) solve the problem by solving:

$min (\boldsymbol{\theta}) = (\hat{\pi}-h(\theta))'\Omega(\hat{\pi}-h(\theta)) $

where $\Omega$ is some weighting function. And again, apparently $\Omega$ should be chosen as the variance of the matrix of parameters, obtained by "the delta method".

...but now it gets confusing, because although it's easy to find info on the delta method online, the matrix of parameters includes the unknowns $\boldsymbol{\theta}$. I guess the way to tackle this is to multiply it all out and plug into optim, using some arbitrary choice on starting values (factor shares of output, if $ valued?). I will update if I get a chance to try that...

In going from equation (1) to (2) it is pretty obvious that:

$$h(\beta_n,\beta_k,\rho) = \begin{pmatrix} \beta_n \\ - \rho \beta_n \\ \beta_k \\ - \rho \beta_k \\ \rho \end{pmatrix}= \begin{pmatrix} \pi_1 \\ \pi_2 \\ \pi_3 \\ \pi_4 \\ \pi_5 \end{pmatrix},$$

and $$\Omega = var(\pi)$$