Toy example: Estimation differences between fect R and Stata

[rondobohrens]

Thank you for this great package!

I am currently collaborating with a coauthor on a project involving Matrix Completion (MC). I work in R, and he works in Stata.
We are trying to reproduce results across both languages, but fail to do so.

Using the simulation data and examples you provide here under Matrix Completion (MC), I am able to reproduce your exact results.

out.mc <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"), force = "two-way", method = "mc", CV = TRUE, se = TRUE, nboots = 200, parallel = TRUE)

Call:
fect.formula(formula = Y ~ D + X1 + X2, data = simdata, index = c("id", 
    "time"), force = "two-way", CV = TRUE, method = "mc", se = TRUE, 
    nboots = 200, parallel = TRUE)

ATT:
                            ATT   S.E. CI.lower CI.upper p.value
Tr obs equally weighted   4.262 0.2933    3.687    4.837       0
Tr units equally weighted 2.826 0.2290    2.377    3.275       0

Covariates:
    Coef    S.E. CI.lower CI.upper p.value
X1 1.017 0.03083   0.9568    1.078       0
X2 2.959 0.03027   2.8999    3.019       0

However, trying to run this in Stata, we get a different ATT estimate, but the same results for covariates.

use "/fect_simdata.dta", clear

fect Y, treat (D) unit(id) time(time) cov(X1 X2) method("mc") force("two-way") cv nboots(200) se 

mat list e(ATT)
mat list e(coefs)

produces

. fect Y, treat (D) unit(id) time(time) cov(X1 X2) method("mc") force("two-way") cv nboots(200) se 
Balanced Panel Data
Treatment has reversals.
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Cross Validation...
mc: lambda=.0131 lambda.norm=1 mspe=7.1809
mc: lambda=.0055 lambda.norm=.4217 mspe=5.5554
mc: lambda=.0023 lambda.norm=.1778 mspe=5.304
mc: lambda=.001 lambda.norm=.075 mspe=5.4196
mc: lambda=.0004 lambda.norm=.0316 mspe=5.5083
mc: lambda=.0002 lambda.norm=.0133 mspe=5.8303
mc: lambda=.0001 lambda.norm=.0056 mspe=7.1136
mc: lambda=0 lambda.norm=.0024 mspe=7.1522
mc: lambda=0 lambda.norm=.001 mspe=7.1687
mc: lambda=0 lambda.norm=.0004 mspe=7.1757
optimal lambda=.0023, lambda.norm=.1778 in mc model
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Bootstrapping...
ATT Estimation: Already Bootstrapped 100 Times
ATT Estimation: Already Bootstrapped 200 Times


e(ATT)[1,6]
            ATT            N           sd  Lower_Bound  Upper_Bound       pvalue
r1    4.0790949         1503     .3387551    3.4314052    4.8718991            0

. mat list e(coefs)

e(coefs)[3,5]
             coef           sd      p_value  lower_bound  upper_bound
cons    11.969906    .28209105            0    11.430512    12.477093
  X1    1.0170004    .02963712            0    .95439136    1.0690237
  X2    2.9594574    .03025816            0    2.8992956    3.0156198

If we force them to use the exact same optimal lambda (as identified in R, lambda.cv = 0.002334298), we still get different estimates for the main effect.

fect Y, treat (D) unit(id) time(time) cov(X1 X2) method("mc") force("two-way") cv nboots(200) se lambda(0.002334298)
mat list e(ATT)
mat list e(coefs)

Balanced Panel Data
Treatment has reversals.
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Cross Validation...
mc: lambda=.0023 lambda.norm=.1778 mspe=5.304
optimal lambda=.0023, lambda.norm=.1778 in mc model
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Bootstrapping...
ATT Estimation: Already Bootstrapped 100 Times
ATT Estimation: Already Bootstrapped 200 Times


e(ATT)[1,6]
            ATT            N           sd  Lower_Bound  Upper_Bound       pvalue
r1    4.0790949         1503    .29853074    3.5759497    4.7775187            0

. mat list e(coefs)

e(coefs)[3,5]
             coef           sd      p_value  lower_bound  upper_bound
cons    11.969906    .29753646            0    11.422481    12.554377
  X1    1.0170004     .0304669            0    .96010947    1.0780883
  X2    2.9594574    .02831529            0    2.9031889    3.0128927

I can produce the example in Stata provided here without any issues.

Can you provide us with any pointers to why this is?

[rondobohrens]

I think we found the root cause for this. The default value for the tolerancetol is different.
In R, it is 0.001 and in Stata it is 0.00001.
We started to see very close agreements with R and Stata when we set 1e-6.

We lost a bit of sanity along the way - but I think it would be good to set those values equal in both languages.