changed solver options --- ECOS is the only option now, as SCS had

trouble enforcing constraints and OSQP didn't solve very accurately
synth-inference · Nov 21, 2019 · 68ccb22 · 68ccb22
1 parent 1cae5ee
commit 68ccb22
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 37 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -10,6 +10,6 @@ Description: A package for estimating average treatment effects in panel data.
     synthdid currently provides methods only for the case that all treated units
     adopt treatment at the same time. This package is currently in beta.
 Depends: 
-    R (>= 3.4.4), osqp 
+    R (>= 3.4.4), CVXR
 License: none
 RoxygenNote: 6.1.1.9000
diff --git a/R/sdid_lib.R b/R/sdid_lib.R
@@ -10,47 +10,23 @@
 #' @param intercept. Defaults to FALSE
 #' @return gamma
 #' @export sc_weight
-sc_weight = function(M, target, zeta = 1, intercept = FALSE, solver = c("OSQP", "ECOS", "SCS"), verbose = FALSE) {
+sc_weight = function(M, target, zeta = 1, intercept = FALSE, solver = c("ECOS"), verbose = FALSE) {
     solver = match.arg(solver)
     if (nrow(M) != length(target)) {
         stop("invalid dimensions")
     }
 
-    if(solver %in% c('ECOS', 'SCS')) {
-        weights = CVXR::Variable(ncol(M))
-        if(intercept) { 
-            theintercept = CVXR::Variable(1)
-            objective = zeta^2 * length(target) * sum(weights^2) + sum((M %*% weights - target - theintercept)^2)
-        } else {
-            objective = zeta^2 * length(target) * sum(weights^2) + sum((M %*% weights - target)^2)
-        }
-        constraints = list(sum(weights) == 1, weights >= 0)  
-        cvx.problem = CVXR::Problem(CVXR::Minimize(objective), constraints)
-        cvx.output = solve(cvx.problem, solver = solver, verbose = verbose)
-        as.numeric(cvx.output$getValue(weights))     
-    } else { 
-        if(intercept) {
-            Dmat = t(cbind(M,1)) %*% cbind(M,1) + zeta^2 * diag(rep(c(1,0),c(ncol(M), 1))) * length(target)
-            dvec = t(target) %*% cbind(M,1)
-            A = cbind(rbind(rep(1, ncol(M)),
-                             diag(1, ncol(M))), 0)
-            lb = c(1, rep(0, ncol(M)))
-            ub = c(1, rep(Inf, ncol(M)))
-            soln = osqp::solve_osqp(Dmat, -dvec, A, lb, ub, osqpSettings(verbose=FALSE, adaptive_rho=FALSE))
-            soln$x[1:ncol(M)]
-
-        } else {
-            Dmat = t(M) %*% M + zeta^2 * diag(ncol(M)) * length(target)
-            dvec = t(target) %*% M
-            meq = 1
-            A = rbind(rep(1, ncol(M)),
-                      diag(1, ncol(M)))
-            lb = c(1, rep(0, ncol(M)))
-            ub = c(1, rep(Inf, ncol(M)))
-            soln = osqp::solve_osqp(Dmat, -dvec, A, lb, ub, osqpSettings(verbose=FALSE, adaptive_rho=FALSE))
-            soln$x
-        }
-    } 
+    weights = CVXR::Variable(ncol(M))
+    if(intercept) { 
+        theintercept = CVXR::Variable(1)
+        objective = zeta^2 * length(target) * sum(weights^2) + sum((M %*% weights - target - theintercept)^2)
+    } else {
+        objective = zeta^2 * length(target) * sum(weights^2) + sum((M %*% weights - target)^2)
+    }
+    constraints = list(sum(weights) == 1, weights >= 0)  
+    cvx.problem = CVXR::Problem(CVXR::Minimize(objective), constraints)
+    cvx.output = solve(cvx.problem, solver = solver, verbose = verbose)
+    as.numeric(cvx.output$getValue(weights))     
 }