Adding univariate Imbens and Wager (2019) min-max optimal estimator

[evanmunro]

I have added an implementation of the Imbens and Wager min-max optimal estimator for a univariate running variable in the sharp RDD setting. I've tested the estimate and confidence interval calculation versus the R code for three different datastes and it matches. I also added a couple of examples to the README for the package.

There's a few things I would like to discuss:

• I still don't like my implementation of the discretization, but I was struggling to integrate it well into RunningVariable logic. We could discuss and potentially refactor here.
• I currently calculate the conditional variance of the outcome using linear regression and do not allow the user to specify it. I calculate a homoskedastic standard error estimate for the RDD estimate. Could add more options here ?
• It might be useful to investigate why Hypatia is so slow vs Mosek for this formulation of the problem when they have benchmarked vs Mosek on other problems. Is there a better way to formulate?
• Would it be worth having some kind of common return type for different SharpRDD estimators?

Once we have discussed and made some changes, I will:

• add docstrings for potential documentation
• add some tests

[codecov-io]

Codecov Report

Merging #13 (43fb7d9) into master (9764a26) will increase coverage by 11.25%. The diff coverage is 92.00%.

@@             Coverage Diff             @@
##           master      #13       +/-   ##
===========================================
+ Coverage   50.52%   61.78%   +11.25%     
===========================================
  Files           6        7        +1     
  Lines         287      348       +61     
===========================================
+ Hits          145      215       +70     
+ Misses        142      133        -9     

Impacted Files	Coverage Δ
src/RegressionDiscontinuity.jl	100.00% <ø> (ø)
src/running_variable.jl	40.44% <87.50%> (+11.11%)	:arrow_up:
src/minmax_optimal.jl	93.10% <93.10%> (ø)
src/local_linear.jl	52.63% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 9764a26...43fb7d9. Read the comment docs.

[nignatiadis]

Thanks a lot for this! The code looks good to me. Let me make some initial remarks and I will look at it more carefully later:

• I think Suppressor is not necessary here. Instead the output level can be controlled by the solver object, for example the following would suppress the Mosek output:

quiet_mosek = optimizer_with_attributes(Mosek.Optimizer, "QUIET" => true)
model = Model(quiet_mosek)
...

• Regarding variance estimation, in local_linear.jl I define the abstract type VarianceEstimator. By dispatching on that we should be able to allow different types of variance estimates. Eventually we could coordinate e.g. with https://github.com/matthieugomez/Vcov.jl.

• I think Hypatia is most efficient when things are specified directly in terms of cones. Here I think the relevant part would be to move the quadratic

  1/4*(sum(n[i]*G[i]^2/σ²[i] for i in 1:d)) + l1^2 

  into a Second order cone constraint

  @constraint(model, x in SecondOrderCone())

• Discretization is always a bit annoying.. So let me think a bit more about that and get back to you. I guess what's special about Sharp RD discretization is that you typically want the discretization to satisfy the following two properties: the cutoff chas to be a part of the grid, and the grid has to be faithful to whether treatment is assigned for Z_i>= c, or Z_i>c and so forth. There already is some code though for doing this so it would be great to reuse it.

• Thinking about it again, I don't like the Centered abstraction I introduced. Probably it's easier if RunningVariable just stores the centered Z as well?

[evanmunro]

Thanks for pointing out the code on the discretization. I hadn't noticed it but will integrate and modify the centering abstraction as well.

Making sure I understand how to explicitly show Hypatia that this is a SOCP. I would do something like

@constraint(model, [x; obj] in SecondOrderCone()) 

where obj is a d+1-length vector with obj[i] = 1/2*sqrt(n[i])*G[i]/sigma for i in 1:d and obj[j] = l1 for j=d+1.

But then in the objective would there not still be a quadratic? e.g. it would be

@objective(model, Min, x^2 - l2 + l3) 

New to SOCP so this may be a very basic question!

Edit: And on the variance estimate. I think that there are two issues here. The first is the conditional mean function variance estimation. And then the second is the variance estimation of the estimate tau, given the estimate of the conditional mean function. We would need abstractions to handle both, I think, right ?

[nignatiadis]

