dynamicALMP_quarter<- function(s,nPeriods,Ydur,S,X){
# dynamic ALMP estimator, reference here:
# VikstrĂ¶m, Johan. "Dynamic treatment assignment and evaluation of active labor market policies." Labour Economics 49 (2017): 42-54.
# equation (2) page 45 (without right-censoring)
# the object of interest in ATET_t(s) which is
# [average effect of treatment at s on the probability of surviving to time point t
# compared with survival throughout the same interval if never treated]
# two main identifying assumptions are:
# (A1) p44: sequential unconfoundedness
# (A2) p44: non-anticipation assumption
# s - (small s), time of program
# nPeriods - number of periods to look ahead
# Ydur - outcomes N x 1 duration of an unemployment spell
# X - covariate N x k matrix (here we assume that X do not change in time)
# change this into a matrix
# S - time between entry into unemployment and the start of the program N x 1 vector
# N - number of observations (individuals)
N <- dim(X)[1]
# T - number of time periods
#T <- 50 #4 years
T <- nPeriods+s #4 years ##CHANGE HERE
#tic
start_time <- Sys.time()
#unemployment duration
#Ydur <- apply(Y==0,1,sum)
#ps is N x T matrix
# the (i,t)-component is hap_p_t for person i
# that is Pr(S=s|X_is-, S>=s, barY_s-1 = 0)
# [estimated probability of obtaining treatment in period s
# given survival until time period s and covariates X]
# [NOTE the typo in the published article! Pr(S=t|X_is-, S>=t, barY_t-1 = 0), see WP]
ps <- matrix(0,N,T)
#ww is N x T matrix
# the (i,t)-component is hat_w_sk from equation (3) page 45
# it is the weight we need to use in order to get the correct
# counterfactual person for our comparisons
ww <- matrix(0,N,T)
#denominator of second term in eq (3)
prodD <- matrix(1,N,1)
#This needs to be run only once
# for given fixed s
# if s==2 it means a person was treated
# 2 months after he/she became unemployed
for (k in 1:T){
index <- which( as.logical( (S>=k)*(Ydur>=k)) )
#TO FIX: careful this index may be numeric(0)
SS <- (S==k)*1
# note that in application, X will be a matrix
ps[index,k] <- predict(glm(SS[index] ~ X[index,], family = binomial(link = "logit")) ,type="response")
#ps[index,k] <- max(0.01,ps[index,k])
for (iI in index){
ps[iI,k] <- min(0.99,ps[iI,k])
}
}
index <- which( as.logical((S>=s)*(Ydur>=s) ) )
print(length(index))
ww[index,s] <- (ps[index,s]/(1-ps[index,s]))
for (k in (s+1):(nPeriods+s)){
#product in the denominator
prodD <- prodD*(1-ps[,k])
ww[index,k] <- (ps[index,s]/(1-ps[index,s]))*(1/prodD[index])
}
#preallocate vector
ATETs <- matrix(0,T,1)
ATETs_unadjusted <- matrix(0,T,1)
#look nPeriods periods ahead
for (t in (s+1):(nPeriods+s)){
prod1 <- 1
prod0 <- 1
prod0_unadjusted <- 1
for (k in s:t){
#num1 <- (sum((Ydur>=k-1)*(S==s)*(Y[,k]==1)))
num1 <- (sum((Ydur>=k)*(S==s)*(Ydur==k)))
den1 <- (sum((Ydur>=k)*(S==s)))
prod1 <- prod1*(1-num1/den1)
#num0 <- (sum((Ydur>=k-1)*(S>k)*(Y[,k]==1)*ww[,k]))
num0 <- (sum((Ydur>=k)*(S>k)*(Ydur==k)*ww[,k]))
den0 <- (sum((Ydur>=k)*(S>k)*ww[,k]))
prod0 <- prod0*(1-num0/den0)
#unadjusted with weights
num0_unadjusted <- (sum((Ydur>=k)*(S>k)*(Ydur==k)))
den0_unadjusted <- (sum((Ydur>=k)*(S>k)))
prod0_unadjusted <- prod0_unadjusted*(1-num0_unadjusted/den0_unadjusted)
}
#print(prod1)
#print(prod0)
ATETs[t] <- prod1 - prod0
ATETs_unadjusted[t] <- prod1 - prod0_unadjusted
}
#plot it
times <- (s+1):(nPeriods+s)
ATETsg <- ATETs[times]
ATETsg_unadjusted <- ATETs_unadjusted[times]
timesg <- times-s
df <- data.frame(timesg,ATETsg)
gg <- ggplot(df)+
geom_line(aes(x=timesg, y=-ATETsg))+
labs(title=paste("ATET after ",toString(s)," months of unempl."), x="Months since start of the program", y="Effect on fraction with employment")+
scale_y_continuous(labels = percent)+
scale_x_continuous(breaks = c(1,2,3,6,12,18,24,30))
#toc
end_time <- Sys.time()
timeElapsed <- end_time - start_time
return(list("figure"=gg,"timeElapsed" = timeElapsed,"ATET" = ATETsg,
"ATET_un"= ATETsg_unadjusted, "PS" = ps))
}