2
$\begingroup$

I would like to use a two stage least squares approach (2SLS), where the first stage would benefit from a Tobit specification.

I cross posted this on stackoverflow because there might be quite some coding involved. I have some sample data as follows:

panelID= c(1:50)
year= c(2005, 2010)
country = c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J")
urban = c("A", "B", "C")
indust = c("D", "E", "F")
sizes = c(1,2,3,4,5)
n <- 2
library(AER)
library(data.table)
library(dplyr)
set.seed(123)
DT <- data.table(   country = rep(sample(country, length(panelID), replace = T), each = n),
                    year = c(replicate(length(panelID), sample(year, n))),
                    sales= round(rnorm(10,10,10),2),
                    industry = rep(sample(indust, length(panelID), replace = T), each = n),
                    urbanisation = rep(sample(urban, length(panelID), replace = T), each = n),
                    size = rep(sample(sizes, length(panelID), replace = T), each = n))
DT <- DT %>%
group_by(country) %>%
mutate(base_rate = as.integer(runif(1, 12.5, 37.5))) %>%
group_by(country, year) %>%
mutate(taxrate = base_rate + as.integer(runif(1,-2.5,+2.5)))
DT <- DT %>%
group_by(country, year) %>%
mutate(vote = sample(c(0,1),1), 
votewon = ifelse(vote==1, sample(c(0,1),1),0))

Say I would like to run:

summary(ivreg(sales ~ taxrate + as.factor(industry) + as.factor(size) + as.factor(urbanisation) + as.factor(vote) | as.factor(votewon) + as.factor(industry) + as.factor(size) + as.factor(urbanisation) + as.factor(vote), data=DT))

But I want the first stage to be a tobit specification? Something like:

tobit_reg <<- censReg(taxrate ~ as.factor(votewon) + as.factor(industry) + as.factor(size) + as.factor(urbanisation) + as.factor(vote), left=3, right=15, data=DT)
summary(tobit_reg)

How can I have the tobit regression as the first stage for the 2SLS?

I would additionally like to test for over identification.

$\endgroup$

1 Answer 1

4
$\begingroup$

You could perform a 2SLS approach by hand, where you estimate a Tobit model of the instrument in the first stage and use the fitted values $\hat y_2$ to estimate an OLS model in the second stage.

The naïve standard errors won't be valid though, because they do not consider that the $\hat y_2$ itself is an estimate; only the variance of the residuals of the second stage are used to calculate the standard errors.

Therefore we have to correct the variance-covariance matrix (VCOV) by applying the correct RMSE using the formula provided in Cameron, A., & Trivedi, P. (2005:43).

$$\hat V[\hat \beta_{2SLS}] = N[X'P_ZX]^{-1}[X'Z(Z'Z)^{-1}\hat S(Z'Z)^{-1}Z'X][X'P_ZX]^{-1}$$

where

$$\hat S=N^{-1}\sum _i \hat u_i^2z_iz'_i$$

and

$$\hat u =y_i-x'_i\hat \beta_{2SLS}$$

Since you appear to using it, here is an implementation of a VCOV function in R code which can be used in lmtest::coeftest. For Stata users there is another relevant answer.

vcov2sls <- function(s1, s2) {
  X <- model.matrix(s2)
  dims <- dim(X)
  n <- dims[1]; p <- dims[2]
  mf <- model.frame(s1)
  y <- model.response(model.frame(s2))
  aux <- X  ## auxilliary model matrix
  aux[, 2L] <- as.numeric(model.response(mf))[seq_len(n)]
  b <- s2$coefficients
  r <- as.vector(y - aux %*% b)  ## residuals
  rss <- sum(r^2)
  rmse2 <- sqrt(mean(s2$residuals^2))  ## RMSE 2nd stage
  Rb <- vcov(s2)  ## biased vcov 2nd stage
  rmse <- sqrt(rss/n)
  corr <- (rmse/rmse2)^2
  R <- corr*Rb  ## corrected vcov
  return(R)
}

Here first the application with OLS models in both stages.

s1 <- lm(taxrate ~ votewon + industry + size + urbanisation + vote, data=DF)
yhat <- s1$fitted.values
s2 <- lm(sales ~ yhat + industry + size + urbanisation + vote, data=DF)