Yeah that was my idea (including keeping a quadratic below). Not fully sure this will speed Hypatia up but that would have been my first attempt to speed it up. Another thing that might work (though again not fully sure) is to not define G explicitly as variable but e.g. do something like (may need broadcasting)

 G = @expression(model, 2*B* f + l2 *(1-W) + l3*W + l4*(X - c) + l5*(W - 0.5)*(X-c))

[evanmunro]

The implementation of the discretized running variable is better now and I got rid of the centering abstraction. Let me know what you think.

That leaves a few more comments/things outstanding:

1. In this setting there are two variance estimates. One for the estimate of the variance of the conditional outcome, and one for the estimate of the treatment effect. For the first, I don't allow the user to specify so far. For the second, I have the homoskedastic estimator implemented the way you suggested.

2. Right now I don't have z and p-value columns in the CoefTable. Should I add that ? I wasn't sure how calculating z and p value interacted with the confidence interval construction that takes into account max bias.

3. Tests and docstrings

Let me know also any other comments you have.

[nignatiadis]

Looks great to me! I think we can postpone 1. for now and add it later if there is interest. 2. is also fine, I can take a stab at it after this PR has been merged, so I think 3. is the most important. Also perhaps we may want to rethink the name for MinMaxOptRD. Maybe just OptRDD? Or ImbensWagerOptRD? (Also in light of your other results, and also the Armstrong/Kolesar work, we could have something like MinMaxOptRD and then a slot with convexclass, e.g. here convexclass would be class of functions with uniformly bounded second derivative, but it could encode additional constraints and dispatch appropriately each time. We can revisit this later.)

By the way did the JuMP model specification changes help with Hypatia's timings?

[evanmunro]

Agreed with your suggestion on the naming for a later iteration of this if we end up adding support for a wider range of convex function classes. For now I have named it ImbensWagerOptRD and added docstrings for minmax_optimal.jl. Still missing docstrings from running_variable.jl since I had one question there. What is the convention for the documenting different constructors for RunningVariable? Do you just document each of them ?

I also added one test which compares the results from this package to a hard coded estimate, bias, and standard error from R. An alternative would be to use RCall to generate this dynamically - would that be better? I also add a confidence interval test to make sure the bias adjusted confidence interval construction is correct. Any ideas for other tests?

Hypatia timing did not meaningfully improve with the two suggested code changes, unfortunately. It is ~40s for the Lee data versus 0.1 seconds for Mosek.

[evanmunro]

One thing I now realize is the test errors on the CI check if MosekTools is not included as a dependency?

[nignatiadis]

Ah this looks great!!

Regarding MosekTools, I think we cannot run Mosek on the CI, so maybe use Hypatia there? (It does not matter if it takes a minute + on the CI). We can have custom dependencies for the "tests" too, basically by creating a Project.toml file within the "tests" folder (https://julialang.github.io/Pkg.jl/v1/creating-packages/#Adding-tests-to-the-package). A user downloading the package then would not be forced to install such dependencies.

I think hard-coding the R results sounds good. It would be helpful to commit the R code used to generate the hard-coded test values though (for reproducibility in the future).

Which constructors do you mean for RunningVariable? I think if they are very closely related then you could even document them in the same docstring, but separately could also work.

Another question: Is quartic_regression_bound the Armstrong-Kolesar proposal for a data-driven choice of B? Also what do you think about making an abstract type and concrete types representing quartic_regression_bound and estimate_second_deriv_bound, and allowing the B field of ImbensWagerOptRD to take on these values? (Similarly to how NaiveLocalLinearRD can take a bandwidth or a bandwidth selector.)

[evanmunro]

I decided to remove the data-driven bound on the second derivative for now. Imbens Wager doesn't suggest doing something like this (optrdd in R doesn't include this as an option). Since I'm working on a better of version of this separately as part of the more general convex constraints I'll wait before implementing some different options for selecting that bound.

I added some doc strings for RunningVariable. I added the R file to the tests folder.

With Hypatia added to the tests, are you ok restricting compatibility to Julia 1.5? That should get the CI checks passing.

[nignatiadis]

Yes sounds good to pin it at 1.5! Two more minor things: the readme example for NaiveLocalLinearRD is missing the fit command.

Also the indentation is off at some spots, can you e.g. run https://github.com/domluna/JuliaFormatter.jl on the files you changed? Otherwise this should be good to be merged!

[evanmunro]

Made those two small changes!

[nignatiadis]

Merged! Thanks @evanmunro!