Skip to content

Commit

Permalink
Merge pull request #50 from asheshrambachan/1period_issue
Browse files Browse the repository at this point in the history 
Fix bug in RM functions with 1 period
  • Loading branch information
@jonathandroth
jonathandroth authored Jan 23, 2024
2 parents 5773ea4 + 31b7bb3 commit bd629c1
Show file tree
Hide file tree
Showing 6 changed files with 110 additions and 61 deletions.
18 changes: 13 additions & 5 deletions R/deltarm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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])
Expand Down Expand Up @@ -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)
}
Expand Down
15 changes: 11 additions & 4 deletions R/deltarmb.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
}
Expand Down
88 changes: 48 additions & 40 deletions R/deltarmm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,38 +20,38 @@
# 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.
Atilde = base::matrix(0, nrow = numPrePeriods+numPostPeriods, ncol = numPrePeriods+numPostPeriods+1)
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)]
Expand All @@ -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))
Expand All @@ -83,53 +83,53 @@

# 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,
mat = A_RMM_s,
dir = dir_RMM,
rhs = d_RMM,
bounds = bounds)

# Create and solve for min
results.min = Rglpk::Rglpk_solve_LP(obj = fDelta,
max = FALSE,
mat = A_RMM_s,
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)])
Expand Down Expand Up @@ -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,
Expand All @@ -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") {
Expand Down Expand Up @@ -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))
Expand All @@ -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,
Expand All @@ -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 ) )
Expand Down
17 changes: 13 additions & 4 deletions R/deltasdrm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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") {
Expand Down
16 changes: 12 additions & 4 deletions R/deltasdrmb.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
}
Expand Down
17 changes: 13 additions & 4 deletions R/deltasdrmm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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") {
Expand Down

0 comments on commit bd629c1

Please sign in to comment.