penalized synthetic control estimation (Abadie & L'Hour, 2021)

penalized_synthetic_control.R

#' Penalized synthetic control estimator
#' 
#' Estimate synthetic control with penalization
#' according to Abadie & L'Hour.
#'
#' @param X1 treated unit covariates
#' @param X0 donor units covariates
#' @param v variable weights
#' @param lambda penalization parameter
#' @param ... osqp settings using osqp::osqpSettings()
#' 
#' @return list of weights and the osqp optimizer object
#' 
#' @export
synthetic_control <- function(X1, X0, v, lambda = 0, opt_pars = osqp::osqpSettings(polish = TRUE)) {
  N_donors <- ncol(X0)
  X0v <- X0*sqrt(v)
  X1v <- X1*sqrt(v)
  
  # components for quadratic program
  # see https://github.com/jeremylhour/pensynth/blob/master/functions/wsoll1.R
  X0VX0 <- crossprod(X0v)
  X1VX0 <- crossprod(X1v, X0v)
  Delta <- apply(X0v - c(X1v), 2, crossprod)
  
  # linear constraint matrix
  Amat   <- rbind(rep(1, N_donors), diag(N_donors))
  lbound <- c(1, rep(0, N_donors))
  ubound <- rep(1, N_donors + 1)
  
  # stop annoying printing with capture.output.
  o <- capture.output({ 
    solver <- osqp::osqp(
      P = X0VX0,
      q = -X1VX0 + lambda*Delta,
      A = Amat,
      l = lbound,
      u = ubound,
      pars = opt_pars
    )
    result <- solver$Solve()
  })
  
  return(list(w = result$x, solution = result, optim = solver))
}

X0 <- matrix(
  c(1, 1.3, 
    0.5, 1.8, 
    1.1, 2.4, 
    1.8, 1.8, 
    1.3, 1.8), 2)
X1 <- matrix(c(0.8, 1.65), 2)

lam_seq <- seq(0, 12, length.out = 100)
res <- sapply(lam_seq, \(l) synthetic_control(X1, X0, rep(1, 2), lambda = l)$w)
round(res, 3)

library(tidyverse)
data.frame(lam = lam_seq, w = t(res)) |> 
  filter(lam != 0) |> 
  pivot_longer(-lam, names_prefix = "w.") |> 
  ggplot(aes(x = lam, y = value, colour = name)) +
  geom_line() +
  scale_x_continuous(trans = "log10") +
  labs(
    x = "Penalty parameter", 
    y = "Weight", 
    colour = "Donor unit",
    title = "Weights path as a function of lambda"
  ) +
  theme_minimal()

solpath <- data.frame(lam = lam_seq, x = t(X0%*%res))

data.frame(id = c("treated", rep("donor", ncol(X0))), x = t(cbind(X1, X0))) |> 
  ggplot(aes(x = x.1, y = x.2)) +
  geom_line(data = solpath, col = "grey", linetype = 1) +
  geom_point(aes(fill = id), size = 3, pch = 21, col = "white") +
  theme_minimal() +
  xlim(0, 2) + 
  ylim(1, 2.75)