#=====================================================================#
# This is the third R code used in the empirical study of the paper:
# "A Distributional Synthetic Control Method for Policy Evaluation"
#=====================================================================#

#--------------------------------------------------------------------------------------------------------------------------------------
#=========================#
# Part I. Define Function # 
#=========================#
#########################################################
# Function: Transfer quantiles to mean & Gini coefficient
#########################################################
MGQ=function(Mtx_Q,tau){  
T=nrow(Mtx_Q)
mean=matrix(0,T,1)                
Gini=matrix(0,T,1)                
Lorenz=matrix(0,T,length(tau))                
s=1
while(s<=T){
 mean[s]=sum(Mtx_Q[s,])*(tau[2]-tau[1])
Lorenz[s,]=t(cumsum(Mtx_Q[s,])*(tau[2]-tau[1]))/mean[s]
 Gini[s]=1-2*(sum(Lorenz[s,])*(tau[2]-tau[1]))
s=s+1}
list(res1=mean,res2=Lorenz,res3=Gini)}

##########################################
# Function: X-axis for time series plots #
##########################################
X_axis=function(Period,q){  
start_Y=substring(Period[1],1,4)     
start_Q=substring(Period[1],5,6)
  end_Y=substring(Period[length(Period)],1,4)
  end_Q=substring(Period[length(Period)],5,6)
if(start_Q=="01"){s_D="/3/30"}
if(start_Q=="02"){s_D="/6/30"}
if(start_Q=="03"){s_D="/9/30"}
if(start_Q=="04"){s_D="/12/31"}
  if(end_Q=="01"){e_D="/3/30"}
  if(end_Q=="02"){e_D="/6/30"}
  if(end_Q=="03"){e_D="/9/30"}
  if(end_Q=="04"){e_D="/12/31"}
start_D=paste0(start_Y,s_D)
  end_D=paste0(end_Y,e_D)
Ptime=seq(as.Date(start_D),as.Date(end_D),by = "quarter")
Paxis=Ptime[seq(1,length(Ptime),q)]
lPaxis=as.numeric(substring(Paxis,1,4))   # X-axis label 
list(res1=Ptime,res2=Paxis,res3=lPaxis)}


#--------------------------------------------------------------------------------------------------------------------------------------
#============================#
# Part II. Empirical Results # 
#============================#

#*****************************************************#
# 1. Load Data & Define Variables                     #
# (See "1.LoadData_and_DefineEvents.R" for this part) #
#*****************************************************#
dirpath="e:/research/SyntheticControl/Empirical/"
Data=read.csv(paste0(dirpath,"Data/","county data.csv"))

# Variable definitions:
       time=Data$time                 # YYYYQQ (Sampling period: 1990Q1-2006Q2) 
      state=Data$statename            # state names
      cfips=Data$countyreal           # county fips, https://en.wikipedia.org/wiki/FIPS_county_code 
     county=Data$countyname           # county names
       cpop=Data$countypop2000        # county population, 2000
      cpopd=Data$cntypopdens          # county population density, 2000
      carea=Data$cntyarea             # county land area
     mw_fed=Data$federalmin           # federal min. wage 
     mw_sta=Data$stminwage            # state min. wage
         mw=Data$minwage              # max(mw_fed,mw_sta) 
    emp_tot=Data$empTOT               # total employment
   emp_rest=Data$emp_rest_both        # restaurant employment (NAICS 7221+7222)
  wage_rest=Data$AWW_rest_both        # restaurant average weekly wage
    emp_722=Data$emp722               # food services and drinking place employment
   wage_722=Data$AWW722               # food services and drinking place average weekly wage
    emp_ret=Data$empRETAIL            # retail employment
   wage_ret=Data$AWWRETAIL            # retail average weekly wage                                           
    emp_mfg=Data$empMFG               # manufacturing employment
   wage_mfg=Data$AWWMFG               # manufacturing average weekly wage

# State indices & names:
NT=length(state)                      # NT: full sample size 
sindx=matrix(0,NT,1)                  # state indices (1,2,3,...)

# Precluding Alaska & Hawaii 
sindx[which(state=="Alabama")]=1;               
sindx[which(state=="Arizona")]=2;               
sindx[which(state=="Arkansas")]=3;              
sindx[which(state=="California")]=4;            
sindx[which(state=="Colorado")]=5;              
sindx[which(state=="Connecticut")]=6;           
sindx[which(state=="Delaware")]=7;              
sindx[which(state=="District of Columbia")]=8;  
sindx[which(state=="Florida")]=9;               
sindx[which(state=="Georgia")]=10;              
sindx[which(state=="Idaho")]=11;                                                                                                              
sindx[which(state=="Illinois")]=12;                                                                                                              
sindx[which(state=="Indiana")]=13;                                                                                                              
sindx[which(state=="Iowa")]=14;                                                                                                              
sindx[which(state=="Kansas")]=15;                                                                                                              
sindx[which(state=="Kentucky")]=16;                                                                                                              
sindx[which(state=="Louisiana")]=17;                                                                                                              
sindx[which(state=="Maine")]=18;                                                                                                              
sindx[which(state=="Maryland")]=19;                                                                                                              
sindx[which(state=="Massachusetts")]=20;                                                                     
sindx[which(state=="Michigan")]=21;                                                                     
sindx[which(state=="Minnesota")]=22;                                                                     
sindx[which(state=="Mississippi")]=23;                                                                     
sindx[which(state=="Missouri")]=24;                                                                     
sindx[which(state=="Montana")]=25;                                                                     
sindx[which(state=="Nebraska")]=26;                                                                     
sindx[which(state=="Nevada")]=27;                                                                     
sindx[which(state=="New Hampshire")]=28;                                                                     
sindx[which(state=="New Jersey")]=29;                                                                     
sindx[which(state=="New Mexico")]=30;                           
sindx[which(state=="New York")]=31;                           
sindx[which(state=="North Carolina")]=32;                           
sindx[which(state=="North Dakota")]=33;                           
sindx[which(state=="Ohio")]=34;                           
sindx[which(state=="Oklahoma")]=35;                           
sindx[which(state=="Oregon")]=36;                           
sindx[which(state=="Pennsylvania")]=37;                           
sindx[which(state=="Rhode Island")]=38;                           
sindx[which(state=="South Carolina")]=39;                           
sindx[which(state=="South Dakota")]=40;         
sindx[which(state=="Tennessee")]=41;            
sindx[which(state=="Texas")]=42;                
sindx[which(state=="Utah")]=43;                 
sindx[which(state=="Vermont")]=44;              
sindx[which(state=="Virginia")]=45;             
sindx[which(state=="Washington")]=46;           
sindx[which(state=="West Virginia")]=47;        
sindx[which(state=="Wisconsin")]=48;            
sindx[which(state=="Wyoming")]=49;              

minci=min(cfips[which(sindx==1)]) 
Time=time[which(cfips==minci)]    # time index (1990Q1-2006Q2: 199001,199002,...,200602)
T=length(Time)                    # number of quarters 
N=49                              # number of states (precluding Alaska & Hawaii)
Y=cbind(wage_rest,emp_rest) 
y1=Y[,1]        # wages
y2=Y[,2]/100    # (rescaled) employments

Sindx=matrix(0,N,1)      # state indices
Sname=matrix(0,N,1)      # state names 
Scode=matrix(0,N,1)      # state codes
 