lmtest::coeftest(s2, vcov.=vcov2sls(s1, s2))
# t test of coefficients:
#   
#                Estimate Std. Error t value Pr(>|t|)
# (Intercept)   -18.45116   58.62518 -0.3147   0.7537
# yhat            1.57784    2.56770  0.6145   0.5405
# industryE       0.98174    4.81772  0.2038   0.8390
# industryF       2.09036    6.84134  0.3055   0.7607
# size2          -8.85327   11.73072 -0.7547   0.4524
# size3          -5.74011    6.74505 -0.8510   0.3970
# size4         -10.79326   12.40129 -0.8703   0.3865
# size5          -3.38280    5.14804 -0.6571   0.5128
# urbanisationB  -1.74588    5.98215 -0.2918   0.7711
# urbanisationC  -2.00370    6.11825 -0.3275   0.7441
# vote1          -1.01661    6.12665 -0.1659   0.8686

Identical with AER::ivreg:

fit1 <- AER::ivreg(sales ~ taxrate + industry + size + urbanisation + vote |
                     votewon + industry + size + urbanisation + vote, data=DF)
cbind(summary(fit1)$coe)
#                  Estimate Std. Error    t value  Pr(>|t|)
# (Intercept)   -18.4511636  58.625184 -0.3147310 0.7537026
# taxrate         1.5778430   2.567702  0.6144961 0.5404553
# industryE       0.9817432   4.817719  0.2037776 0.8389925
# industryF       2.0903627   6.841342  0.3055486 0.7606618
# size2          -8.8532690  11.730721 -0.7547080 0.4524166
# size3          -5.7401121   6.745046 -0.8510115 0.3970468
# size4         -10.7932611  12.401290 -0.8703337 0.3864595
# size5          -3.3828028   5.148043 -0.6571047 0.5128093
# urbanisationB  -1.7458832   5.982153 -0.2918486 0.7710816
# urbanisationC  -2.0037006   6.118249 -0.3274958 0.7440620
# vote1          -1.0166078   6.126650 -0.1659321 0.8685868

