From 4052885e9db04923b71b434030f9bf01ae76b4d2 Mon Sep 17 00:00:00 2001 From: Jonathan Roth Date: Tue, 16 Jan 2024 13:00:06 -0500 Subject: [PATCH 1/2] fixed bug in deltarm with 1 postperiod where we werent' constructing A right --- R/deltarm.R | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/R/deltarm.R b/R/deltarm.R index 3ba1164..f5d3a57 100644 --- a/R/deltarm.R +++ b/R/deltarm.R @@ -91,7 +91,7 @@ dir_RM = base::rep("<=", base::length(d_RM)) # Add equality constraint for pre-period coefficients - prePeriodEqualityMat = base::cbind(base::diag(numPrePeriods), + prePeriodEqualityMat = base::cbind(base::diag(numPrePeriods), base::matrix(data = 0, nrow = numPrePeriods, ncol = numPostPeriods)) A_RM_s = base::rbind(A_RM_s, prePeriodEqualityMat) d_RM = base::c(d_RM, trueBeta[1:numPrePeriods]) @@ -203,10 +203,18 @@ d_RM = .create_d_RM(numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods) # If only use post period moments, construct indices for the post period moments only. - if (postPeriodMomentsOnly & numPostPeriods > 1){ - postPeriodIndices <- (numPrePeriods +1):base::NCOL(A_RM_s) - postPeriodRows <- base::which( base::rowSums( A_RM_s[ , postPeriodIndices] != 0 ) > 0 ) - rowsForARP <- postPeriodRows + if (postPeriodMomentsOnly){ + if(numPostPeriods > 1){ + postPeriodIndices <- (numPrePeriods +1):base::NCOL(A_RM_s) + postPeriodRows <- base::which( base::rowSums( A_RM_s[ , postPeriodIndices] != 0 ) > 0 ) + rowsForARP <- postPeriodRows + }else{ + #If only one post-period, then it is the last column + postPeriodRows <- base::which(A_RM_s[ ,NCOL(A_RM_s)] != 0 ) + A_RM_s <- A_RM_s[postPeriodRows, ] + d_RM <- d_RM[postPeriodRows] + + } } else{ rowsForARP <- 1:base::NROW(A_RM_s) } From 31b7bb34cc92d21e127b53e3f519edeb9eb43fb3 Mon Sep 17 00:00:00 2001 From: Jonathan Roth Date: Tue, 16 Jan 2024 13:38:47 -0500 Subject: [PATCH 2/2] updated other rm functions to fix issue with 1 post period --- R/deltarmb.R | 15 ++++++--- R/deltarmm.R | 88 +++++++++++++++++++++++++++----------------------- R/deltasdrm.R | 17 +++++++--- R/deltasdrmb.R | 16 ++++++--- R/deltasdrmm.R | 17 +++++++--- 5 files changed, 97 insertions(+), 56 deletions(-) diff --git a/R/deltarmb.R b/R/deltarmb.R index 953dd68..ae6a8e2 100644 --- a/R/deltarmb.R +++ b/R/deltarmb.R @@ -219,10 +219,17 @@ d_RMB = .create_d_RMB(numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods) # If only use post period moments, construct indices for the post period moments only. - if (postPeriodMomentsOnly & numPostPeriods > 1){ - postPeriodIndices <- (numPrePeriods +1):base::NCOL(A_RMB_s) - postPeriodRows <- base::which( base::rowSums( A_RMB_s[ , postPeriodIndices] != 0 ) > 0 ) - rowsForARP <- postPeriodRows + if (postPeriodMomentsOnly){ + if(numPostPeriods > 1){ + postPeriodIndices <- (numPrePeriods +1):base::NCOL(A_RMB_s) + postPeriodRows <- base::which( base::rowSums( A_RMB_s[ , postPeriodIndices] != 0 ) > 0 ) + rowsForARP <- postPeriodRows + }else{ + #If only one post-period, then it is the last column + postPeriodRows <- base::which(A_RMB_s[ ,NCOL(A_RMB_s)] != 0 ) + A_RMB_s <- A_RMB_s[postPeriodRows, ] + d_RMB <- d_RMB[postPeriodRows] + } } else{ rowsForARP <- 1:base::NROW(A_RMB_s) } diff --git a/R/deltarmm.R b/R/deltarmm.R index c95cacb..1c4bdcb 100644 --- a/R/deltarmm.R +++ b/R/deltarmm.R @@ -20,7 +20,7 @@ # numPrePeriods = number of pre-periods. This is an element of resultsObjects. # numPostPeriods = number of post-periods. This is an element of resultsObjects. # monotonicityDirection = "increasing" or "decreasing" - + # First construct matrix Atilde that takes first diffrences -- (numPrePeriods+numPostPeriods) x (numPrePeriods+numPostPeriods+1) # Note Atilde is just the positive moments; is not related to Atilde, the rotate matrix, in the paper # Note: Atilde initially includes t = 0. We then drop it. @@ -28,30 +28,30 @@ for (r in 1:(numPrePeriods+numPostPeriods)) { Atilde[r, r:(r+1)] = c(-1, 1) } - + # Create a vector to extract the max first dif, which corresponds with the first dif for period s, or minus this if max_positive == FALSE v_max_dif <- base::matrix(0, nrow = 1, ncol = numPrePeriods + numPostPeriods + 1) v_max_dif[(numPrePeriods+s):(numPrePeriods+1+s)] <- c(-1,1) - + if (max_positive == FALSE){ v_max_dif <- -v_max_dif } - + # The bounds for the first dif starting with period t are 1*v_max_dif if t<=0 and M*v_max_dif if t>0 A_UB <- base::rbind( pracma::repmat(v_max_dif, n=numPrePeriods, m = 1), pracma::repmat(Mbar*v_max_dif, n=numPostPeriods, m = 1)) - + # Construct A that imposes |Atilde * delta | <= A_UB * delta A = base::rbind(Atilde - A_UB, -Atilde - A_UB) - + # Remove all-zero rows of the matrix Atilde, corresponding with the constraint (delta_s - delta_s-1) - (delta_s - delta_s-1) <= (delta_s - delta_s-1) - (delta_s - delta_s-1) zerorows <- base::apply(X = A, MARGIN = 1, FUN = function(x) base::t(x) %*% x) <= 10^-10 A <- A[!zerorows, ] - + # Create matrix for shape restriction A_M <- .create_A_M(numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods, monotonicityDirection = monotonicityDirection, postPeriodMomentsOnly = FALSE) - + # Remove the period corresponding with t=0 if (dropZero) { A = A[, -(numPrePeriods+1)] @@ -72,7 +72,7 @@ # Inputs: # numPrePeriods = number of pre-periods. This is an element of resultsObjects. # numPostPeriods = number of post-periods. This is an element of resultsObjects. - + A_RM = .create_A_RM(numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods, Mbar = 0, s = 0, dropZero = dropZero) # d doesn't depend on Mbar or s; we just use this to get the dims right d_RM = base::rep(0, base::NROW(A_RM)) @@ -83,37 +83,37 @@ # DELTA^{RMM}(Mbar) Identified Set Helper Functions -------------------- .compute_IDset_DeltaRMM_fixedS <- function(s, Mbar, max_positive, - trueBeta, l_vec, numPrePeriods, numPostPeriods, + trueBeta, l_vec, numPrePeriods, numPostPeriods, monotonicityDirection) { # This helperf function computes the upper and lower bound of the identified set # given the event study coefficients, lvec and Mbar. It computes the identified # set for a user-specified choice of s and (+), (-). This is used by # the function compute_IDset_DeltaRMM below. - + # Create objective function: Wish to min/max l'delta_post fDelta = base::c(base::rep(0, numPrePeriods), l_vec) - + # Create A_RMM, d_RMM for this choice of s, max_positive A_RMM_s = .create_A_RMM(numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods, Mbar = Mbar, s = s, max_positive = max_positive, monotonicityDirection = monotonicityDirection) d_RMM = .create_d_RMM(numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods) - + # Create vector for direction of inequalities associated with RM dir_RMM = base::rep("<=", base::length(d_RMM)) - + # Add equality constraint for pre-period coefficients prePeriodEqualityMat = base::cbind(base::diag(numPrePeriods), base::matrix(data = 0, nrow = numPrePeriods, ncol = numPostPeriods)) A_RMM_s = base::rbind(A_RMM_s, prePeriodEqualityMat) d_RMM = base::c(d_RMM, trueBeta[1:numPrePeriods]) dir_RMM = base::c(dir_RMM, base::rep("==", base::NROW(prePeriodEqualityMat))) - + # Specify variables between (-inf, inf) bounds = base::list(lower = base::list(ind = 1:(numPrePeriods + numPostPeriods), val = base::rep(-Inf, numPrePeriods+numPostPeriods)), upper = base::list(ind = 1:(numPrePeriods + numPostPeriods), val = base::rep(Inf, numPrePeriods+numPostPeriods))) - + # Create and solve for max results.max = Rglpk::Rglpk_solve_LP(obj = fDelta, max = TRUE, @@ -121,7 +121,7 @@ dir = dir_RMM, rhs = d_RMM, bounds = bounds) - + # Create and solve for min results.min = Rglpk::Rglpk_solve_LP(obj = fDelta, max = FALSE, @@ -129,7 +129,7 @@ dir = dir_RMM, rhs = d_RMM, bounds = bounds) - + if (results.max$status != 0 & results.min$status != 0) { # If the solver does not return solution, we just return the l_vec'trueBeta. id.ub = (base::t(l_vec) %*% trueBeta[(numPrePeriods+1):(numPrePeriods+numPostPeriods)]) @@ -166,28 +166,28 @@ # dataframe with columns # id.ub = upper bound of ID set # id.lb = lower bound of ID set - + # Construct identified sets for (+) at each value of s min_s = -(numPrePeriods - 1) id_bounds_plus = purrr::map_dfr( .x = min_s:0, .f = ~.compute_IDset_DeltaRMM_fixedS(s = .x, Mbar = Mbar, max_positive = TRUE, trueBeta = trueBeta, l_vec = l_vec, - numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods, + numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods, monotonicityDirection = monotonicityDirection) ) id_bounds_minus = purrr::map_dfr( .x = min_s:0, .f = ~.compute_IDset_DeltaRMM_fixedS(s = .x, Mbar = Mbar, max_positive = FALSE, trueBeta = trueBeta, l_vec = l_vec, - numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods, + numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods, monotonicityDirection = monotonicityDirection) ) - + # Construct the identified set by taking the max of the upper bound and the min of the lower bound id.lb = base::min(base::min(id_bounds_plus$id.lb), base::min(id_bounds_minus$id.lb)) id.ub = base::max(base::max(id_bounds_plus$id.ub), base::max(id_bounds_minus$id.ub)) - + # Return identified set base::return(tibble::tibble( id.lb = id.lb, @@ -204,29 +204,37 @@ # This function computes the ARP CI that includes nuisance parameters # for Delta^{RMM}(Mbar) for a fixed s and (+),(-). This functions uses ARP_computeCI for all # of its computations. It is used as a helper function in computeConditionalCS_DeltaRM below. - + # Check that hybrid_flag equals LF or ARP if (hybrid_flag != "LF" & hybrid_flag != "ARP") { base::stop("hybrid_flag must equal 'ARP' or 'LF'") } - + # Create hybrid_list object hybrid_list = base::list(hybrid_kappa = hybrid_kappa) - + # Create matrix A_RMM_s, and vector d_RMM A_RMM_s = .create_A_RMM(numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods, Mbar = Mbar, s = s, max_positive = max_positive, monotonicityDirection = monotonicityDirection) d_RMM = .create_d_RMM(numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods) - + # If only use post period moments, construct indices for the post period moments only. - if (postPeriodMomentsOnly & numPostPeriods > 1){ - postPeriodIndices <- (numPrePeriods +1):base::NCOL(A_RMM_s) - postPeriodRows <- base::which( base::rowSums( A_RMM_s[ , postPeriodIndices] != 0 ) > 0 ) - rowsForARP <- postPeriodRows + if (postPeriodMomentsOnly){ + if(numPostPeriods > 1){ + postPeriodIndices <- (numPrePeriods +1):base::NCOL(A_RMM_s) + postPeriodRows <- base::which( base::rowSums( A_RMM_s[ , postPeriodIndices] != 0 ) > 0 ) + rowsForARP <- postPeriodRows + }else{ + #If only one post-period, then it is the last column + postPeriodRows <- base::which(A_RMM_s[ ,NCOL(A_RMM_s)] != 0 ) + A_RMM_s <- A_RMM_s[postPeriodRows, ] + d_RMM <- d_RMM[postPeriodRows] + + } } else{ rowsForARP <- 1:base::NROW(A_RMM_s) } - + # if there is only one post-period, we use the no-nuisance parameter functions if (numPostPeriods == 1) { if (hybrid_flag == "LF") { @@ -279,17 +287,17 @@ computeConditionalCS_DeltaRMM <- function(betahat, sigma, numPrePeriods, numPost # # Outputs: # data_frame containing upper and lower bounds of CI. - + # Create minimal s index for looping. min_s = -(numPrePeriods - 1) s_indices = min_s:0 - + # If grid.ub, grid.lb is not specified, we set these bounds to be equal to the id set under parallel trends # {0} +- 20*sdTheta (i.e. [-20*sdTheta, 20*sdTheta]. sdTheta <- base::c(base::sqrt(base::t(l_vec) %*% sigma[(numPrePeriods+1):(numPrePeriods+numPostPeriods), (numPrePeriods+1):(numPrePeriods+numPostPeriods)] %*% l_vec)) if (base::is.na(grid.ub)) { grid.ub = 20*sdTheta } if (base::is.na(grid.lb)) { grid.lb = -20*sdTheta } - + # Loop over s values for (+), (-), left join the resulting CIs based on the grid CIs_RMM_plus_allS = base::matrix(0, nrow = gridPoints, ncol = base::length(s_indices)) CIs_RMM_minus_allS = base::matrix(0, nrow = gridPoints, ncol = base::length(s_indices)) @@ -301,8 +309,8 @@ computeConditionalCS_DeltaRMM <- function(betahat, sigma, numPrePeriods, numPost alpha = alpha, hybrid_flag = hybrid_flag, hybrid_kappa = hybrid_kappa, postPeriodMomentsOnly = postPeriodMomentsOnly, monotonicityDirection = monotonicityDirection, gridPoints = gridPoints, grid.ub = grid.ub, grid.lb = grid.lb) - CIs_RM_plus_allS[,s_i] = CI_s_plus$accept - + CIs_RMM_plus_allS[,s_i] = CI_s_plus$accept + # Compute CI for s, (-) and bind it to all CI's for (-) CI_s_minus = .computeConditionalCS_DeltaRMM_fixedS(s = s_indices[s_i], max_positive = FALSE, Mbar = Mbar, betahat = betahat, sigma = sigma, numPrePeriods = numPrePeriods, @@ -312,14 +320,14 @@ computeConditionalCS_DeltaRMM <- function(betahat, sigma, numPrePeriods, numPost gridPoints = gridPoints, grid.ub = grid.ub, grid.lb = grid.lb) CIs_RMM_minus_allS[,s_i] = CI_s_minus$accept } - + CIs_RMM_plus_maxS = base::apply(CIs_RMM_plus_allS, MARGIN = 1, FUN = base::max) CIs_RMM_minus_maxS = base::apply(CIs_RMM_minus_allS, MARGIN = 1, FUN = base::max) - + # Take the max between (+), (-) and Construct grid containing theta points and whether any CI accepted CI_RMM = tibble::tibble(grid = base::seq(grid.lb, grid.ub, length.out = gridPoints), accept = base::pmax(CIs_RMM_plus_maxS, CIs_RMM_minus_maxS)) - + # Compute length if returnLength == TRUE, else return grid if (returnLength) { gridLength <- 0.5 * ( base::c(0, base::diff(CI_RMM$grid)) + base::c(base::diff(CI_RMM$grid), 0 ) ) diff --git a/R/deltasdrm.R b/R/deltasdrm.R index 8f31942..01f125d 100644 --- a/R/deltasdrm.R +++ b/R/deltasdrm.R @@ -199,14 +199,23 @@ d_SDRM = .create_d_SDRM(numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods) # If only use post period moments, construct indices for the post period moments only. - if (postPeriodMomentsOnly & numPostPeriods > 1){ - postPeriodIndices <- (numPrePeriods +1):base::NCOL(A_SDRM_s) - postPeriodRows <- base::which( base::rowSums( A_SDRM_s[ , postPeriodIndices] != 0 ) > 0 ) - rowsForARP <- postPeriodRows + if (postPeriodMomentsOnly){ + if(numPostPeriods > 1){ + postPeriodIndices <- (numPrePeriods +1):base::NCOL(A_SDRM_s) + postPeriodRows <- base::which( base::rowSums( A_SDRM_s[ , postPeriodIndices] != 0 ) > 0 ) + rowsForARP <- postPeriodRows + }else{ + #If only one post-period, then it is the last column + postPeriodRows <- base::which(A_SDRM_s[ ,NCOL(A_SDRM_s)] != 0 ) + A_SDRM_s <- A_SDRM_s[postPeriodRows, ] + d_SDRM <- d_SDRM[postPeriodRows] + + } } else{ rowsForARP <- 1:base::NROW(A_SDRM_s) } + # if there is only one post-period, we use the no-nuisance parameter functions if (numPostPeriods == 1) { if (hybrid_flag == "LF") { diff --git a/R/deltasdrmb.R b/R/deltasdrmb.R index e0073c5..c1653e9 100644 --- a/R/deltasdrmb.R +++ b/R/deltasdrmb.R @@ -213,10 +213,18 @@ d_SDRMB = .create_d_SDRMB(numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods) # If only use post period moments, construct indices for the post period moments only. - if (postPeriodMomentsOnly & numPostPeriods > 1){ - postPeriodIndices <- (numPrePeriods +1):base::NCOL(A_SDRMB_s) - postPeriodRows <- base::which( base::rowSums( A_SDRMB_s[ , postPeriodIndices] != 0 ) > 0 ) - rowsForARP <- postPeriodRows + if (postPeriodMomentsOnly){ + if(numPostPeriods > 1){ + postPeriodIndices <- (numPrePeriods +1):base::NCOL(A_SDRMB_s) + postPeriodRows <- base::which( base::rowSums( A_SDRMB_s[ , postPeriodIndices] != 0 ) > 0 ) + rowsForARP <- postPeriodRows + }else{ + #If only one post-period, then it is the last column + postPeriodRows <- base::which(A_SDRMB_s[ ,NCOL(A_SDRMB_s)] != 0 ) + A_SDRMB_s <- A_SDRMB_s[postPeriodRows, ] + d_SDRMB <- d_SDRMB[postPeriodRows] + + } } else{ rowsForARP <- 1:base::NROW(A_SDRMB_s) } diff --git a/R/deltasdrmm.R b/R/deltasdrmm.R index 77fb8a4..a6b7eca 100644 --- a/R/deltasdrmm.R +++ b/R/deltasdrmm.R @@ -215,14 +215,23 @@ d_SDRMM = .create_d_SDRMM(numPrePeriods = numPrePeriods, numPostPeriods = numPostPeriods) # If only use post period moments, construct indices for the post period moments only. - if (postPeriodMomentsOnly & numPostPeriods > 1){ - postPeriodIndices <- (numPrePeriods +1):base::NCOL(A_SDRMM_s) - postPeriodRows <- base::which( base::rowSums( A_SDRMM_s[ , postPeriodIndices] != 0 ) > 0 ) - rowsForARP <- postPeriodRows + if (postPeriodMomentsOnly){ + if(numPostPeriods > 1){ + postPeriodIndices <- (numPrePeriods +1):base::NCOL(A_SDRMM_s) + postPeriodRows <- base::which( base::rowSums( A_SDRMM_s[ , postPeriodIndices] != 0 ) > 0 ) + rowsForARP <- postPeriodRows + }else{ + #If only one post-period, then it is the last column + postPeriodRows <- base::which(A_SDRMM_s[ ,NCOL(A_SDRMM_s)] != 0 ) + A_SDRMM_s <- A_SDRMM_s[postPeriodRows, ] + d_SDRMM <- d_SDRMM[postPeriodRows] + + } } else{ rowsForARP <- 1:base::NROW(A_SDRMM_s) } + # if there is only one post-period, we use the no-nuisance parameter functions if (numPostPeriods == 1) { if (hybrid_flag == "LF") {