i=1
while(i<=N){
s_indx=sindx[which(sindx==i)]
s_name=state[which(sindx==i)]
Sindx[i]=s_indx[i]
Sname[i]=as.character(s_name[i])
if(Sname[i]=="Alabama"){Scode[i]="AL"}                                  
if(Sname[i]=="Arizona"){Scode[i]="AZ"}                                  
if(Sname[i]=="Arkansas"){Scode[i]="AR"}                                 
if(Sname[i]=="California"){Scode[i]="CA"}                               
if(Sname[i]=="Colorado"){Scode[i]="CO"}                                 
if(Sname[i]=="Connecticut"){Scode[i]="CT"}                              
if(Sname[i]=="Delaware"){Scode[i]="DE"}                                 
if(Sname[i]=="District of Columbia"){Scode[i]="DC"}                     
if(Sname[i]=="Florida"){Scode[i]="FL"}                                  
if(Sname[i]=="Georgia"){Scode[i]="GA"}                                  
if(Sname[i]=="Idaho"){Scode[i]="ID"}                                    
if(Sname[i]=="Illinois"){Scode[i]="IL"}                                 
if(Sname[i]=="Indiana"){Scode[i]="IN"}                                  
if(Sname[i]=="Iowa"){Scode[i]="IA"}                                     
if(Sname[i]=="Kansas"){Scode[i]="KS"}                                   
if(Sname[i]=="Kentucky"){Scode[i]="KY"}                                 
if(Sname[i]=="Louisiana"){Scode[i]="LA"}                                
if(Sname[i]=="Maine"){Scode[i]="ME"}                                    
if(Sname[i]=="Maryland"){Scode[i]="MD"}                                 
if(Sname[i]=="Massachusetts"){Scode[i]="MA"}                            
if(Sname[i]=="Michigan"){Scode[i]="MI"}                                 
if(Sname[i]=="Minnesota"){Scode[i]="MN"}                                
if(Sname[i]=="Mississippi"){Scode[i]="MS"}                              
if(Sname[i]=="Missouri"){Scode[i]="MO"}                                 
if(Sname[i]=="Montana"){Scode[i]="MT"}                                  
if(Sname[i]=="Nebraska"){Scode[i]="NE"}                                 
if(Sname[i]=="Nevada"){Scode[i]="NV"}                                   
if(Sname[i]=="New Hampshire"){Scode[i]="NH"}                            
if(Sname[i]=="New Jersey"){Scode[i]="NJ"}                               
if(Sname[i]=="New Mexico"){Scode[i]="NM"}                               
if(Sname[i]=="New York"){Scode[i]="NY"}                                 
if(Sname[i]=="North Carolina"){Scode[i]="NC"}                           
if(Sname[i]=="North Dakota"){Scode[i]="ND"}                             
if(Sname[i]=="Ohio"){Scode[i]="OH"}                                     
if(Sname[i]=="Oklahoma"){Scode[i]="OK"}                                 
if(Sname[i]=="Oregon"){Scode[i]="OR"}                                   
if(Sname[i]=="Pennsylvania"){Scode[i]="PA"}                             
if(Sname[i]=="Rhode Island"){Scode[i]="RI"}                             
if(Sname[i]=="South Carolina"){Scode[i]="SC"}                           
if(Sname[i]=="South Dakota"){Scode[i]="SD"}                             
if(Sname[i]=="Tennessee"){Scode[i]="TN"}                                
if(Sname[i]=="Texas"){Scode[i]="TX"}                                    
if(Sname[i]=="Utah"){Scode[i]="UT"}                                     
if(Sname[i]=="Vermont"){Scode[i]="VT"}                                  
if(Sname[i]=="Virginia"){Scode[i]="VA"}                                 
if(Sname[i]=="Washington"){Scode[i]="WA"}                               
if(Sname[i]=="West Virginia"){Scode[i]="WV"}                            
if(Sname[i]=="Wisconsin"){Scode[i]="WI"}                                
if(Sname[i]=="Wyoming"){Scode[i]="WY"}                                  
minci=min(cfips[which(sindx==i)])
nmw_sta=mw_sta
nmw_sta[is.na(nmw_sta)]=0
i=i+1}

BSname=Sname
BSname[8]="D. of Columbia"
BSname[32]="N. Carolina"
BSname[33]="N. Dakota"
BSname[39]="S. Carolina"
BSname[40]="S. Dakota"
BSname[47]="W. Virginia"

# Set quantile probability indices
tau=seq(0.01,0.99,0.01)     

#*********************************************************************#
# 3. Four Case Studies                                                #
# (See "2.Distributional_Synthetic_Control_Analysis.R" for this part) #
#*********************************************************************#
case=1
while(case<=4){
if(case==1){trt_name="California"; Pre_Smw=4.25; Size_Smwch=0.75; Pre_Emw=4.75; Size_Emwch=0.25; PreStart=199001; PreEnd=199604; PostStart=199701; PostEnd=199702; PreDura=28; PostDura=2 } 
if(case==2){trt_name="Oregon";     Pre_Smw=4.75; Size_Smwch=0.75; Pre_Emw=4.75; Size_Emwch=0.75; PreStart=199101; PreEnd=199604; PostStart=199701; PostEnd=199704; PreDura=24; PostDura=4 } 
if(case==3){trt_name="Washington"; Pre_Smw=4.9 ; Size_Smwch=0.8 ; Pre_Emw=5.15; Size_Emwch=0.55; PreStart=199401; PreEnd=199804; PostStart=199901; PostEnd=199904; PreDura=20; PostDura=4 } 
if(case==4){trt_name="Wisconsin";  Pre_Smw=3.65; Size_Smwch=2.05; Pre_Emw=5.15; Size_Emwch=0.55; PreStart=199001; PreEnd=200501; PostStart=200502; PostEnd=200602; PreDura=61; PostDura=4 } 
cat("case=",case,"Treated=",trt_name,"\n")
# Define time indices
a=which(Time==PreStart)   # time index of the start of a pre-intervention period  
b=which(Time==PreEnd)     # time index of the end of a pre-intervention period
c=which(Time==PostStart)  # time index of the start of a post-intervention period (policy-change point) 
d=which(Time==PostEnd)    # time index of the end of a post-intervention period    
Pre=Time[a:b]             # pre-intervention period
Post=Time[c:d]            # post-intervention period
Period=c(Pre,Post)        # case-study period
res=X_axis(Period,4)
Xtime=res[[1]]
Xaxis=res[[2]]
x_lab=res[[3]]

###################################
# I. Results for the Wage Effects # 
###################################

#--------------------------------------------------------------#
# 1. Load empirical results (generated from the second R code) #
#--------------------------------------------------------------#
pcu_indx0 = read.csv(paste0(dirpath, "Data/Results/", "pcu_indx_", case, ".csv"))
TRT0=read.csv(paste0(dirpath,"Data/Results/","TRT_",case,".csv"))
SYN0=read.csv(paste0(dirpath,"Data/Results/","SYN_",case,".csv"))
SYN00=read.csv(paste0(dirpath,"Data/Results/","SYN0_",case,".csv"))
wSYN0=read.csv(paste0(dirpath,"Data/Results/","wSYN_",case,".csv"))
wSYN00=read.csv(paste0(dirpath,"Data/Results/","wSYN0_",case,".csv"))
trtPre_fit0=read.csv(paste0(dirpath,"Data/Results/","trtPre_fit_",case,".csv"))
pcuPre_fit0=read.csv(paste0(dirpath,"Data/Results/","pcuPre_fit_",case,".csv"))
pcu_indx=as.matrix(pcu_indx0[,-1])  # Remove the 1st column (ordering number)
TRT=as.matrix(TRT0[,-1])
SYN=as.matrix(SYN0[,-1])
SYN0=as.matrix(SYN00[,-1])
wSYN=as.matrix(wSYN0[,-1])
wSYN0=as.matrix(wSYN00[,-1])
trtPre_fit=as.matrix(trtPre_fit0[,-1])
pcuPre_fit=as.matrix(pcuPre_fit0[,-1])
J=length(pcu_indx)
j=1
while(j<=J){
var=paste0("plb_pcu_indx_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plb_pcu_indx_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_TRT_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plb_TRT_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_SYN_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plb_SYN_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_SYN0_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plb_SYN0_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_wSYN_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plb_wSYN_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_wSYN0_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plb_wSYN0_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_trtPre_fit_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plb_trtPre_fit_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_pcuPre_fit_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plb_pcuPre_fit_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
j=j+1}

#--------------------#
# 2. Optimal weights #
#--------------------#
### Figure 2  (Panel A)
# mean weight
mt=paste0("Case ",case,":"," ",trt_name) 
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/MeanWeights",case,".pdf",sep=""))
optw=rowMeans(wSYN)
bar=barplot(as.vector(optw), main=mt, cex.main=1.5, ylim=c(0,1), col="grey", angle=45, 
border=TRUE,las=2,cex.names=1.2)
text(bar,as.vector(optw)+0.05,round(as.vector(optw),2),cex=1.2)
text(bar,BSname[pcu_indx],y = -0.01,cex = 1.1,srt = 60, xpd = TRUE, pos = 2 )
dev.off()}

### Figure A.12 (Panel A)
# quantile weights
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/QuantileWeights",case,".pdf",sep=""))
cr=rainbow(length(tau),s=1,v=1,start=0,end=0.9)
lg=tau
lcr=c(cr[10],cr[20],cr[30],cr[40],cr[50],cr[60],cr[70],cr[80],cr[90])
llg=c(lg[10],lg[20],lg[30],lg[40],lg[50],lg[60],lg[70],lg[80],lg[90])
optw=rowMeans(wSYN)
bar=barplot(as.vector(optw), main=mt, cex.main=1.5, ylim=c(0,1), col=NA, angle=45, 
border=FALSE,las=2,cex.names=1.2)
text(bar,BSname[pcu_indx],y = -0.01,cex = 1.1,srt = 60, xpd = TRUE, pos = 2 )
i=1
while(i<=length(tau)){
points(x=bar,wSYN[,i],col=cr[i],type="l")
i=i+1}
if(case==1){legend("topleft",legend=llg,lty=1,col=lcr,bty="n",lwd=1,cex=1.2,y.intersp=1)}
dev.off()}