Now we could calculate a tobit model in the first stage. (I use AER::tobit since censReg doesn't seem to provide the fitted values.)

s1.tobit <- AER::tobit(taxrate ~ votewon + industry + size + urbanisation + vote,
                       left=12, right=33, data=DF)
yhat <- fitted(s1.tobit)
s2.tobit <- lm(sales ~ yhat + industry + size + urbanisation + vote, data=DF)

lmtest::coeftest(s2.tobit, vcov.=vcov2sls(s1.tobit, s2.tobit))
# t test of coefficients:
#   
#               Estimate Std. Error t value Pr(>|t|)
# (Intercept)   -6.71785   35.49067 -0.1893   0.8503
# yhat           1.08249    1.57840  0.6858   0.4946
# industryE      0.40529    3.77675  0.1073   0.9148
# industryF      1.18311    5.03625  0.2349   0.8148
# size2         -7.28430    8.43989 -0.8631   0.3904
# size3         -5.38410    5.78586 -0.9306   0.3546
# size4         -9.17192    9.06884 -1.0114   0.3146
# size5         -3.95675    4.27762 -0.9250   0.3575
# urbanisationB -2.50969    4.57805 -0.5482   0.5849
# urbanisationC -2.66487    4.80591 -0.5545   0.5806
# vote1         -0.60525    4.97518 -0.1217   0.9034

There might be objections to using a Tobit model in the first stage, but I am not aware of any.


Data:

DF <- structure(list(country = c("C", "C", "C", "C", "J", "J", "B", 
"B", "F", "F", "E", "E", "D", "D", "F", "F", "I", "I", "J", "J", 
"E", "E", "C", "C", "I", "I", "I", "I", "I", "I", "C", "C", "H", 
"H", "J", "J", "G", "G", "J", "J", "I", "I", "C", "C", "D", "D", 
"A", "A", "G", "G", "E", "E", "J", "J", "G", "G", "I", "I", "I", 
"I", "J", "J", "G", "G", "E", "E", "G", "G", "E", "E", "F", "F", 
"I", "I", "B", "B", "E", "E", "H", "H", "B", "B", "A", "A", "I", 
"I", "I", "I", "F", "F", "E", "E", "I", "I", "J", "J", "D", "D", 
"F", "F"), year = c(2005, 2010, 2010, 2005, 2005, 2010, 2010, 
2005, 2010, 2005, 2005, 2010, 2010, 2005, 2005, 2010, 2005, 2010, 
2005, 2010, 2010, 2005, 2010, 2005, 2005, 2010, 2005, 2010, 2010, 
2005, 2010, 2005, 2005, 2010, 2010, 2005, 2005, 2010, 2005, 2010, 
2005, 2010, 2005, 2010, 2010, 2005, 2005, 2010, 2010, 2005, 2010, 
2005, 2010, 2005, 2010, 2005, 2010, 2005, 2010, 2005, 2010, 2005, 
2010, 2005, 2010, 2005, 2010, 2005, 2005, 2010, 2005, 2010, 2005, 
2010, 2005, 2010, 2005, 2010, 2005, 2010, 2010, 2005, 2005, 2010, 
2005, 2010, 2010, 2005, 2010, 2005, 2010, 2005, 2005, 2010, 2005, 
2010, 2010, 2005, 2010, 2005), sales = c(15.48, 12.39, 3.72, 
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72, 
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72, 
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72, 
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72, 
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72, 
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72, 
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72, 
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72, 
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72, 
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9), industry = c("D", 
"D", "E", "E", "F", "F", "F", "F", "D", "D", "E", "E", "D", "D", 
"E", "E", "F", "F", "F", "F", "D", "D", "F", "F", "E", "E", "D", 
"D", "D", "D", "E", "E", "F", "F", "D", "D", "E", "E", "E", "E", 
"D", "D", "E", "E", "D", "D", "D", "D", "E", "E", "D", "D", "F", 
"F", "D", "D", "D", "D", "E", "E", "D", "D", "E", "E", "D", "D", 
"D", "D", "D", "D", "F", "F", "F", "F", "E", "E", "D", "D", "E", 
"E", "F", "F", "E", "E", "F", "F", "E", "E", "F", "F", "D", "D", 
"D", "D", "D", "D", "D", "D", "F", "F"), urbanisation = c("B", 
"B", "A", "A", "B", "B", "A", "A", "C", "C", "C", "C", "A", "A", 
"B", "B", "C", "C", "A", "A", "C", "C", "B", "B", "A", "A", "A", 
"A", "A", "A", "A", "A", "A", "A", "C", "C", "B", "B", "B", "B", 
"B", "B", "C", "C", "A", "A", "B", "B", "B", "B", "A", "A", "B", 
"B", "A", "A", "A", "A", "B", "B", "C", "C", "A", "A", "C", "C", 
"A", "A", "B", "B", "A", "A", "B", "B", "B", "B", "B", "B", "C", 
"C", "A", "A", "A", "A", "A", "A", "A", "A", "C", "C", "A", "A", 
"B", "B", "A", "A", "B", "B", "B", "B"), size = c(1, 1, 5, 5, 
5, 5, 1, 1, 1, 1, 5, 5, 5, 5, 2, 2, 2, 2, 5, 5, 1, 1, 1, 1, 5, 
5, 5, 5, 4, 4, 5, 5, 5, 5, 4, 4, 2, 2, 5, 5, 1, 1, 1, 1, 2, 2, 
1, 1, 2, 2, 5, 5, 1, 1, 3, 3, 2, 2, 2, 2, 5, 5, 4, 4, 1, 1, 5, 
5, 2, 2, 5, 5, 2, 2, 2, 2, 4, 4, 3, 3, 4, 4, 3, 3, 3, 3, 3, 3, 
5, 5, 3, 3, 2, 2, 3, 3, 1, 1, 5, 5), base_rate = c(14L, 14L, 
14L, 14L, 19L, 19L, 30L, 30L, 20L, 20L, 29L, 29L, 20L, 20L, 20L, 
20L, 24L, 24L, 19L, 19L, 29L, 29L, 14L, 14L, 24L, 24L, 24L, 24L, 
24L, 24L, 14L, 14L, 17L, 17L, 19L, 19L, 33L, 33L, 19L, 19L, 24L, 
24L, 14L, 14L, 20L, 20L, 23L, 23L, 33L, 33L, 29L, 29L, 19L, 19L, 
33L, 33L, 24L, 24L, 24L, 24L, 19L, 19L, 33L, 33L, 29L, 29L, 33L, 
33L, 29L, 29L, 20L, 20L, 24L, 24L, 30L, 30L, 29L, 29L, 17L, 17L, 
30L, 30L, 23L, 23L, 24L, 24L, 24L, 24L, 20L, 20L, 29L, 29L, 24L, 
24L, 19L, 19L, 20L, 20L, 20L, 20L), taxrate = c(12L, 14L, 14L, 
12L, 21L, 18L, 30L, 30L, 20L, 20L, 29L, 30L, 20L, 20L, 20L, 20L, 
24L, 24L, 21L, 18L, 30L, 29L, 14L, 12L, 24L, 24L, 24L, 24L, 24L, 
24L, 14L, 12L, 18L, 19L, 18L, 21L, 33L, 32L, 21L, 18L, 24L, 24L, 
12L, 14L, 20L, 20L, 22L, 25L, 32L, 33L, 30L, 29L, 18L, 21L, 32L, 
33L, 24L, 24L, 24L, 24L, 18L, 21L, 32L, 33L, 30L, 29L, 32L, 33L, 
29L, 30L, 20L, 20L, 24L, 24L, 30L, 30L, 29L, 30L, 18L, 19L, 30L, 
30L, 22L, 25L, 24L, 24L, 24L, 24L, 20L, 20L, 30L, 29L, 24L, 24L, 
21L, 18L, 20L, 20L, 20L, 20L), vote = c(0, 0, 0, 0, 1, 1, 1, 
0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 
1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 
1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 
1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 
1, 0, 1, 1, 1, 1, 0, 1, 1), votewon = c(0, 0, 0, 0, 1, 0, 1, 
0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 
1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 
0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 
1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 
0, 0, 1, 1, 0, 1, 0, 1, 1)), class = "data.frame", row.names = c(NA, 
-100L))

## convert variables to factors beforehand
DF[c(1, 2, 4, 5, 6, 9, 10)] <- lapply(DF[c(1, 2, 4, 5, 6, 9, 10)], factor)
$\endgroup$
27
  • 1
    $\begingroup$ Thank you so much! This is absolutely amazing. If there are any objections to the first stage being a tobit model I guess I can always report the normal and tobit-2SLS alongside each other. Once again, thank you very much, I am very grateful for it! $\endgroup$
    – Tom
    Commented Oct 10, 2020 at 17:59
  • $\begingroup$ Hi jay.sf. I have a small question which is a bit of a long shot.. I was running your code on my actual data, and I noticed that the line DATA <- as.data.frame(model.matrix(phantom ~ ., transform(data, phantom=0))) from vcov2sls, gives me the Error in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]) : contrasts can be applied only to factors with 2 or more levels error. Do you, from the top of your head have any idea why it might show up there, while s1.tobit and s2.tobit do not have that problem? $\endgroup$
    – Tom
    Commented Oct 11, 2020 at 9:15
  • $\begingroup$ @TomKisters Your issue seems to be caused by factor variables that have no variance. I could replicate the problem by doing DF$vote <- factor(1) befor running the regressions or your cited line. It might also be possible that they're all NA, i.e. DF$vote <- factor(NA). Hope that helps. Cheers. $\endgroup$
    – jay.sf
    Commented Oct 11, 2020 at 13:55
  • 1
    $\begingroup$ @TomKisters I'm delighted:) Unfortunately I don't know a source. However, you could try to censor the predictions using survival::Surv(time=ifelse(yhat >= 33, 33, ifelse(yhat <= 12, 12, yhat)), time2=33, event=ifelse(yhat >= 33, 0, ifelse(yhat <= 12, 2, 1)), type="interval"), I'm not sure how to make vcov2sls work then, though. AER::ivreg can be considered as a convenience function, because S.taxrate <- with(DF, survival::Surv(time=ifelse(taxrate >= 33, 33, ifelse(taxrate <= 12, 12, taxrate)), time2=33, event=ifelse(taxrate >= 33, 0, ifelse(taxrate <= 12, 2, 1)), type="interval")); ... $\endgroup$
    – jay.sf
    Commented Jan 6, 2021 at 10:45
  • 1
    $\begingroup$ @TomKisters ... survival::survreg(S.taxrate ~ votewon + industry + size + urbanisation + vote, data=DF, dist="gaussian") gives an identical result. May this help! $\endgroup$
    – jay.sf
    Commented Jan 6, 2021 at 10:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.