Import
library(readxl)
library(ggplot2)
library(readr)
library(dplyr)
library(stringr)
library(plm)
library(lmtest)
library(synthdid)
RENT.merge<-read_xlsx("MSA-RENT merge All.xlsx")
RENT.merge$treat<-ifelse(RENT.merge$state=="LA",1,0)
colnames(RENT.merge)
## [1] "FirstNAME" "AREANAME" "RENT_0"
## [4] "RENT_1" "RENT_2" "RENT_3"
## [7] "RENT_4" "YEAR" "NAME.x"
## [10] "state" "Unemployee.Rate" "Personal.Income"
## [13] "Resident.Population" "nyear" "MSA"
## [16] "MSA.Code" "Price" "Change"
## [19] "NAME.y" "treat"
# Remove NA
RENT.merge<-RENT.merge%>%
group_by(MSA)%>%
filter(!any(is.na(Unemployee.Rate)))
dloop<-data.frame(YEAR=unique(RENT.merge$YEAR))
## YEAR DUMMY
for (i in 1:nrow(dloop)) {
RENT.merge[paste0("d.",
dloop$YEAR[i])] <- as.numeric(RENT.merge$YEAR == unique(RENT.merge$YEAR)[i])
}
#KEEP NEW ORLEAN
RENT.merge1<-RENT.merge
RENT.merge1$NS<-paste0(RENT.merge1$FirstNAME,RENT.merge1$state)
RENT.NEWLA<-RENT.merge1%>%
filter(NS=="newLA")
RENTUS<-RENT.merge1%>%
filter(state!="LA")
RENT.merge1<-RENT.merge1%>%
filter(NS!="newLA")
## WITH NEW ORLEAN ONLY
RENT.merge2<-rbind(RENT.NEWLA,RENTUS)
RENT.merge2$YEAR<-as.numeric(RENT.merge2$YEAR)
RENT 0
RENT.merge2$s.treat<-ifelse(RENT.merge2$state=="LA" &
RENT.merge2$YEAR>2005,1,0)
##BALLANCE BY NS
unq<-data.frame(NS=rep(unique(RENT.merge2$NS),15))
unq<-arrange(unq,NS)
unq$YEAR<-rep(2001:2015,length(unique(RENT.merge2$NS)))
Ballance.RENT<-left_join(unq,RENT.merge2,by=c("YEAR","NS"))
Ballance.RENT<-Ballance.RENT%>%
group_by(NS)%>%
filter(!any(is.na(RENT_0)))
Ballance.RENT$lg<-log(Ballance.RENT$RENT_0)
#Ballance.RENT$NS<-factor(Ballance.RENT$NS)
Ballance.RENT<-data.frame(Ballance.RENT)
panelp<-panel.matrices(Ballance.RENT,unit = "NS",time = "YEAR",
outcome = "lg",treatment = "s.treat")
## SDID Estimate
r0<-synthdid_estimate(Y=panelp$Y,N0 = panelp$N0,T0 = panelp$T0)
plot(r0,overlay=1)+theme(legend.position = "bottom",
legend.direction = "horizontal",
legend.box.background = element_blank())+
ylab("RENT 0")+xlab("Year")
## Sc estimate
sc0<-sc_estimate(Y=panelp$Y,N0 = panelp$N0,T0 = panelp$T0)
plot(sc0)+theme(legend.position = "bottom",
legend.direction = "horizontal",
legend.box.background = element_blank())+
ylab("RENT 0")+xlab("Year")
Dot
weights = attr(r0, "weights")
Y = panelp$Y
T0 = panelp$T0; T1 = ncol(Y) - T0
N0 = panelp$N0; N1 = nrow(Y) - N0
#omega.synth = c(weights$omega, rep(0, N1))
omega.synth = c(weights$omega, rep(0, N1))
omega.target = c(rep(0, N0), rep(1 / N1, N1))
lambda.synth = c(weights$lambda, rep(0, T1)) #lambda sdid
lambda.target = c(rep(0, T0), rep(1 / T1, T1))
#over = attr(r0, 'overlay')
over=1
intercept.offset = over * c((omega.target - omega.synth) %*% Y %*% lambda.synth)
obs.trajectory = as.numeric(omega.target %*% Y) #Ytr Treat
#syn.trajectory.i = as.numeric(omega.synth %*% Y) + intercept.offset
syn.trajectory=as.numeric(omega.synth %*% Y) #Yco
# [Yt(tr)-Yt(co)]-[Yb(tr)-Yb(co)] # Yt is ploted and Yb is Baseline
##USED after treatment data
dfw<-data.frame("lambda"=lambda.synth[5:15],
"ytr"=obs.trajectory[5:15],
"yco"=syn.trajectory[5:15])
dfw$aux=dfw$lambda*(dfw$ytr-dfw$yco)
head(dfw)
## lambda ytr yco aux
## 1 0.8310792 6.300786 6.316231 -0.01283596
## 2 0.0000000 6.645091 6.353434 0.00000000
## 3 0.0000000 6.684612 6.384448 0.00000000
## 4 0.0000000 6.694562 6.416096 0.00000000
## 5 0.0000000 6.656727 6.457127 0.00000000
## 6 0.0000000 6.690842 6.483999 0.00000000
meanpre_0<-dfw$aux[1] ## Export Average pre treatment
differ<-data.frame("year"=2001:2015,"d"=obs.trajectory- syn.trajectory) #Ytr-Yco
differ$d<-differ$d- meanpre_0 #Equation 8
ggplot(differ,aes(x=year,y=d))+geom_point()+
geom_line()+theme_bw()+
ylab("RENT 0")
Loop
ltm<-Ballance.RENT%>%
filter(NS!="newLA")
ltm<-unique(ltm$NS)
lno<-Ballance.RENT%>%
filter(NS=="newLA")
lno<-unique(lno$NS)
#Number of sample
ns=round(length(ltm)*.8)
# Create Null Data Frame
dfs<-data.frame("Year"=2001:2015)
for (i in 1:100){ #Loop 100 times
sa.c<-sample(x = ltm,size = ns)
sa.c<-c(sa.c,lno)
sa.Ballance<-Ballance.RENT%>%
filter(NS %in% sa.c)
sa.panel<-panel.matrices(sa.Ballance,
unit = "NS",
time = "YEAR",
outcome = "lg",
treatment = "s.treat")
sa.sdid<-synthdid_estimate(Y = sa.panel$Y,
N0 = sa.panel$N0,
T0 = sa.panel$T0)
weights = attr(sa.sdid, "weights")
Y = sa.panel$Y
T0 = sa.panel$T0; T1 = ncol(Y) - T0
N0 = sa.panel$N0; N1 = nrow(Y) - N0
#omega.synth = c(weights$omega, rep(0, N1))
omega.synth = c(weights$omega, rep(0, N1))
omega.target = c(rep(0, N0), rep(1 / N1, N1))
lambda.synth = c(weights$lambda, rep(0, T1)) #lambda sdid
lambda.target = c(rep(0, T0), rep(1 / T1, T1))
#over = attr(sa.sdid, 'overlay')
over=1
intercept.offset = over * c((omega.target - omega.synth) %*% Y %*% lambda.synth)
obs.trajectory = as.numeric(omega.target %*% Y) #Ytr Treat
#syn.trajectory.i = as.numeric(omega.synth %*% Y) + intercept.offset
syn.trajectory=as.numeric(omega.synth %*% Y) #Yco
# [Yt(tr)-Yt(co)]-[Yb(tr)-Yb(co)] # Yt is ploted and Yb is Baseline
##USED after treatment data
dfw<-data.frame("lambda"=lambda.synth[5:15],
"ytr"=obs.trajectory[5:15],
"yco"=syn.trajectory[5:15])
dfw$aux=dfw$lambda*(dfw$ytr-dfw$yco)
#head(dfw)
meanpre_0=dfw$aux[1] ## Export Average pre treatment
dfs[,i+1]=(obs.trajectory- syn.trajectory)- meanpre_0 #Save Looping Output
colnames(dfs)[i+1]=paste0("d",i) #Rename Looping output
}
dfs<-dfs%>%
mutate(average=apply(dfs[,-1],1,FUN = mean),
sd=apply(dfs[,-1],1,FUN = sd))%>%
mutate(up=average + sd*qnorm(0.975), lo=average - sd*qnorm(0.975)) #95% Confident Interval
ggplot(dfs, aes(x = Year, y = average)) +
geom_point(shape=23, fill="blue", size=1.5)+
geom_line() + theme_bw()+ylab("RENT 0")+
geom_ribbon(aes(ymin = lo, ymax = up), alpha = 0.2)