#--------------------------------------------------------------------------------------#
# 3. Compare the economic features of the treated unit and the synthetic control units #
#--------------------------------------------------------------------------------------#
if(case==1){XData=read.csv(paste0(dirpath,"Data/","US_state_comparison1996.csv"))}
if(case==2){XData=read.csv(paste0(dirpath,"Data/","US_state_comparison1996.csv"))}
if(case==3){XData=read.csv(paste0(dirpath,"Data/","US_state_comparison1998.csv"))}
if(case==4){XData=read.csv(paste0(dirpath,"Data/","US_state_comparison2004.csv"))}
XData=as.matrix(XData[,-1])
X0_pop=XData[1,]       # Resident population
X0_white=XData[2,]     # White alone, percentage
X0_old=XData[3,]       # 65 years old and over, percentage
X0_land=XData[4,]      # Land area
X0_pden=XData[5,]      # Population per square mile (Preclude this from the matching condition because it is redundant)
X0_pov=XData[6,]       # Person below poverty level, percentage
X0_pinc=XData[7,]      # Personal income per capita
X0_hinc=XData[8,]      # Real median household income
X0_emp=XData[9,]       # Annual nonfarm employment
X0_nhu=XData[10,]      # New private housing units authorized by buliding permits
X0econ=as.matrix(rbind(X0_pop,X0_white,X0_old,X0_land,X0_pov,X0_pinc,X0_hinc,X0_emp,X0_nhu))

trt_indx=Sindx[which(Sname==trt_name)]
Xecon_trt=X0econ[,trt_indx]
Xecon_pcu=X0econ[,pcu_indx]
optw=as.matrix(rowMeans(wSYN))
Xecon_syn=Xecon_pcu%*%optw
optw01=as.matrix(wSYN[,10])
Xecon_syn01=Xecon_pcu%*%optw01
optw05=as.matrix(wSYN[,50])
Xecon_syn05=Xecon_pcu%*%optw05
optw09=as.matrix(wSYN[,90])
Xecon_syn09=Xecon_pcu%*%optw09
Distance=matrix(0,length(pcu_indx),1)
i=1
while(i<=nrow(Distance)){
Distance[i]=(mean((Xecon_trt-Xecon_pcu[,i])^2))^0.5
i=i+1}
singleindx=pcu_indx[which.min(Distance)]
singlename=Sname[singleindx]
Xecon_single=X0econ[,singleindx]

### Table 2 (Panel A)
print("wage effects")
print(singlename)
print(cbind(round(Xecon_trt,2),round(Xecon_syn,2),round(Xecon_syn01,2),round(Xecon_syn05,2),round(Xecon_syn09,2),round(Xecon_single,2)))


#---------------------------#
# 4. Performance evaluation #
#---------------------------#
# Mean sqaured (prediction) errors
preMSE_SYN=matrix(0,length(tau),1)
preMSE_SYN0=matrix(0,length(tau),1)
postMSE_SYN=matrix(0,length(tau),1)
postMSE_SYN0=matrix(0,length(tau),1)
q=1
while(q<=length(tau)){
tauq=tau[q]
preMSE_SYN[q]=mean((TRT[1:(b-a+1),q]-SYN[1:(b-a+1),q])^2)
preMSE_SYN0[q]=mean((TRT[1:(b-a+1),q]-SYN0[1:(b-a+1),q])^2)
postMSE_SYN[q]=mean((TRT[(c-a+1):(d-a+1),q]-SYN[(c-a+1):(d-a+1),q])^2)
postMSE_SYN0[q]=mean((TRT[(c-a+1):(d-a+1),q]-SYN0[(c-a+1):(d-a+1),q])^2)
q=q+1}

### Figure 3 (Panel A)
if(case==1){
mt=paste0("Case ",case,":"," ",trt_name) 
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/preMSE",case,".pdf",sep=""))
matplot(tau,cbind(preMSE_SYN,preMSE_SYN0),main=mt, cex.main=1.5, ylim=c(0,300),xlab=" ",ylab=" ",type="l",col=c(2,4),lty=c(1,2),lwd=c(3,3),xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
legend("topleft",legend=c("proposed","conventional"),col=c(2,4),lty=c(1,2),lwd=c(3,3),bty="n",cex=1.8,y.intersp=1)
dev.off()}}
if(case>1 & case<4){
mt=paste0("Case ",case,":"," ",trt_name) 
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/preMSE",case,".pdf",sep=""))
matplot(tau,cbind(preMSE_SYN,preMSE_SYN0),main=mt, cex.main=1.5, ylim=c(0,100),xlab=" ",ylab=" ",type="l",col=c(2,4),lty=c(1,2),lwd=c(3,3),xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
dev.off()}}
if(case==4){
mt=paste0("Case ",case,":"," ",trt_name) 
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/preMSE",case,".pdf",sep=""))
matplot(tau,cbind(preMSE_SYN,preMSE_SYN0),main=mt, cex.main=1.5, ylim=c(0,100),xlab=" ",ylab=" ",type="l",col=c(2,4),lty=c(1,2),lwd=c(3,3),xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
dev.off()}}

#-----------------------------------------#
# 5. Mean & quantile intervention effects #
#-----------------------------------------#
T=nrow(TRT)
res=MGQ(TRT,tau)
mean_trt=res[[1]]
res=MGQ(SYN,tau)
mean_syn=res[[1]]
res=MGQ(SYN0,tau)
mean_syn0=res[[1]]

### Figure 4
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/TimeSeriesEffect",case,".pdf",sep=""))
par(mfrow=c(1,4),mar=c(2,2.5,2.5,2))
ymin=min(cbind(TRT[,10],SYN[,10],SYN0[,10]))
ymax=max(cbind(TRT[,90],SYN[,90],SYN0[,90]))

