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 RENT 0
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) #80% Sample ratio
# Create Null Data Frame
dfs<-data.frame("Year"=2001:2015)
Ballance.RENT<-Ballance.RENT%>%
mutate(lg=log(RENT_0))%>%
group_by(NS)%>%
filter(!any(is.na(RENT_0)))%>%data.frame()
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
dfs%>%
select(Year,average,sd,up,lo)%>%
head()
## Year average sd up lo
## 1 2001 -0.069269686 0.009037496 -0.051556520 -0.086982853
## 2 2002 0.044433867 0.003284384 0.050871141 0.037996592
## 3 2003 0.012232292 0.004044023 0.020158431 0.004306154
## 4 2004 -0.007147967 0.005277076 0.003194911 -0.017490845
## 5 2005 -0.002429225 0.004421067 0.006235907 -0.011094358
## 6 2006 0.309167321 0.007122833 0.323127817 0.295206824
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)
png("RENT0 SDID EVENT STUDY.png",width = 10,height = 5.6,units = "in",res = 400)
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)
dev.off()
## png
## 2
RENT 1
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) #80% Sample ratio
# Create Null Data Frame
dfs<-data.frame("Year"=2001:2015)
Ballance.RENT<-Ballance.RENT%>%
mutate(lg=log(RENT_1))%>%
group_by(NS)%>%
filter(!any(is.na(RENT_1)))%>%data.frame()
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
dfs%>%
select(Year,average,sd,up,lo)%>%
head()
## Year average sd up lo
## 1 2001 -0.076875650 0.009769244 -0.057728283 -0.09602302
## 2 2002 0.024104807 0.006055283 0.035972944 0.01223667
## 3 2003 -0.010683247 0.004331579 -0.002193509 -0.01917299
## 4 2004 -0.028104566 0.005962213 -0.016418843 -0.03979029
## 5 2005 -0.003694921 0.008754240 0.013463074 -0.02085292
## 6 2006 0.307560121 0.008235553 0.323701508 0.29141873
ggplot(dfs, aes(x = Year, y = average)) +
geom_point(shape=23, fill="blue", size=1.5)+
geom_line() + theme_bw()+ylab("RENT 1")+
geom_ribbon(aes(ymin = lo, ymax = up), alpha = 0.2)
png("RENT1 SDID EVENT STUDY.png",width = 10,height = 5.6,units = "in",res = 400)
ggplot(dfs, aes(x = Year, y = average)) +
geom_point(shape=23, fill="blue", size=1.5)+
geom_line() + theme_bw()+ylab("RENT 1")+
geom_ribbon(aes(ymin = lo, ymax = up), alpha = 0.2)
dev.off()
## png
## 2
RENT 2
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) #80% Sample ratio
# Create Null Data Frame
dfs<-data.frame("Year"=2001:2015)
Ballance.RENT<-Ballance.RENT%>%
mutate(lg=log(RENT_2))%>%
group_by(NS)%>%
filter(!any(is.na(RENT_2)))%>%data.frame()
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
dfs%>%
select(Year,average,sd,up,lo)%>%
head()
## Year average sd up lo
## 1 2001 -0.057859783 0.014816683 -0.028819619 -0.08689995
## 2 2002 0.025308792 0.004396113 0.033925016 0.01669257
## 3 2003 -0.002939130 0.006745612 0.010282027 -0.01616029
## 4 2004 -0.020352748 0.009375013 -0.001978061 -0.03872744
## 5 2005 -0.004809596 0.009316724 0.013450847 -0.02307004
## 6 2006 0.297595061 0.009266355 0.315756784 0.27943334
ggplot(dfs, aes(x = Year, y = average)) +
geom_point(shape=23, fill="blue", size=1.5)+
geom_line() + theme_bw()+ylab("RENT 2")+
geom_ribbon(aes(ymin = lo, ymax = up), alpha = 0.2)
png("RENT2 SDID EVENT STUDY.png",width = 10,height = 5.6,units = "in",res = 400)
ggplot(dfs, aes(x = Year, y = average)) +
geom_point(shape=23, fill="blue", size=1.5)+
geom_line() + theme_bw()+ylab("RENT 2")+
geom_ribbon(aes(ymin = lo, ymax = up), alpha = 0.2)
dev.off()
## png
## 2
RENT 3
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) #80% Sample ratio
# Create Null Data Frame
dfs<-data.frame("Year"=2001:2015)
Ballance.RENT<-Ballance.RENT%>%
mutate(lg=log(RENT_3))%>%
group_by(NS)%>%
filter(!any(is.na(RENT_3)))%>%data.frame()
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
dfs%>%
select(Year,average,sd,up,lo)%>%
head()
## Year average sd up lo
## 1 2001 -0.05202237 0.011194270 -0.0300820036 -0.073962735
## 2 2002 0.03282861 0.007389434 0.0473116315 0.018345583
## 3 2003 0.00762068 0.004698323 0.0168292244 -0.001587864
## 4 2004 -0.01150850 0.004271266 -0.0031369716 -0.019880029
## 5 2005 -0.01042528 0.005677094 0.0007016244 -0.021552177
## 6 2006 0.29554533 0.011293598 0.3176803703 0.273410281
ggplot(dfs, aes(x = Year, y = average)) +
geom_point(shape=23, fill="blue", size=1.5)+
geom_line() + theme_bw()+ylab("RENT 3")+
geom_ribbon(aes(ymin = lo, ymax = up), alpha = 0.2)
png("RENT3 SDID EVENT STUDY.png",width = 10,height = 5.6,units = "in",res = 400)
ggplot(dfs, aes(x = Year, y = average)) +
geom_point(shape=23, fill="blue", size=1.5)+
geom_line() + theme_bw()+ylab("RENT 3")+
geom_ribbon(aes(ymin = lo, ymax = up), alpha = 0.2)
dev.off()
## png
## 2
RENT 4
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) #80% Sample ratio
# Create Null Data Frame
dfs<-data.frame("Year"=2001:2015)
Ballance.RENT<-Ballance.RENT%>%
mutate(lg=log(RENT_4))%>%
group_by(NS)%>%
filter(!any(is.na(RENT_4)))%>%data.frame()
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
dfs%>%
select(Year,average,sd,up,lo)%>%
head()
## Year average sd up lo
## 1 2001 -0.045912908 0.012924172 -0.020581997 -0.07124382
## 2 2002 0.056821207 0.007572320 0.071662683 0.04197973
## 3 2003 0.040003293 0.004694016 0.049203395 0.03080319
## 4 2004 0.023558451 0.004488407 0.032355568 0.01476133
## 5 2005 -0.006294281 0.005501365 0.004488196 -0.01707676
## 6 2006 0.288854976 0.011087768 0.310586602 0.26712335
ggplot(dfs, aes(x = Year, y = average)) +
geom_point(shape=23, fill="blue", size=1.5)+
geom_line() + theme_bw()+ylab("RENT 4")+
geom_ribbon(aes(ymin = lo, ymax = up), alpha = 0.2)
png("RENT4 SDID EVENT STUDY.png",width = 10,height = 5.6,units = "in",res = 400)
ggplot(dfs, aes(x = Year, y = average)) +
geom_point(shape=23, fill="blue", size=1.5)+
geom_line() + theme_bw()+ylab("RENT 4")+
geom_ribbon(aes(ymin = lo, ymax = up), alpha = 0.2)
dev.off()
## png
## 2