if(case==1){mtm=expression("Case1: mean")}
if(case==2){mtm=expression("Case2: mean")}
if(case==3){mtm=expression("Case3: mean")}
if(case==4){mtm=expression("Case4: mean")}
matplot(Xtime,cbind(mean_trt,mean_syn,mean_syn0),main=mtm, cex.main=1.8,ylim=c(ymin,ymax),type="l",lty=c(1,2,4),col=c(2,4,1),lwd=c(3,3,3),pch=19,cex=0.5,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
if(case==1){legend("topleft",legend=c("actual","proposed","conventional"),lty=c(1,2,4),col=c(2,4,1),lwd=c(3,3,3),bty="n",cex=1.2,seg.len=0.5,y.intersp=1.3)}
abline(v=Xtime[(c-a+1)],lty=5,lwd=2,col="darkblue")
axis(1,cex.axis=1.5,at=Xaxis,labels=x_lab)
axis(2,cex.axis=1.5)

if(case==1){mt1=expression("Case1:" ~ tau ~ "=0.1")}
if(case==2){mt1=expression("Case2:" ~ tau ~ "=0.1")}
if(case==3){mt1=expression("Case3:" ~ tau ~ "=0.1")}
if(case==4){mt1=expression("Case4:" ~ tau ~ "=0.1")}
matplot(Xtime,cbind(TRT[,10],SYN[,10],SYN0[,10]),main=mt1, cex.main=1.8,ylim=c(ymin,ymax),type="l",lty=c(1,2,4),col=c(2,4,1),lwd=c(3,3,3),pch=19,cex=0.5,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
abline(v=Xtime[(c-a+1)],lty=5,lwd=2,col="darkblue")
axis(1,cex.axis=1.5,at=Xaxis,labels=x_lab)
axis(2,cex.axis=1.5)

if(case==1){mt5=expression("Case1:" ~ tau ~ "=0.5")}
if(case==2){mt5=expression("Case2:" ~ tau ~ "=0.5")}
if(case==3){mt5=expression("Case3:" ~ tau ~ "=0.5")}
if(case==4){mt5=expression("Case4:" ~ tau ~ "=0.5")}
matplot(Xtime,cbind(TRT[,50],SYN[,50],SYN0[,50]),main=mt5, cex.main=1.8,ylim=c(ymin,ymax),type="l",lty=c(1,2,4),col=c(2,4,1),lwd=c(3,3,3),pch=19,cex=0.5,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
abline(v=Xtime[(c-a+1)],lty=5,lwd=2,col="darkblue")
axis(1,cex.axis=1.5,at=Xaxis,labels=x_lab)
axis(2,cex.axis=1.5)

if(case==1){mt9=expression("Case1:" ~ tau ~ "=0.9")}
if(case==2){mt9=expression("Case2:" ~ tau ~ "=0.9")}
if(case==3){mt9=expression("Case3:" ~ tau ~ "=0.9")}
if(case==4){mt9=expression("Case4:" ~ tau ~ "=0.9")}
matplot(Xtime,cbind(TRT[,90],SYN[,90],SYN0[,90]),main=mt9, cex.main=1.8,ylim=c(ymin,ymax),type="l",lty=c(1,2,4),col=c(2,4,1),lwd=c(3,3,3),pch=19,cex=0.5,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
abline(v=Xtime[(c-a+1)],lty=5,lwd=2,col="darkblue")
axis(1,cex.axis=1.5,at=Xaxis,labels=x_lab)
axis(2,cex.axis=1.5)
par(mfrow=c(1,1))
dev.off()}

#-------------------------------#
# 6. Mean effect & placebo test #
#-------------------------------#
Mdiff=mean_trt-mean_syn
Rmean_SYN=mean(Mdiff[(c-a+1):(d-a+1)]^2)/mean(Mdiff[1:(b-a+1)]^2)
plb_Mdiff=matrix(0,T,J)
plb_Rmean_SYN=matrix(0,J,1)
j=1
while(j<=J){
plb_TRT=get(paste0("plb_TRT_",case,"_",j))
plb_SYN=get(paste0("plb_SYN_",case,"_",j))
res=MGQ(plb_TRT,tau)
mean_plbtrt=res[[1]]
res=MGQ(plb_SYN,tau)
mean_plbsyn=res[[1]]
plb_Mdiff[,j]=mean_plbtrt-mean_plbsyn
plb_Rmean_SYN[j]=mean(plb_Mdiff[(c-a+1):(d-a+1),j]^2)/mean(plb_Mdiff[1:(b-a+1),j]^2)
j=j+1}
Meffect=Mdiff[(c-a+1):(d-a+1)]
Delta_Mdiff=Meffect/mean_syn[(c-a+1):(d-a+1)]
Delta_Emw=Size_Emwch/Pre_Emw
Elast_Mdiff=Delta_Mdiff/Delta_Emw
Rmean=c(Rmean_SYN,plb_Rmean_SYN)
Rrank=rank(Rmean)
Pvalue=1-Rrank[1]/(J+1)

########################
### Table 2: mean effects (part A)
s=1
while(s<=length(Meffect)){
cat("case=",case,"mean effect on earnings","s=",s,"level=",round(Meffect[s],3),"elasticity=",round(Elast_Mdiff[s],3), "\n")
s=s+1}
cat("case=",case,"mean effect on earnings","p-value=",round(Pvalue,3),"\n")
# mean wages at To,To+1,...
cat("case=",case,"actual mean earnings","\n")
print(round(mean_trt[(b-a+1):(d-a+1)],3))


### Figure 6 (Panel A)
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/MeanPlacebo",case,".pdf",sep=""))
matplot(Xtime,plb_Mdiff,main=mt, cex.main=1.5,col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
lines(Xtime,Mdiff,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
abline(v=Xtime[(c-a+1)],lty=5,lwd=2,col="darkblue")
abline(h=0,lty=5,lwd=2,col="darkblue")
axis(1,cex.axis=1.5,at=Xaxis,labels=x_lab)
axis(2,cex.axis=1.5)
if(case==1){legend("topleft",legend=c("treated","placebo"),lty=c(1,1),col=c(2,gray(0.5)),lwd=c(3,1),bty="n",cex=1.8,y.intersp=1)}
dev.off()}

#---------------------------#
# 7. Distributional effects #
#---------------------------#
y1max=400
S=d-b
s=1
while(s<=S){
res=density(TRT[(b+s-a+1),],from=0,to=y1max)
x_axis=res[[1]]
f_trt=res[[2]]
res=density(SYN[(b+s-a+1),],from=0,to=y1max)
f_syn=res[[2]]
var=paste0("x_axis",s)
assign(var,x_axis)
var=paste0("f_trt",s)
assign(var,f_trt)
var=paste0("f_syn",s)
assign(var,f_syn)
s=s+1}

#####
##if(case==1){
##{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/DistEffects",case,".pdf",sep=""))
##par(mfrow=c(1,4),mar=c(2,2,2,2))
##mt1=paste0("Case ",case,":"," ",Time[b+1]) 
##matplot(x_axis1,cbind(f_trt1,f_syn1),main=mt1, cex.main=1.5,type="l",lty=c(1,2),col=c(2,4),lwd=c(3,3),xlab=" ",ylab=" ",xaxt="n",yaxt="n")
##axis(1,cex.axis=1.5)
##axis(2,cex.axis=1.5)
##legend("topleft",legend=c("actual","counterfactual"),lty=c(1,2),col=c(2,4),lwd=c(3,3),bty="n",cex=1,y.intersp=2)
##mt2=paste0("Case ",case,":"," ",Time[b+2]) 
##matplot(x_axis2,cbind(f_trt2,f_syn2),main=mt2, cex.main=1.5,type="l",lty=c(1,2),col=c(2,4),lwd=c(3,3),xlab=" ",ylab=" ",xaxt="n",yaxt="n")
##axis(1,cex.axis=1.5)
##axis(2,cex.axis=1.5)
##matplot(0,type="n",axes=F,xlab="",ylab="")
##matplot(0,type="n",axes=F,xlab="",ylab="")
##par(mfrow=c(1,1))
##dev.off()}}
##
##if(case>1){
##{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/DistEffects",case,".pdf",sep=""))
##par(mfrow=c(1,4),mar=c(2,2,2,2))
##mt1=paste0("Case ",case,":"," ",Time[b+1]) 
##matplot(x_axis1,cbind(f_trt1,f_syn1),main=mt1, cex.main=1.5,type="l",lty=c(1,2),col=c(2,4),lwd=c(3,3),xlab=" ",ylab=" ",xaxt="n",yaxt="n")
##axis(1,cex.axis=1.5)
##axis(2,cex.axis=1.5)
##mt2=paste0("Case ",case,":"," ",Time[b+2]) 
##matplot(x_axis2,cbind(f_trt2,f_syn2),main=mt2, cex.main=1.5,type="l",lty=c(1,2),col=c(2,4),lwd=c(3,3),xlab=" ",ylab=" ",xaxt="n",yaxt="n")
##axis(1,cex.axis=1.5)
##axis(2,cex.axis=1.5)
##mt3=paste0("Case ",case,":"," ",Time[b+3]) 
##matplot(x_axis2,cbind(f_trt3,f_syn3),main=mt3, cex.main=1.5,type="l",lty=c(1,2),col=c(2,4),lwd=c(3,3),xlab=" ",ylab=" ",xaxt="n",yaxt="n")
##axis(1,cex.axis=1.5)
##axis(2,cex.axis=1.5)
##mt4=paste0("Case ",case,":"," ",Time[b+4]) 
##matplot(x_axis2,cbind(f_trt4,f_syn4),main=mt4, cex.main=1.5,type="l",lty=c(1,2),col=c(2,4),lwd=c(3,3),xlab=" ",ylab=" ",xaxt="n",yaxt="n")
##axis(1,cex.axis=1.5)
##axis(2,cex.axis=1.5)
##par(mfrow=c(1,1))
##dev.off()}}

#------------------------------------#
# 8. Quantile effects & placebo test #
#------------------------------------#
#############################
### Table 2: quantile effects (part A)
s=1
while(s<=S){
Qdiff=TRT[(b+s-a+1),]-SYN[(b+s-a+1),]      # quantile effects
plb_Qdiff=matrix(0,length(tau),J)  # placebo quantile effects
j=1
while(j<=J){
plb_TRT=get(paste0("plb_TRT_",case,"_",j))
plb_SYN=get(paste0("plb_SYN_",case,"_",j))
plb_Qdiff[,j]=plb_TRT[(b+s-a+1),]-plb_SYN[(b+s-a+1),]
j=j+1}
var=paste0("Qdiff",s)
assign(var,Qdiff)
var=paste0("plb_Qdiff",s)
assign(var,plb_Qdiff)
Delta_Qdiff10=Qdiff[10]/SYN[(b+s-a+1),10]
Delta_Qdiff50=Qdiff[50]/SYN[(b+s-a+1),50]
Delta_Qdiff90=Qdiff[90]/SYN[(b+s-a+1),90]
Delta_Emw=Size_Emwch/Pre_Emw
Elast_Qdiff10=Delta_Qdiff10/Delta_Emw
Elast_Qdiff50=Delta_Qdiff50/Delta_Emw
Elast_Qdiff90=Delta_Qdiff90/Delta_Emw
cat("case=",case,"s=",s,"quantile effect on earnings","level01=",round(Qdiff[10],3),"elas01=",round(Elast_Qdiff10,3),"level05=",round(Qdiff[50],3),"elas05=",round(Elast_Qdiff50,3),"level09=",round(Qdiff[90],3),"elas09=",round(Elast_Qdiff90,3), "\n")
s=s+1}
cat("case=",case,"actual 10% quantile earnings","\n")
print(round(TRT[(b-a+1):(d-a+1),10],3))
cat("case=",case,"actual 50% quantile earnings","\n")
print(round(TRT[(b-a+1):(d-a+1),50],3))
cat("case=",case,"actual 90% quantile earnings","\n")
print(round(TRT[(b-a+1):(d-a+1),90],3))

### Figure 7 (Panel A)
if(case==1){
xtau=cbind(matrix(1,length(tau),1),tau,tau^2,tau^3)
ytau=Qdiff1
beta=solve(t(xtau)%*%xtau)%*%(t(xtau)%*%ytau)
fit1=xtau%*%beta
ytau=Qdiff2
beta=solve(t(xtau)%*%xtau)%*%(t(xtau)%*%ytau)
fit2=xtau%*%beta
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/QuantileEffects",case,".pdf",sep=""))
ymin=min(cbind(Qdiff1,Qdiff2))-10
ymax=max(cbind(Qdiff1,Qdiff2))+10
par(mfrow=c(1,4),mar=c(2,2,2,2))
mt1=paste0("Case ",case,":"," ",Time[b+1]) 
matplot(tau,plb_Qdiff1,main=mt1, cex.main=1.8,ylim=c(ymin,ymax),col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
lines(tau,Qdiff1,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
lines(tau,fit1,col=4,type="l",lty=2,lwd=4,xlab=" ",ylab=" ")
abline(h=Meffect[1],col=1,lty=1,lwd=2,xlab=" ",ylab=" ")
abline(h=0,lty=5,lwd=2,col="darkblue")
legend("bottomleft",legend=c("treated","placebo","fit","mean"),lty=c(1,1,2,1),col=c(2,gray(0.7),4,1),lwd=c(3,1,4,2),bty="n",cex=1.7,seg.len=0.5,y.intersp=1.2)
mt2=paste0("Case ",case,":"," ",Time[b+2]) 
matplot(tau,plb_Qdiff2,main=mt2, cex.main=1.5,ylim=c(ymin,ymax),col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
lines(tau,Qdiff2,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
lines(tau,fit2,col=4,type="l",lty=2,lwd=4,xlab=" ",ylab=" ")
abline(h=Meffect[2],col=1,lty=1,lwd=2,xlab=" ",ylab=" ")
abline(h=0,lty=5,lwd=2,col="darkblue")
matplot(0,type='n',axes=FALSE,ann=FALSE)
matplot(0,type='n',axes=FALSE,ann=FALSE)
par(mfrow=c(1,1))
dev.off()}}

if(case>1){
xtau=cbind(matrix(1,length(tau),1),tau,tau^2,tau^3)
ytau=Qdiff1
beta=solve(t(xtau)%*%xtau)%*%(t(xtau)%*%ytau)
fit1=xtau%*%beta
ytau=Qdiff2
beta=solve(t(xtau)%*%xtau)%*%(t(xtau)%*%ytau)
fit2=xtau%*%beta
ytau=Qdiff3
beta=solve(t(xtau)%*%xtau)%*%(t(xtau)%*%ytau)
fit3=xtau%*%beta
ytau=Qdiff4
beta=solve(t(xtau)%*%xtau)%*%(t(xtau)%*%ytau)
fit4=xtau%*%beta
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/QuantileEffects",case,".pdf",sep=""))
ymin=min(cbind(Qdiff1,Qdiff2,Qdiff3,Qdiff4))-10
ymax=max(cbind(Qdiff1,Qdiff2,Qdiff3,Qdiff4))+10
par(mfrow=c(1,4),mar=c(2,2,2,2))
mt1=paste0("Case ",case,":"," ",Time[b+1]) 
matplot(tau,plb_Qdiff1,main=mt1, cex.main=1.5,ylim=c(ymin,ymax),col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
lines(tau,Qdiff1,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
lines(tau,fit1,col=4,type="l",lty=2,lwd=4,xlab=" ",ylab=" ")
abline(h=Meffect[1],col=1,lty=1,lwd=2,xlab=" ",ylab=" ")
abline(h=0,lty=5,lwd=2,col="darkblue")
mt2=paste0("Case ",case,":"," ",Time[b+2]) 
matplot(tau,plb_Qdiff2,main=mt2, cex.main=1.5,ylim=c(ymin,ymax),col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
lines(tau,Qdiff2,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
lines(tau,fit2,col=4,type="l",lty=2,lwd=4,xlab=" ",ylab=" ")
abline(h=Meffect[2],col=1,lty=1,lwd=2,xlab=" ",ylab=" ")
abline(h=0,lty=5,lwd=2,col="darkblue")
mt3=paste0("Case ",case,":"," ",Time[b+3]) 
matplot(tau,plb_Qdiff3,main=mt3, cex.main=1.5,ylim=c(ymin,ymax),col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
lines(tau,Qdiff3,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
lines(tau,fit3,col=4,type="l",lty=2,lwd=4,xlab=" ",ylab=" ")
abline(h=Meffect[3],col=1,lty=1,lwd=2,xlab=" ",ylab=" ")
abline(h=0,lty=5,lwd=2,col="darkblue")
mt4=paste0("Case ",case,":"," ",Time[b+4]) 
matplot(tau,plb_Qdiff4,main=mt4, cex.main=1.5,ylim=c(ymin,ymax),col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
lines(tau,Qdiff4,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
lines(tau,fit4,col=4,type="l",lty=2,lwd=4,xlab=" ",ylab=" ")
abline(h=Meffect[4],col=1,lty=1,lwd=2,xlab=" ",ylab=" ")
abline(h=0,lty=5,lwd=2,col="darkblue")
par(mfrow=c(1,1))
dev.off()}}

# Placebo tests for the quantile effects
plb_preMSE_SYN=matrix(0,length(tau),J)
plb_postMSE_SYN=matrix(0,length(tau),J)
j=1
while(j<=J){
plb_TRT=get(paste0("plb_TRT_",case,"_",j))
plb_SYN=get(paste0("plb_SYN_",case,"_",j))
plb_trt_name=Sname[pcu_indx[j]]
q=1
while(q<=length(tau)){
plb_preMSE_SYN[q,j]=mean((plb_TRT[1:(b-a+1),q]-plb_SYN[1:(b-a+1),q])^2)
plb_postMSE_SYN[q,j]=mean((plb_TRT[(c-a+1):(d-a+1),q]-plb_SYN[(c-a+1):(d-a+1),q])^2)
q=q+1}
j=j+1}

# P-value of the rank statistic of the postMSE_SYN/preMSE_SYN
R_SYN=postMSE_SYN/preMSE_SYN
plb_R_SYN=plb_postMSE_SYN/plb_preMSE_SYN
R=cbind(R_SYN,plb_R_SYN)
Rrank=matrix(0,length(tau),1)
i=1
while(i<=length(tau)){
Rrank[i]=rank(R[i,])
i=i+1}
Pvalue=1-Rrank/(J+1)
cat("case=",case,"quantile effects on earnings","p-value=",round(Pvalue[10],3),round(Pvalue[50],3),round(Pvalue[90],3),"\n")

### Figure A.13 (Panel A)
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/placebotest",case,".pdf",sep=""))
plot(tau,Pvalue,main=mt, cex.main=1.5,xlab=" ",ylab=" ",type="l",ylim=c(0,1),col=2,lty=c(1,2),lwd=2,xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
abline(h=0.1,lty=5,lwd=1,col="blue")
abline(h=0.05,lty=5,lwd=1,col="darkblue")
dev.off()}

#####################################################################

##########################################
# II. Results for the Employment Effects # 
##########################################

#---------------------------#
# 1. Load empirical results #
#---------------------------#
pcu_indx0 = read.csv(paste0(dirpath, "Data/Results/", "pcu_indx_", case, ".csv"))
TRT0=read.csv(paste0(dirpath,"Data/Results/","TRTemp_",case,".csv"))
SYN0=read.csv(paste0(dirpath,"Data/Results/","SYNemp_",case,".csv"))
SYN00=read.csv(paste0(dirpath,"Data/Results/","SYN0emp_",case,".csv"))
wSYN0=read.csv(paste0(dirpath,"Data/Results/","wSYNemp_",case,".csv"))
wSYN00=read.csv(paste0(dirpath,"Data/Results/","wSYN0emp_",case,".csv"))
trtPre_fit0=read.csv(paste0(dirpath,"Data/Results/","trtPreemp_fit_",case,".csv"))
pcuPre_fit0=read.csv(paste0(dirpath,"Data/Results/","pcuPreemp_fit_",case,".csv"))
pcu_indx=as.matrix(pcu_indx0[,-1])  # Remove the 1st column (ordering number)
TRT=as.matrix(TRT0[,-1])
SYN=as.matrix(SYN0[,-1])
SYN0=as.matrix(SYN00[,-1])
wSYN=as.matrix(wSYN0[,-1])
wSYN0=as.matrix(wSYN00[,-1])
trtPre_fit=as.matrix(trtPre_fit0[,-1])
pcuPre_fit=as.matrix(pcuPre_fit0[,-1])
J=length(pcu_indx)
j=1
while(j<=J){
var=paste0("plb_pcu_indx_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plbemp_pcu_indx_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_TRT_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plbemp_TRT_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_SYN_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plbemp_SYN_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_SYN0_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plbemp_SYN0_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_wSYN_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plbemp_wSYN_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_wSYN0_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plbemp_wSYN0_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_trtPre_fit_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plbemp_trtPre_fit_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
var=paste0("plb_pcuPre_fit_",case,"_",j)
var0=read.csv(paste0(dirpath,"Data/Results/","plbemp_pcuPre_fit_",case,"_",j,".csv"))
assign(var,as.matrix(var0[,-1]))
j=j+1}


#--------------------#
# 2. Optimal weights #
#--------------------#
### Figure 2  (Panel B)
# mean weight
mt=paste0("Case ",case,":"," ",trt_name) 
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/emp_MeanWeights",case,".pdf",sep=""))
optw=rowMeans(wSYN)
bar=barplot(as.vector(optw), main=mt, cex.main=1.5, ylim=c(0,1), col="grey", angle=45, 
border=TRUE,las=2,cex.names=1.2)
text(bar,as.vector(optw)+0.05,round(as.vector(optw),2),cex=1.2)
text(bar,BSname[pcu_indx],y = -0.01,cex = 1.1,srt = 60, xpd = TRUE, pos = 2 )
dev.off()}



### Figure A.12 (Panel B)
# quantile weights
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/emp_QuantileWeights",case,".pdf",sep=""))
cr=rainbow(length(tau),s=1,v=1,start=0,end=0.9)
lg=tau
lcr=c(cr[10],cr[20],cr[30],cr[40],cr[50],cr[60],cr[70],cr[80],cr[90])
llg=c(lg[10],lg[20],lg[30],lg[40],lg[50],lg[60],lg[70],lg[80],lg[90])
optw=rowMeans(wSYN)
bar=barplot(as.vector(optw), main=mt, cex.main=1.5, ylim=c(0,1), col=NA, angle=45, 
border=FALSE,las=2,cex.names=1.2)
text(bar,BSname[pcu_indx],y = -0.01,cex = 1.1,srt = 60, xpd = TRUE, pos = 2 )
i=1
while(i<=length(tau)){
points(x=bar,wSYN[,i],col=cr[i],type="l")
i=i+1}
#legend("topleft",legend=llg,lty=1,col=lcr,bty="n",lwd=1,cex=1.2,y.intersp=1)
dev.off()}

#--------------------------------------------------------------------------------------#
# 3. Compare the economic features of the treated unit and the synthetic control units #
#--------------------------------------------------------------------------------------#
if(case==1){XData=read.csv(paste0(dirpath,"Data/","US_state_comparison1996.csv"))}
if(case==2){XData=read.csv(paste0(dirpath,"Data/","US_state_comparison1996.csv"))}
if(case==3){XData=read.csv(paste0(dirpath,"Data/","US_state_comparison1998.csv"))}
if(case==4){XData=read.csv(paste0(dirpath,"Data/","US_state_comparison2004.csv"))}
XData=as.matrix(XData[,-1])
X0_pop=XData[1,]       # Resident population
X0_white=XData[2,]     # White alone, percentage
X0_old=XData[3,]       # 65 years old and over, percentage
X0_land=XData[4,]      # Land area
X0_pden=XData[5,]      # Population per square mile (Preclude this from the matching condition because it is redundant)
X0_pov=XData[6,]       # Person below poverty level, percentage
X0_pinc=XData[7,]      # Personal income per capita
X0_hinc=XData[8,]      # Real median household income
X0_emp=XData[9,]       # Annual nonfarm employment
X0_nhu=XData[10,]      # New private housing units authorized by buliding permits
X0econ=as.matrix(rbind(X0_pop,X0_white,X0_old,X0_land,X0_pov,X0_pinc,X0_hinc,X0_emp,X0_nhu))

trt_indx=Sindx[which(Sname==trt_name)]
Xecon_trt=X0econ[,trt_indx]
Xecon_pcu=X0econ[,pcu_indx]
optw=as.matrix(rowMeans(wSYN))
Xecon_syn=Xecon_pcu%*%optw
optw01=as.matrix(wSYN[,10])
Xecon_syn01=Xecon_pcu%*%optw01
optw05=as.matrix(wSYN[,50])
Xecon_syn05=Xecon_pcu%*%optw05
optw09=as.matrix(wSYN[,90])
Xecon_syn09=Xecon_pcu%*%optw09
Distance=matrix(0,length(pcu_indx),1)
i=1
while(i<=nrow(Distance)){
Distance[i]=(mean((Xecon_trt-Xecon_pcu[,i])^2))^0.5
i=i+1}
singleindx=pcu_indx[which.min(Distance)]
singlename=Sname[singleindx]
Xecon_single=X0econ[,singleindx]

### Table 2 (Panel B)
print("wage effects")
print(singlename)
print(cbind(round(Xecon_trt,2),round(Xecon_syn,2),round(Xecon_syn01,2),round(Xecon_syn05,2),round(Xecon_syn09,2),round(Xecon_single,2)))


#---------------------------#
# 4. Performance evaluation #
#---------------------------#
# Mean sqaured (prediction) errors
preMSE_SYN=matrix(0,length(tau),1)
preMSE_SYN0=matrix(0,length(tau),1)
postMSE_SYN=matrix(0,length(tau),1)
postMSE_SYN0=matrix(0,length(tau),1)
q=1
while(q<=length(tau)){
tauq=tau[q]
preMSE_SYN[q]=mean((TRT[1:(b-a+1),q]-SYN[1:(b-a+1),q])^2)
preMSE_SYN0[q]=mean((TRT[1:(b-a+1),q]-SYN0[1:(b-a+1),q])^2)
postMSE_SYN[q]=mean((TRT[(c-a+1):(d-a+1),q]-SYN[(c-a+1):(d-a+1),q])^2)
postMSE_SYN0[q]=mean((TRT[(c-a+1):(d-a+1),q]-SYN0[(c-a+1):(d-a+1),q])^2)
q=q+1}

### Figure 3 (Panel B)
if(case==1){
mt=paste0("Case ",case,":"," ",trt_name) 
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/emp_preMSE",case,".pdf",sep=""))
matplot(tau,cbind(preMSE_SYN,preMSE_SYN0),main=mt, cex.main=1.5, ylim=c(0,300),xlab=" ",ylab=" ",type="l",col=c(2,4),lty=c(1,2),lwd=c(3,3),xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
#legend("topleft",legend=c("proposed","conventional"),col=c(2,4),lty=c(1,2),lwd=c(3,3),bty="n",cex=1.5,y.intersp=2)
dev.off()}}
if(case>1 & case<4){
mt=paste0("Case ",case,":"," ",trt_name) 
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/emp_preMSE",case,".pdf",sep=""))
matplot(tau,cbind(preMSE_SYN,preMSE_SYN0),main=mt, cex.main=1.5, ylim=c(0,100),xlab=" ",ylab=" ",type="l",col=c(2,4),lty=c(1,2),lwd=c(3,3),xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
#legend("topleft",legend=c("proposed","conventional"),col=c(2,4),lty=c(1,2),lwd=c(3,3),bty="n",cex=1.5,y.intersp=2)
dev.off()}}
if(case==4){
mt=paste0("Case ",case,":"," ",trt_name) 
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/emp_preMSE",case,".pdf",sep=""))
matplot(tau,cbind(preMSE_SYN,preMSE_SYN0),main=mt, cex.main=1.5, ylim=c(0,100),xlab=" ",ylab=" ",type="l",col=c(2,4),lty=c(1,2),lwd=c(3,3),xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
#legend("topleft",legend=c("proposed","conventional"),col=c(2,4),lty=c(1,2),lwd=c(3,3),bty="n",cex=1.5,y.intersp=2)
dev.off()}}


#-----------------------------------------#
# 5. Mean & quantile intervention effects #
#-----------------------------------------#
T=nrow(TRT)
res=MGQ(TRT,tau)
mean_trt=res[[1]]
res=MGQ(SYN,tau)
mean_syn=res[[1]]
res=MGQ(SYN0,tau)
mean_syn0=res[[1]]

### Figure 5
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/emp_TimeSeriesEffect",case,".pdf",sep=""))
par(mfrow=c(1,4),mar=c(2,2.5,2.5,2))
ymin=min(cbind(TRT[,10],SYN[,10],SYN0[,10]))
ymax=max(cbind(TRT[,90],SYN[,90],SYN0[,90]))

if(case==1){mtm=expression("Case1: mean")}
if(case==2){mtm=expression("Case2: mean")}
if(case==3){mtm=expression("Case3: mean")}
if(case==4){mtm=expression("Case4: mean")}
matplot(Xtime,cbind(mean_trt,mean_syn,mean_syn0),main=mtm, cex.main=1.5,ylim=c(ymin,ymax),type="l",lty=c(1,2,4),col=c(2,4,1),lwd=c(3,3,3),pch=19,cex=0.5,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
if(case==1){legend("topleft",legend=c("actual","proposed","conventional"),lty=c(1,2,4),col=c(2,4,1),lwd=c(3,3,3),bty="n",cex=1.2,seg.len=0.5,y.intersp=1.3)}
abline(v=Xtime[(c-a+1)],lty=5,lwd=2,col="darkblue")
axis(1,cex.axis=1.5,at=Xaxis,labels=x_lab)
axis(2,cex.axis=1.5)

if(case==1){mt1=expression("Case1:" ~ tau ~ "=0.1")}
if(case==2){mt1=expression("Case2:" ~ tau ~ "=0.1")}
if(case==3){mt1=expression("Case3:" ~ tau ~ "=0.1")}
if(case==4){mt1=expression("Case4:" ~ tau ~ "=0.1")}
matplot(Xtime,cbind(TRT[,10],SYN[,10],SYN0[,10]),main=mt1, cex.main=1.5,ylim=c(ymin,ymax),type="l",lty=c(1,2,4),col=c(2,4,1),lwd=c(3,3,3),pch=19,cex=0.5,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
abline(v=Xtime[(c-a+1)],lty=5,lwd=2,col="darkblue")
axis(1,cex.axis=1.5,at=Xaxis,labels=x_lab)
axis(2,cex.axis=1.5)

if(case==1){mt5=expression("Case1:" ~ tau ~ "=0.5")}
if(case==2){mt5=expression("Case2:" ~ tau ~ "=0.5")}
if(case==3){mt5=expression("Case3:" ~ tau ~ "=0.5")}
if(case==4){mt5=expression("Case4:" ~ tau ~ "=0.5")}
matplot(Xtime,cbind(TRT[,50],SYN[,50],SYN0[,50]),main=mt5, cex.main=1.5,ylim=c(ymin,ymax),type="l",lty=c(1,2,4),col=c(2,4,1),lwd=c(3,3,3),pch=19,cex=0.5,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
abline(v=Xtime[(c-a+1)],lty=5,lwd=2,col="darkblue")
axis(1,cex.axis=1.5,at=Xaxis,labels=x_lab)
axis(2,cex.axis=1.5)

if(case==1){mt9=expression("Case1:" ~ tau ~ "=0.9")}
if(case==2){mt9=expression("Case2:" ~ tau ~ "=0.9")}
if(case==3){mt9=expression("Case3:" ~ tau ~ "=0.9")}
if(case==4){mt9=expression("Case4:" ~ tau ~ "=0.9")}
matplot(Xtime,cbind(TRT[,90],SYN[,90],SYN0[,90]),main=mt9, cex.main=1.5,ylim=c(ymin,ymax),type="l",lty=c(1,2,4),col=c(2,4,1),lwd=c(3,3,3),pch=19,cex=0.5,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
abline(v=Xtime[(c-a+1)],lty=5,lwd=2,col="darkblue")
axis(1,cex.axis=1.5,at=Xaxis,labels=x_lab)
axis(2,cex.axis=1.5)
par(mfrow=c(1,1))
dev.off()}




#-------------------------------#
# 6. Mean effect & placebo test #
#-------------------------------#
Mdiff=mean_trt-mean_syn
Rmean_SYN=mean(Mdiff[(c-a+1):(d-a+1)]^2)/mean(Mdiff[1:(b-a+1)]^2)
plb_Mdiff=matrix(0,T,J)
plb_Rmean_SYN=matrix(0,J,1)
j=1
while(j<=J){
plb_TRT=get(paste0("plb_TRT_",case,"_",j))
plb_SYN=get(paste0("plb_SYN_",case,"_",j))
res=MGQ(plb_TRT,tau)
mean_plbtrt=res[[1]]
res=MGQ(plb_SYN,tau)
mean_plbsyn=res[[1]]
plb_Mdiff[,j]=mean_plbtrt-mean_plbsyn
plb_Rmean_SYN[j]=mean(plb_Mdiff[(c-a+1):(d-a+1),j]^2)/mean(plb_Mdiff[1:(b-a+1),j]^2)
j=j+1}
Meffect=Mdiff[(c-a+1):(d-a+1)]
Delta_Mdiff=Meffect/mean_syn[(c-a+1):(d-a+1)]
Delta_Emw=Size_Emwch/Pre_Emw
Elast_Mdiff=Delta_Mdiff/Delta_Emw
Rmean=c(Rmean_SYN,plb_Rmean_SYN)
Rrank=rank(Rmean)
Pvalue=1-Rrank[1]/(J+1)

########################
### Table 2: mean effects (part B)
s=1
while(s<=length(Meffect)){
cat("case=",case,"mean effect on employments","s=",s,"level=",round(100*Meffect[s],0),"elasticity=",round(Elast_Mdiff[s],3), "\n")
s=s+1}
cat("case=",case,"mean effect on employments","p-value=",round(Pvalue,3),"\n")
# mean wages at To,To+1,...
cat("case=",case,"actual mean employments","\n")
print(round(100*mean_trt[(b-a+1):(d-a+1)],0))

### Figure 6 (Panel B)
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/emp_MeanPlacebo",case,".pdf",sep=""))
matplot(Xtime,plb_Mdiff,main=mt, cex.main=1.5,col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
lines(Xtime,Mdiff,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
abline(v=Xtime[(c-a+1)],lty=5,lwd=2,col="darkblue")
abline(h=0,lty=5,lwd=2,col="darkblue")
axis(1,cex.axis=1.5,at=Xaxis,labels=x_lab)
axis(2,cex.axis=1.5)
#if(case==1){legend("topleft",legend=c("treated","placebo"),lty=c(1,1),col=c(2,gray(0.5)),lwd=c(3,1),bty="n",cex=1.5,y.intersp=2)}
dev.off()}


#---------------------------#
# 7. Distributional effects #
#---------------------------#
y1max=300
if(case==1){y1max=400}
S=d-b
s=1
while(s<=S){
res=density(TRT[(b+s-a+1),],from=0,to=y1max)
x_axis=res[[1]]
f_trt=res[[2]]
res=density(SYN[(b+s-a+1),],from=0,to=y1max)
f_syn=res[[2]]
var=paste0("x_axis",s)
assign(var,x_axis)
var=paste0("f_trt",s)
assign(var,f_trt)
var=paste0("f_syn",s)
assign(var,f_syn)
s=s+1}

##### 
##if(case==1){
##{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/emp_DistEffects",case,".pdf",sep=""))
##par(mfrow=c(1,4),mar=c(2,2,2,2))
##mt1=paste0("Case ",case,":"," ",Time[b+1]) 
##matplot(x_axis1,cbind(f_trt1,f_syn1),main=mt1, cex.main=1.5,type="l",lty=c(1,2),col=c(2,4),lwd=c(3,3),xlab=" ",ylab=" ",xaxt="n",yaxt="n")
##axis(1,cex.axis=1.5)
##axis(2,cex.axis=1.5)
###legend("topleft",legend=c("actual","counterfactual"),lty=c(1,2),col=c(2,4),lwd=c(3,3),bty="n",cex=1,y.intersp=2)
##mt2=paste0("Case ",case,":"," ",Time[b+2]) 
##matplot(x_axis2,cbind(f_trt2,f_syn2),main=mt2, cex.main=1.5,type="l",lty=c(1,2),col=c(2,4),lwd=c(3,3),xlab=" ",ylab=" ",xaxt="n",yaxt="n")
##axis(1,cex.axis=1.5)
##axis(2,cex.axis=1.5)
##matplot(0,type="n",axes=F,xlab="",ylab="")
##matplot(0,type="n",axes=F,xlab="",ylab="")
##par(mfrow=c(1,1))
##dev.off()}}
##
##if(case>1){
##{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/emp_DistEffects",case,".pdf",sep=""))
##par(mfrow=c(1,4),mar=c(2,2,2,2))
##mt1=paste0("Case ",case,":"," ",Time[b+1]) 
##matplot(x_axis1,cbind(f_trt1,f_syn1),main=mt1, cex.main=1.5,type="l",lty=c(1,2),col=c(2,4),lwd=c(3,3),xlab=" ",ylab=" ",xaxt="n",yaxt="n")
##axis(1,cex.axis=1.5)
##axis(2,cex.axis=1.5)
##mt2=paste0("Case ",case,":"," ",Time[b+2]) 
##matplot(x_axis2,cbind(f_trt2,f_syn2),main=mt2, cex.main=1.5,type="l",lty=c(1,2),col=c(2,4),lwd=c(3,3),xlab=" ",ylab=" ",xaxt="n",yaxt="n")
##axis(1,cex.axis=1.5)
##axis(2,cex.axis=1.5)
##mt3=paste0("Case ",case,":"," ",Time[b+3]) 
##matplot(x_axis2,cbind(f_trt3,f_syn3),main=mt3, cex.main=1.5,type="l",lty=c(1,2),col=c(2,4),lwd=c(3,3),xlab=" ",ylab=" ",xaxt="n",yaxt="n")
##axis(1,cex.axis=1.5)
##axis(2,cex.axis=1.5)
##mt4=paste0("Case ",case,":"," ",Time[b+4]) 
##matplot(x_axis2,cbind(f_trt4,f_syn4),main=mt4, cex.main=1.5,type="l",lty=c(1,2),col=c(2,4),lwd=c(3,3),xlab=" ",ylab=" ",xaxt="n",yaxt="n")
##axis(1,cex.axis=1.5)
##axis(2,cex.axis=1.5)
##par(mfrow=c(1,1))
##dev.off()}}

#------------------------------------#
# 8. Quantile effects & placebo test #
#------------------------------------#
#############################
### Table 2: quantile effects (part B)
s=1
while(s<=S){
Qdiff=TRT[(b+s-a+1),]-SYN[(b+s-a+1),]      # quantile effects
plb_Qdiff=matrix(0,length(tau),J)  # placebo quantile effects
j=1
while(j<=J){
plb_TRT=get(paste0("plb_TRT_",case,"_",j))
plb_SYN=get(paste0("plb_SYN_",case,"_",j))
plb_Qdiff[,j]=plb_TRT[(b+s-a+1),]-plb_SYN[(b+s-a+1),]
j=j+1}
var=paste0("Qdiff",s)
assign(var,Qdiff)
var=paste0("plb_Qdiff",s)
assign(var,plb_Qdiff)
Delta_Qdiff10=Qdiff[10]/SYN[(b+s-a+1),10]
Delta_Qdiff50=Qdiff[50]/SYN[(b+s-a+1),50]
Delta_Qdiff90=Qdiff[90]/SYN[(b+s-a+1),90]
Delta_Emw=Size_Emwch/Pre_Emw
Elast_Qdiff10=Delta_Qdiff10/Delta_Emw
Elast_Qdiff50=Delta_Qdiff50/Delta_Emw
Elast_Qdiff90=Delta_Qdiff90/Delta_Emw
cat("case=",case,"s=",s,"quantile effect on employments","level01=",round(100*Qdiff[10],0),"elas01=",round(Elast_Qdiff10,3),"level05=",round(100*Qdiff[50],0),"elas05=",round(Elast_Qdiff50,3),"level09=",round(100*Qdiff[90],0),"elas09=",round(Elast_Qdiff90,3), "\n")
s=s+1}
cat("case=",case,"actual 10% quantile employments","\n")
print(round(100*TRT[(b-a+1):(d-a+1),10],0))
cat("case=",case,"actual 50% quantile employments","\n")
print(round(100*TRT[(b-a+1):(d-a+1),50],0))
cat("case=",case,"actual 90% quantile employments","\n")
print(round(100*TRT[(b-a+1):(d-a+1),90],0))


### Figure 7 (Panel B)
if(case==1){
xtau=cbind(matrix(1,length(tau),1),tau,tau^2,tau^3)
ytau=Qdiff1
beta=solve(t(xtau)%*%xtau)%*%(t(xtau)%*%ytau)
fit1=xtau%*%beta
ytau=Qdiff2
beta=solve(t(xtau)%*%xtau)%*%(t(xtau)%*%ytau)
fit2=xtau%*%beta
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/emp_QuantileEffects",case,".pdf",sep=""))
ymin=min(cbind(Qdiff1,Qdiff2))-10
ymax=max(cbind(Qdiff1,Qdiff2))+10
par(mfrow=c(1,4),mar=c(2,2,2,2))
mt1=paste0("Case ",case,":"," ",Time[b+1]) 
matplot(tau,plb_Qdiff1,main=mt1, cex.main=1.5,ylim=c(ymin,ymax),col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
lines(tau,Qdiff1,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
lines(tau,fit1,col=4,type="l",lty=2,lwd=4,xlab=" ",ylab=" ")
abline(h=Meffect[1],col=1,lty=1,lwd=2,xlab=" ",ylab=" ")
abline(h=0,lty=5,lwd=2,col="darkblue")
#legend("topleft",legend=c("treated","placebo","fit","mean"),lty=c(1,1,2,3),col=c(2,gray(0.7),4,1),lwd=c(3,1,4,4),bty="n",cex=1,y.intersp=2)
mt2=paste0("Case ",case,":"," ",Time[b+2]) 
matplot(tau,plb_Qdiff2,main=mt2, cex.main=1.5,ylim=c(ymin,ymax),col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
lines(tau,Qdiff2,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
lines(tau,fit2,col=4,type="l",lty=2,lwd=4,xlab=" ",ylab=" ")
abline(h=Meffect[2],col=1,lty=1,lwd=2,xlab=" ",ylab=" ")
abline(h=0,lty=5,lwd=2,col="darkblue")
matplot(0,type='n',axes=FALSE,ann=FALSE)
matplot(0,type='n',axes=FALSE,ann=FALSE)
par(mfrow=c(1,1))
dev.off()}}

if(case>1){
xtau=cbind(matrix(1,length(tau),1),tau,tau^2,tau^3)
ytau=Qdiff1
beta=solve(t(xtau)%*%xtau)%*%(t(xtau)%*%ytau)
fit1=xtau%*%beta
ytau=Qdiff2
beta=solve(t(xtau)%*%xtau)%*%(t(xtau)%*%ytau)
fit2=xtau%*%beta
ytau=Qdiff3
beta=solve(t(xtau)%*%xtau)%*%(t(xtau)%*%ytau)
fit3=xtau%*%beta
ytau=Qdiff4
beta=solve(t(xtau)%*%xtau)%*%(t(xtau)%*%ytau)
fit4=xtau%*%beta
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/emp_QuantileEffects",case,".pdf",sep=""))
ymin=min(cbind(Qdiff1,Qdiff2,Qdiff3,Qdiff4))-10
ymax=max(cbind(Qdiff1,Qdiff2,Qdiff3,Qdiff4))+10
par(mfrow=c(1,4),mar=c(2,2,2,2))
mt1=paste0("Case ",case,":"," ",Time[b+1]) 
matplot(tau,plb_Qdiff1,main=mt1, cex.main=1.5,ylim=c(ymin,ymax),col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
lines(tau,Qdiff1,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
lines(tau,fit1,col=4,type="l",lty=2,lwd=4,xlab=" ",ylab=" ")
abline(h=Meffect[1],col=1,lty=1,lwd=2,xlab=" ",ylab=" ")
abline(h=0,lty=5,lwd=2,col="darkblue")
mt2=paste0("Case ",case,":"," ",Time[b+2]) 
matplot(tau,plb_Qdiff2,main=mt2, cex.main=1.5,ylim=c(ymin,ymax),col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
lines(tau,Qdiff2,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
lines(tau,fit2,col=4,type="l",lty=2,lwd=4,xlab=" ",ylab=" ")
abline(h=Meffect[2],col=1,lty=1,lwd=2,xlab=" ",ylab=" ")
abline(h=0,lty=5,lwd=2,col="darkblue")
mt3=paste0("Case ",case,":"," ",Time[b+3]) 
matplot(tau,plb_Qdiff3,main=mt3, cex.main=1.5,ylim=c(ymin,ymax),col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
lines(tau,Qdiff3,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
lines(tau,fit3,col=4,type="l",lty=2,lwd=4,xlab=" ",ylab=" ")
abline(h=Meffect[3],col=1,lty=1,lwd=2,xlab=" ",ylab=" ")
abline(h=0,lty=5,lwd=2,col="darkblue")
mt4=paste0("Case ",case,":"," ",Time[b+4]) 
matplot(tau,plb_Qdiff4,main=mt4, cex.main=1.5,ylim=c(ymin,ymax),col=gray(0.7),type="l",lwd=1,lty=1,xlab=" ",ylab=" ",xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
lines(tau,Qdiff4,col=2,type="l",lwd=3,xlab=" ",ylab=" ")
lines(tau,fit4,col=4,type="l",lty=2,lwd=4,xlab=" ",ylab=" ")
abline(h=Meffect[4],col=1,lty=1,lwd=2,xlab=" ",ylab=" ")
abline(h=0,lty=5,lwd=2,col="darkblue")
par(mfrow=c(1,1))
dev.off()}}


# Placebo tests for the quantile effects
plb_preMSE_SYN=matrix(0,length(tau),J)
plb_postMSE_SYN=matrix(0,length(tau),J)
j=1
while(j<=J){
plb_TRT=get(paste0("plb_TRT_",case,"_",j))
plb_SYN=get(paste0("plb_SYN_",case,"_",j))
plb_trt_name=Sname[pcu_indx[j]]
q=1
while(q<=length(tau)){
plb_preMSE_SYN[q,j]=mean((plb_TRT[1:(b-a+1),q]-plb_SYN[1:(b-a+1),q])^2)
plb_postMSE_SYN[q,j]=mean((plb_TRT[(c-a+1):(d-a+1),q]-plb_SYN[(c-a+1):(d-a+1),q])^2)
q=q+1}
j=j+1}

# P-value of the rank statistic of the postMSE_SYN/preMSE_SYN
R_SYN=postMSE_SYN/preMSE_SYN
plb_R_SYN=plb_postMSE_SYN/plb_preMSE_SYN
R=cbind(R_SYN,plb_R_SYN)
Rrank=matrix(0,length(tau),1)
i=1
while(i<=length(tau)){
Rrank[i]=rank(R[i,])
i=i+1}
Pvalue=1-Rrank/(J+1)
cat("case=",case,"quantile effects on earnings","p-value=",round(Pvalue[10],3),round(Pvalue[50],3),round(Pvalue[90],3),"\n")

### Figure A.13 (Panel B)
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/emp_placebotest",case,".pdf",sep=""))
plot(tau,Pvalue,main=mt, cex.main=1.5,xlab=" ",ylab=" ",type="l",ylim=c(0,1),col=2,lty=c(1,2),lwd=2,xaxt="n",yaxt="n")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5)
abline(h=0.1,lty=5,lwd=1,col="blue")
abline(h=0.05,lty=5,lwd=1,col="darkblue")
dev.off()}

case=case+1}


#################################
########### END HERE ############
#################################
