#====================================================================#
# This is the first R code used in the empirical part of the paper:
# "A Distributional Synthetic Control Method for Policy Evaluation"
#====================================================================#
rm(list=ls())
#library(Synth)
library(Matrix)
set.seed(123)

#--------------------------------------------------------------------------------------------------------------------------------------
#==========================#
# Part I. Define Functions # 
#==========================#
##########################################
# 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. Data Exploration and Events Identification # 
#=====================================================#

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

### Footnote 1:  Data source: Dube, Lester and Reich(2010, REStat)
# Minimum Wage Effects across State Borders: Estimates Using Contiguous Counties
# https://dataverse.harvard.edu/dataset.xhtml?persistentId=hdl:1902.1/15969

# 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;              

# Distribution evolution to be explored 
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)

#-------------------------------------------#
# 1. Define the policy-influencing variable #
#-------------------------------------------#
Y=cbind(wage_rest,emp_rest) 
y1=Y[,1]        # wages
y2=Y[,2]/100    # (rescaled) employments

#---------------------------------------------#
# 2. State indicies, names, min. wage changes #
#---------------------------------------------#
 Sindx=matrix(0,N,1)      # state indices
 Sname=matrix(0,N,1)      # state names 
 Scode=matrix(0,N,1)      # state codes
   Smw=matrix(0,T,N)      # state min.wage time series
 SmwCh=matrix(0,T,N)      # state min.wage change time series

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
Smw[,i]=nmw_sta[which(cfips==minci)]
t=2
while(t<=T){
if(Smw[t,i]==0 & Smw[(t-1),i]>0){Smw[t,i]=Smw[(t-1),i]}
if(Smw[t,i]>Smw[(t-1),i]){SmwCh[t,i]=1}
t=t+1}
i=i+1}

# Federal min-wage series
minci=min(cfips[which(sindx==1)])
Fmw=mw_fed[which(cfips==minci)]

#--------------------------------------#
# 3. Number of counties for each state #
#--------------------------------------#
Ncf=matrix(0,T,N)           
Ncfi=matrix(0,N,1)          
i=1
while(i<=N){
t=1
while(t<=T){
cf=cfips[which(sindx==i & time==Time[t])]
Ncf[t,i]=length(cf)
t=t+1}
Ncfi[i]=max(Ncf[,i])
cat(Sname[i],Ncfi[i],"\n")
i=i+1}

### Figure
Sname[8]="D. of Columbia"
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/number_counties.pdf",sep=""))
bar=barplot(as.vector(Ncfi), main=" ", ylim=c(0,275), col="grey",border=TRUE,names.arg=Sname,las=3,cex.names=0.7)
text(bar,as.vector(Ncfi)+10,as.vector(Ncfi),cex=0.6)
dev.off()}

#-------------------------------------------------------#
# 4. State sample sizes (after removing missing values) #
#-------------------------------------------------------#
Nit=matrix(0,T,N)             # (state,time)-specific sample size
Ni=matrix(0,N,1)              # state sample size, precluding missing values  
i=1
while(i<=N){
t=1
while(t<=T){
w=na.omit(Y[which(sindx==i & time==Time[t])])
Nit[t,i]=length(w)
t=t+1}
cat(Sname[i],max(Nit[,i]),"\n")
Ni[i]=min(Nit[,i])
i=i+1}

### Figure A.8
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/number_completeobs.pdf",sep=""))
bar=barplot(as.vector(Ni), main=" ", ylim=c(0,115), col="grey",border=TRUE, names.arg=Sname,las=3,cex.names=0.7)
text(bar,as.vector(Ni)+5,as.vector(Ni),cex=0.6)
abline(h=20,lty=5,lwd=2,col="darkblue")
dev.off()}

# set (minN+1) to be the min. sample size for estimation
minN=20         

#-----------------------------------------------------#
# 4-1. Plot: States' wage mean & quantile time series #
#-----------------------------------------------------#
tau=seq(0.01,0.99,0.01)     # quantile probability indices
res=X_axis(Time,12)         # X-axis for the time series plot  
Xtime=res[[1]]
Xaxis=res[[2]]
x_lab=res[[3]]
cr=rainbow(length(tau),s=1,v=1,start=0,end=0.9) # Rainbow colors for the quantile sequences

### Figure A.9
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/MQevol_wage.pdf",sep=""))
par(mfrow=c(8,4),mar=c(2,2,2,2))
i=1
while(i<=N){
wm=matrix(0,T,1)
wq=matrix(0,T,length(tau))
if(i==1){Sys.sleep(1)}
if(Ni[i]>minN){
mt=as.character(Sname[i])
t=1
while(t<=T){
w=y1[which(sindx==i & time==Time[t])]
wm[t]=mean(na.omit(w))               # mean
wq[t,]=quantile(na.omit(w),tau)      # quantile 
t=t+1}
matplot(Xtime,wq,main=Sname[i],type="l",lty=1,lwd=0.5,col=cr,cex=0.5,xlab=" ",ylab=" ",xaxt="n")
lines(Xtime,wm,type="l",col=1,lty=1,lwd=1)
axis(1,at=Xaxis,labels=x_lab)
}
write.csv(wq,file=paste0(dirpath,"Data/Results/","Q_",Sname[i],".csv"))
i=i+1}
par(mfrow=c(1,1))
dev.off()}

#-----------------------------------------------------------#
# 4-2. Plot: States' employment mean & quantile time series #
#-----------------------------------------------------------#
tau=seq(0.01,0.99,0.01)     # quantile probability indices
res=X_axis(Time,12)         # X-axis for the time series plot  
Xtime=res[[1]]
Xaxis=res[[2]]
x_lab=res[[3]]
cr=rainbow(length(tau),s=1,v=1,start=0,end=0.9) # Rainbow colors for the quantile sequences

### Figure A.10
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/MQevol_emp.pdf",sep=""))
par(mfrow=c(8,4),mar=c(2,2,2,2))
i=1
while(i<=N){
wm=matrix(0,T,1)
wq=matrix(0,T,length(tau))
if(i==1){Sys.sleep(1)}
if(Ni[i]>minN){
mt=as.character(Sname[i])
t=1
while(t<=T){
w=y2[which(sindx==i & time==Time[t])]
wm[t]=mean(na.omit(w))               # mean
wq[t,]=quantile(na.omit(w),tau)      # quantile 
t=t+1}
matplot(Xtime,wq,main=Sname[i],type="l",lty=1,lwd=0.5,col=cr,cex=0.5,xlab=" ",ylab=" ",xaxt="n")
lines(Xtime,wm,type="l",col=1,lwd=1)
axis(1,at=Xaxis,labels=x_lab)
}
write.csv(wq,file=paste0(dirpath,"Data/Results/","Qemp_",Sname[i],".csv"))
i=i+1}
par(mfrow=c(1,1))
dev.off()}

#--------------------------------------------------#
# 4-3. Plot: Fed and states' min. wage time series #
#--------------------------------------------------#
### Figure A.11 
{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/min_wage.pdf",sep=""))
par(mfrow=c(8,4),mar=c(2,2,2,2))
i=1
while(i<=N){
if(i==1){Sys.sleep(1)}
if(Ni[i]>minN){
plot(Xtime,Smw[,i],main=Sname[i],ylim=c(0,9),type="p",pch=19,col=4,cex=0.2,xlab=" ",ylab=" ",xaxt="n")
lines(Xtime,Fmw,type="l",lty=1,col=2,lwd=1)
axis(1,at=Xaxis,labels=x_lab)
#Sys.sleep(1)
}
i=i+1}
par(mfrow=c(1,1))
dev.off()}

#*******************************************************#
# 3. Identification of the state min-wage-change events #
#*******************************************************#

#--------------------------------------------#
# 1. Define the state min-wage-change events #
#--------------------------------------------#
statesel=matrix(0,N,1) # state-selection dummies (Ni[i]>minN)
i=1
while(i<=N){
if(Ni[i]>minN){
statesel[i]=1                      # the ith state satisfies the min.# of obs. restriction      
Smwchtime=Time[1]
t=2
while(t<=T-1){
if(SmwCh[t,i]==1){     
Smwchtime=c(Smwchtime,Time[t])}
t=t+1}
Smwchtime=c(Smwchtime,Time[T])     # state min. wage change times if the dim of this vector > 2
if(length(Smwchtime)>2){
j=1
while(j<=(length(Smwchtime)-2)){
# Define an event of state min. wage change
Stime=Smwchtime[j:(j+2)]
 PreChange=Stime[1]                           # change time of the previous event
    Change=Stime[2]                           # change time of this event
NextChange=Stime[3]                           # change time of the next event        
    chindx=which(Time==Change)                # change time index        
   Nchindx=which(Time==NextChange)            # next change time index        
  FmwchPre=Fmw[(chindx-1)]                    # federal min. wage before change
  FmwchPost=Fmw[chindx]                       # federal min. wage after change
  SmwchPre=Smw[(chindx-1),i]                  # state min. wage before change
  SmwchPost=Smw[chindx,i]                     # state min. wage after change
  EmwchPre=max(SmwchPre,FmwchPre)             # effective min. wage before change 
  EmwchPost=max(SmwchPost,FmwchPost)          # effective min. wage after change 
  Fmwchsize=FmwchPost-FmwchPre                # federal min. wage change size
  Smwchsize=SmwchPost-SmwchPre                # state min. wage change size
  Emwchsize=SmwchPost-EmwchPre                # state-level effective min. wage change size
  PreStart=PreChange                          # start of the pre-intervention period
    PreEnd=Time[(chindx-1)]                   # end of the pre-intervention period 
 PostStart=Change                             # start of the post-intervention period   
   if(NextChange<Time[T]){                    # end of the post-intervention period
   PostEnd=Time[(Nchindx)-1]}
   if(NextChange==Time[T]){
   PostEnd=Time[T]}
    Year1=as.numeric(substring(PreChange,1,4))
 Quarter1=as.numeric(substring(PreChange,6,6))
    Yeart=as.numeric(substring(Change,1,4))
 Quartert=as.numeric(substring(Change,6,6))
    Year2=as.numeric(substring(NextChange,1,4))
 Quarter2=as.numeric(substring(NextChange,6,6))
  Pre_Duration=(4*Yeart+Quartert)-(4*Year1+Quarter1)  # Duration of the pre-intervention period
 Post_Duration=(4*Year2+Quarter2)-(4*Yeart+Quartert)  # Duration of the post-intervention period
 # # Print the events of minimum wage increases
 #cat(Sname[i],j,"PreStart=",PreStart,"PreEnd=",PreEnd,"PreDura=",Pre_Duration,"Pre_Fmw=",FmwchPre,"Pre_Smw=",SmwchPre,"Pre_Emw=",EmwchPre,"PostStart=",PostStart,"PostEnd=",PostEnd,"PostDura=",Post_Duration,"Post_Fmw=",FmwchPost,"Post_Smw=",SmwchPost,"Post_Emw=",EmwchPost,"Size_Fmwch=",Fmwchsize,"Size_Smwch=",Smwchsize,"Size_Emwch=",Emwchsize,"\n")
 j=j+1}}}
i=i+1}
SSindx=Sindx[which(statesel==1)] # state-selection index (Ni[i]>minN)

#### Table A.3:
# California 1 PreStart= 199001 PreEnd= 199604 PreDura= 28 Pre_Fmw= 4.75 Pre_Smw= 4.25 Pre_Emw= 4.75 PostStart= 199701 PostEnd= 199702 PostDura= 2  Post_Fmw= 4.75 Post_Smw= 5    Post_Emw= 5    Size_Fmwch= 0   Size_Smwch= 0.75 Size_Emwch= 0.25California 1 Pre_mw= 4.25 Size_mwch= 0.75 PreStart= 199001 PreEnd= 199604 PostStart= 199701 PostEnd= 199702 PreDura= 28 PostDura= 2 
# California 2 PreStart= 199701 PreEnd= 199702 PreDura= 2  Pre_Fmw= 4.75 Pre_Smw= 5    Pre_Emw= 5    PostStart= 199703 PostEnd= 199704 PostDura= 2  Post_Fmw= 5.15 Post_Smw= 5.15 Post_Emw= 5.15 Size_Fmwch= 0.4 Size_Smwch= 0.15 Size_Emwch= 0.15California 2 Pre_mw= 5    Size_mwch= 0.15 PreStart= 199701 PreEnd= 199702 PostStart= 199703 PostEnd= 199704 PreDura= 2 PostDura= 2 
# California 3 PreStart= 199703 PreEnd= 199704 PreDura= 2  Pre_Fmw= 5.15 Pre_Smw= 5.15 Pre_Emw= 5.15 PostStart= 199801 PostEnd= 200004 PostDura= 12 Post_Fmw= 5.15 Post_Smw= 5.75 Post_Emw= 5.75 Size_Fmwch= 0   Size_Smwch= 0.6  Size_Emwch= 0.6 California 3 Pre_mw= 5.15 Size_mwch= 0.6  PreStart= 199703 PreEnd= 199704 PostStart= 199801 PostEnd= 200004 PreDura= 2 PostDura= 12 
# California 4 PreStart= 199801 PreEnd= 200004 PreDura= 12 Pre_Fmw= 5.15 Pre_Smw= 5.75 Pre_Emw= 5.75 PostStart= 200101 PostEnd= 200104 PostDura= 4  Post_Fmw= 5.15 Post_Smw= 6.25 Post_Emw= 6.25 Size_Fmwch= 0   Size_Smwch= 0.5  Size_Emwch= 0.5 California 4 Pre_mw= 5.75 Size_mwch= 0.5  PreStart= 199801 PreEnd= 200004 PostStart= 200101 PostEnd= 200104 PreDura= 12 PostDura= 4 
# California 5 PreStart= 200101 PreEnd= 200104 PreDura= 4  Pre_Fmw= 5.15 Pre_Smw= 6.25 Pre_Emw= 6.25 PostStart= 200201 PostEnd= 200602 PostDura= 17 Post_Fmw= 5.15 Post_Smw= 6.75 Post_Emw= 6.75 Size_Fmwch= 0   Size_Smwch= 0.5  Size_Emwch= 0.5 California 5 Pre_mw= 6.25 Size_mwch= 0.5  PreStart= 200101 PreEnd= 200104 PostStart= 200201 PostEnd= 200602 PreDura= 4 PostDura= 17 
# Florida 1    PreStart= 199001 PreEnd= 200501 PreDura= 61 Pre_Fmw= 5.15 Pre_Smw= 0    Pre_Emw= 5.15 PostStart= 200502 PostEnd= 200504 PostDura= 3  Post_Fmw= 5.15 Post_Smw= 6.15 Post_Emw= 6.15 Size_Fmwch= 0   Size_Smwch= 6.15 Size_Emwch= 1   Florida    1 Pre_mw= 0    Size_mwch= 6.15 PreStart= 199001 PreEnd= 200501 PostStart= 200502 PostEnd= 200504 PreDura= 61 PostDura= 3 
# Florida 2    PreStart= 200502 PreEnd= 200504 PreDura= 3  Pre_Fmw= 5.15 Pre_Smw= 6.15 Pre_Emw= 6.15 PostStart= 200601 PostEnd= 200602 PostDura= 1  Post_Fmw= 5.15 Post_Smw= 6.4  Post_Emw= 6.4  Size_Fmwch= 0   Size_Smwch= 0.25 Size_Emwch= 0.25Florida    2 Pre_mw= 6.15 Size_mwch= 0.25 PreStart= 200502 PreEnd= 200504 PostStart= 200601 PostEnd= 200602 PreDura= 3 PostDura= 1 
# Illinois 1   PreStart= 199001 PreEnd= 200304 PreDura= 56 Pre_Fmw= 5.15 Pre_Smw= 0    Pre_Emw= 5.15 PostStart= 200401 PostEnd= 200404 PostDura= 4  Post_Fmw= 5.15 Post_Smw= 5.5  Post_Emw= 5.5  Size_Fmwch= 0   Size_Smwch= 5.5  Size_Emwch= 0.35Illinois   1 Pre_mw= 0    Size_mwch= 5.5  PreStart= 199001 PreEnd= 200304 PostStart= 200401 PostEnd= 200404 PreDura= 56 PostDura= 4 
# Illinois 2   PreStart= 200401 PreEnd= 200404 PreDura= 4  Pre_Fmw= 5.15 Pre_Smw= 5.5  Pre_Emw= 5.5  PostStart= 200501 PostEnd= 200602 PostDura= 5  Post_Fmw= 5.15 Post_Smw= 6.5  Post_Emw= 6.5  Size_Fmwch= 0   Size_Smwch= 1    Size_Emwch= 1   Illinois   2 Pre_mw= 5.5  Size_mwch= 1    PreStart= 200401 PreEnd= 200404 PostStart= 200501 PostEnd= 200602 PreDura= 4 PostDura= 5 
# Iowa 1       PreStart= 199001 PreEnd= 199004 PreDura= 4  Pre_Fmw= 3.8  Pre_Smw= 3.85 Pre_Emw= 3.85 PostStart= 199101 PostEnd= 199104 PostDura= 4  Post_Fmw= 3.8  Post_Smw= 4.25 Post_Emw= 4.25 Size_Fmwch= 0   Size_Smwch= 0.4  Size_Emwch= 0.4 Iowa       1 Pre_mw= 3.85 Size_mwch= 0.4  PreStart= 199001 PreEnd= 199004 PostStart= 199101 PostEnd= 199104 PreDura= 4 PostDura= 4 
# Iowa 2       PreStart= 199101 PreEnd= 199104 PreDura= 4  Pre_Fmw= 4.25 Pre_Smw= 4.25 Pre_Emw= 4.25 PostStart= 199201 PostEnd= 200602 PostDura= 57 Post_Fmw= 4.25 Post_Smw= 4.65 Post_Emw= 4.65 Size_Fmwch= 0   Size_Smwch= 0.4  Size_Emwch= 0.4 Iowa       2 Pre_mw= 4.25 Size_mwch= 0.4  PreStart= 199101 PreEnd= 199104 PostStart= 199201 PostEnd= 200602 PreDura= 4 PostDura= 57 
# Maryland 1   PreStart= 199001 PreEnd= 200504 PreDura= 64 Pre_Fmw= 5.15 Pre_Smw= 0    Pre_Emw= 5.15 PostStart= 200601 PostEnd= 200602 PostDura= 1  Post_Fmw= 5.15 Post_Smw= 6.15 Post_Emw= 6.15 Size_Fmwch= 0   Size_Smwch= 6.15 Size_Emwch= 1   Maryland   1 Pre_mw= 0    Size_mwch= 6.15 PreStart= 199001 PreEnd= 200504 PostStart= 200601 PostEnd= 200602 PreDura= 64 PostDura= 1 
# Minnesota 1  PreStart= 199001 PreEnd= 199004 PreDura= 4  Pre_Fmw= 3.8  Pre_Smw= 3.95 Pre_Emw= 3.95 PostStart= 199101 PostEnd= 200502 PostDura= 58 Post_Fmw= 3.8  Post_Smw= 4.25 Post_Emw= 4.25 Size_Fmwch= 0   Size_Smwch= 0.3  Size_Emwch= 0.3 Minnesota  1 Pre_mw= 3.95 Size_mwch= 0.3  PreStart= 199001 PreEnd= 199004 PostStart= 199101 PostEnd= 200502 PreDura= 4 PostDura= 58 
# Minnesota 2  PreStart= 199101 PreEnd= 200502 PreDura= 58 Pre_Fmw= 5.15 Pre_Smw= 4.25 Pre_Emw= 5.15 PostStart= 200503 PostEnd= 200602 PostDura= 3  Post_Fmw= 5.15 Post_Smw= 6.15 Post_Emw= 6.15 Size_Fmwch= 0   Size_Smwch= 1.9  Size_Emwch= 1   Minnesota  2 Pre_mw= 4.25 Size_mwch= 1.9  PreStart= 199101 PreEnd= 200502 PostStart= 200503 PostEnd= 200602 PreDura= 58 PostDura= 3 
# New York 1   PreStart= 199001 PreEnd= 200404 PreDura= 60 Pre_Fmw= 5.15 Pre_Smw= 0    Pre_Emw= 5.15 PostStart= 200501 PostEnd= 200504 PostDura= 4  Post_Fmw= 5.15 Post_Smw= 6    Post_Emw= 6    Size_Fmwch= 0   Size_Smwch= 6    Size_Emwch= 0.85New York   1 Pre_mw= 0    Size_mwch= 6    PreStart= 199001 PreEnd= 200404 PostStart= 200501 PostEnd= 200504 PreDura= 60 PostDura= 4 
# New York 2   PreStart= 200501 PreEnd= 200504 PreDura= 4  Pre_Fmw= 5.15 Pre_Smw= 6    Pre_Emw= 6    PostStart= 200601 PostEnd= 200602 PostDura= 1  Post_Fmw= 5.15 Post_Smw= 6.75 Post_Emw= 6.75 Size_Fmwch= 0   Size_Smwch= 0.75 Size_Emwch= 0.75New York   2 Pre_mw= 6    Size_mwch= 0.75 PreStart= 200501 PreEnd= 200504 PostStart= 200601 PostEnd= 200602 PreDura= 4 PostDura= 1 
# Oregon 1     PreStart= 199001 PreEnd= 199004 PreDura= 4  Pre_Fmw= 3.8  Pre_Smw= 4.25 Pre_Emw= 4.25 PostStart= 199101 PostEnd= 199604 PostDura= 24 Post_Fmw= 3.8  Post_Smw= 4.75 Post_Emw= 4.75 Size_Fmwch= 0   Size_Smwch= 0.5  Size_Emwch= 0.5 Oregon     1 Pre_mw= 4.25 Size_mwch= 0.5  PreStart= 199001 PreEnd= 199004 PostStart= 199101 PostEnd= 199604 PreDura= 4 PostDura= 24 
# Oregon 2     PreStart= 199101 PreEnd= 199604 PreDura= 24 Pre_Fmw= 4.75 Pre_Smw= 4.75 Pre_Emw= 4.75 PostStart= 199701 PostEnd= 199704 PostDura= 4  Post_Fmw= 4.75 Post_Smw= 5.5  Post_Emw= 5.5  Size_Fmwch= 0   Size_Smwch= 0.75 Size_Emwch= 0.75Oregon     2 Pre_mw= 4.75 Size_mwch= 0.75 PreStart= 199101 PreEnd= 199604 PostStart= 199701 PostEnd= 199704 PreDura= 24 PostDura= 4 
# Oregon 3     PreStart= 199701 PreEnd= 199704 PreDura= 4  Pre_Fmw= 5.15 Pre_Smw= 5.5  Pre_Emw= 5.5  PostStart= 199801 PostEnd= 199804 PostDura= 4  Post_Fmw= 5.15 Post_Smw= 6    Post_Emw= 6    Size_Fmwch= 0   Size_Smwch= 0.5  Size_Emwch= 0.5 Oregon     3 Pre_mw= 5.5  Size_mwch= 0.5  PreStart= 199701 PreEnd= 199704 PostStart= 199801 PostEnd= 199804 PreDura= 4 PostDura= 4 
# Oregon 4     PreStart= 199801 PreEnd= 199804 PreDura= 4  Pre_Fmw= 5.15 Pre_Smw= 6    Pre_Emw= 6    PostStart= 199901 PostEnd= 200204 PostDura= 16 Post_Fmw= 5.15 Post_Smw= 6.5  Post_Emw= 6.5  Size_Fmwch= 0   Size_Smwch= 0.5  Size_Emwch= 0.5 Oregon     4 Pre_mw= 6    Size_mwch= 0.5  PreStart= 199801 PreEnd= 199804 PostStart= 199901 PostEnd= 200204 PreDura= 4 PostDura= 16 
# Oregon 5     PreStart= 199901 PreEnd= 200204 PreDura= 16 Pre_Fmw= 5.15 Pre_Smw= 6.5  Pre_Emw= 6.5  PostStart= 200301 PostEnd= 200504 PostDura= 12 Post_Fmw= 5.15 Post_Smw= 6.9  Post_Emw= 6.9  Size_Fmwch= 0   Size_Smwch= 0.4  Size_Emwch= 0.4 Oregon     5 Pre_mw= 6.5  Size_mwch= 0.4  PreStart= 199901 PreEnd= 200204 PostStart= 200301 PostEnd= 200504 PreDura= 16 PostDura= 12 
# Oregon 6     PreStart= 200301 PreEnd= 200504 PreDura= 12 Pre_Fmw= 5.15 Pre_Smw= 6.9  Pre_Emw= 6.9  PostStart= 200601 PostEnd= 200602 PostDura= 1  Post_Fmw= 5.15 Post_Smw= 7.5  Post_Emw= 7.5  Size_Fmwch= 0   Size_Smwch= 0.6  Size_Emwch= 0.6 Oregon     6 Pre_mw= 6.9  Size_mwch= 0.6  PreStart= 200301 PreEnd= 200504 PostStart= 200601 PostEnd= 200602 PreDura= 12 PostDura= 1 
# Washington 1 PreStart= 199001 PreEnd= 199304 PreDura= 16 Pre_Fmw= 4.25 Pre_Smw= 4.25 Pre_Emw= 4.25 PostStart= 199401 PostEnd= 199804 PostDura= 20 Post_Fmw= 4.25 Post_Smw= 4.9  Post_Emw= 4.9  Size_Fmwch= 0   Size_Smwch= 0.65 Size_Emwch= 0.65Washington 1 Pre_mw= 4.25 Size_mwch= 0.65 PreStart= 199001 PreEnd= 199304 PostStart= 199401 PostEnd= 199804 PreDura= 16 PostDura= 20 
# Washington 2 PreStart= 199401 PreEnd= 199804 PreDura= 20 Pre_Fmw= 5.15 Pre_Smw= 4.9  Pre_Emw= 5.15 PostStart= 199901 PostEnd= 199904 PostDura= 4  Post_Fmw= 5.15 Post_Smw= 5.7  Post_Emw= 5.7  Size_Fmwch= 0   Size_Smwch= 0.8  Size_Emwch= 0.55Washington 2 Pre_mw= 4.9  Size_mwch= 0.8  PreStart= 199401 PreEnd= 199804 PostStart= 199901 PostEnd= 199904 PreDura= 20 PostDura= 4 
# Washington 3 PreStart= 199901 PreEnd= 199904 PreDura= 4  Pre_Fmw= 5.15 Pre_Smw= 5.7  Pre_Emw= 5.7  PostStart= 200001 PostEnd= 200004 PostDura= 4  Post_Fmw= 5.15 Post_Smw= 6.5  Post_Emw= 6.5  Size_Fmwch= 0   Size_Smwch= 0.8  Size_Emwch= 0.8 Washington 3 Pre_mw= 5.7  Size_mwch= 0.8  PreStart= 199901 PreEnd= 199904 PostStart= 200001 PostEnd= 200004 PreDura= 4 PostDura= 4 
# Washington 4 PreStart= 200001 PreEnd= 200004 PreDura= 4  Pre_Fmw= 5.15 Pre_Smw= 6.5  Pre_Emw= 6.5  PostStart= 200101 PostEnd= 200104 PostDura= 4  Post_Fmw= 5.15 Post_Smw= 6.72 Post_Emw= 6.72 Size_Fmwch= 0   Size_Smwch= 0.22 Size_Emwch= 0.22Washington 4 Pre_mw= 6.5  Size_mwch= 0.22 PreStart= 200001 PreEnd= 200004 PostStart= 200101 PostEnd= 200104 PreDura= 4 PostDura= 4 
# Washington 5 PreStart= 200101 PreEnd= 200104 PreDura= 4  Pre_Fmw= 5.15 Pre_Smw= 6.72 Pre_Emw= 6.72 PostStart= 200201 PostEnd= 200204 PostDura= 4  Post_Fmw= 5.15 Post_Smw= 6.9  Post_Emw= 6.9  Size_Fmwch= 0   Size_Smwch= 0.18 Size_Emwch= 0.18Washington 5 Pre_mw= 6.72 Size_mwch= 0.18 PreStart= 200101 PreEnd= 200104 PostStart= 200201 PostEnd= 200204 PreDura= 4 PostDura= 4 
# Washington 6 PreStart= 200201 PreEnd= 200204 PreDura= 4  Pre_Fmw= 5.15 Pre_Smw= 6.9  Pre_Emw= 6.9  PostStart= 200301 PostEnd= 200304 PostDura= 4  Post_Fmw= 5.15 Post_Smw= 7.01 Post_Emw= 7.01 Size_Fmwch= 0   Size_Smwch= 0.11 Size_Emwch= 0.11Washington 6 Pre_mw= 6.9  Size_mwch= 0.11 PreStart= 200201 PreEnd= 200204 PostStart= 200301 PostEnd= 200304 PreDura= 4 PostDura= 4 
# Washington 7 PreStart= 200301 PreEnd= 200304 PreDura= 4  Pre_Fmw= 5.15 Pre_Smw= 7.01 Pre_Emw= 7.01 PostStart= 200401 PostEnd= 200404 PostDura= 4  Post_Fmw= 5.15 Post_Smw= 7.16 Post_Emw= 7.16 Size_Fmwch= 0   Size_Smwch= 0.15 Size_Emwch= 0.15Washington 7 Pre_mw= 7.01 Size_mwch= 0.15 PreStart= 200301 PreEnd= 200304 PostStart= 200401 PostEnd= 200404 PreDura= 4 PostDura= 4 
# Washington 8 PreStart= 200401 PreEnd= 200404 PreDura= 4  Pre_Fmw= 5.15 Pre_Smw= 7.16 Pre_Emw= 7.16 PostStart= 200501 PostEnd= 200504 PostDura= 4  Post_Fmw= 5.15 Post_Smw= 7.35 Post_Emw= 7.35 Size_Fmwch= 0   Size_Smwch= 0.19 Size_Emwch= 0.19Washington 8 Pre_mw= 7.16 Size_mwch= 0.19 PreStart= 200401 PreEnd= 200404 PostStart= 200501 PostEnd= 200504 PreDura= 4 PostDura= 4 
# Washington 9 PreStart= 200501 PreEnd= 200504 PreDura= 4  Pre_Fmw= 5.15 Pre_Smw= 7.35 Pre_Emw= 7.35 PostStart= 200601 PostEnd= 200602 PostDura= 1  Post_Fmw= 5.15 Post_Smw= 7.63 Post_Emw= 7.63 Size_Fmwch= 0   Size_Smwch= 0.28 Size_Emwch= 0.28Washington 9 Pre_mw= 7.35 Size_mwch= 0.28 PreStart= 200501 PreEnd= 200504 PostStart= 200601 PostEnd= 200602 PreDura= 4 PostDura= 1 
# Wisconsin 1  PreStart= 199001 PreEnd= 200501 PreDura= 61 Pre_Fmw= 5.15 Pre_Smw= 3.65 Pre_Emw= 5.15 PostStart= 200502 PostEnd= 200602 PostDura= 4  Post_Fmw= 5.15 Post_Smw= 5.7  Post_Emw= 5.7  Size_Fmwch= 0   Size_Smwch= 2.05 Size_Emwch= 0.55Wisconsin  1 Pre_mw= 3.65 Size_mwch= 2.05 PreStart= 199001 PreEnd= 200501 PostStart= 200502 PostEnd= 200602 PreDura= 61 PostDura= 4 

#---------------------------------------------------------#
# 2. Check outliers for the case studies with PreDura>=20 #
#---------------------------------------------------------#
case=1
while(case<=9){
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="Florida";    Pre_Smw=0   ; Size_Smwch=6.15; Pre_Emw=5.15; Size_Emwch=1.00; PreStart=199001; PreEnd=200501; PostStart=200502; PostEnd=200504; PreDura=61; PostDura=3 } 
if(case== 3){trt_name="Illinois";   Pre_Smw=0   ; Size_Smwch=5.5 ; Pre_Emw=5.15; Size_Emwch=0.35; PreStart=199001; PreEnd=200304; PostStart=200401; PostEnd=200404; PreDura=56; PostDura=4 } 
if(case== 4){trt_name="Maryland";   Pre_Smw=0   ; Size_Smwch=6.15; Pre_Emw=5.15; Size_Emwch=1.00; PreStart=199001; PreEnd=200504; PostStart=200601; PostEnd=200602; PreDura=64; PostDura=1 } 
if(case== 5){trt_name="Minnesota";  Pre_Smw=4.25; Size_Smwch=1.9 ; Pre_Emw=5.15; Size_Emwch=1.00; PreStart=199101; PreEnd=200502; PostStart=200503; PostEnd=200602; PreDura=58; PostDura=3 } 
if(case== 6){trt_name="New York";   Pre_Smw=0   ; Size_Smwch=6   ; Pre_Emw=5.15; Size_Emwch=0.85; PreStart=199001; PreEnd=200404; PostStart=200501; PostEnd=200504; PreDura=60; PostDura=4 } 
if(case== 7){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== 8){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== 9){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 } 

# 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    

outlier=100  # Threshold for defining an outlier
cr=rainbow(length(tau),s=1,v=1,start=0,end=0.9) # Rainbow colors for the quantile sequences
par(mfrow=c(8,4),mar=c(2,2,2,2))
S=d-a+1
Oindx0=0
i=which(Sname==trt_name)
wm=matrix(0,S,1)
wq=matrix(0,S,length(tau))
wQ=matrix(0,S,1)
s=1
while(s<=S){
t=a+s-1
w=na.omit(y1[which(sindx==i & time==Time[t])])
wm[s]=mean(w)                        # mean
wq[s,]=quantile(w,tau)      # quantile 
wQ[s]=max(w)
s=s+1}
tt=which.max(wQ)
if(tt==1){DwQ=wQ[tt]-wQ[tt+1]}
if(tt==S){DwQ=wQ[tt]-wQ[tt-1]}
if(tt>1 & tt<S){DwQ=max((wQ[tt]-wQ[tt-1]),(wQ[tt]-wQ[tt+1]))}
if(DwQ>outlier){Oindx0=c(Oindx0,i)             # Definition of outliers
print("Outliers detected")
cat(case,Sname[i],Sindx[i],Time[tt],DwQ,"\n")}
case=case+1}

### Section 2.2.3 of the Appendix: (Detected outliers)
#2 Florida 9 199403 146.8748 
#3 Illinois 12 199302 168.8267 
#4 Maryland 19 199001 459.6809 
#5 Minnesota 22 199402 743.9656 
#6 New York 31 199103 444.1349 

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


##-------------------------------------------#
## Plot: States' wage distribution evolution #
##-------------------------------------------#
#y1max=400
#cr=rainbow(T,s=1,v=1,start=0,end=0.9) # Color: rainbow
#
#{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/distevol_wage.pdf",sep=""))
#par(mfrow=c(8,4),mar=c(2,2,2,2))
#i=1
#while(i<=N){
#if(i==1){Sys.sleep(1)}
#if(Ni[i]>minN){
#w=y1[which(sindx==i & time==Time[1])]
#mt=as.character(Sname[i])
#res=density(na.omit(w),from=0,to=y1max)
#x_axis=res[[1]]
#f1=res[[2]]
#plot(x_axis,f1,type="l",,lty=1,lwd=0.5,col=cr[1],ylim=c(0,0.05),xlab=" ",ylab=" ",main=mt)
#t=2
#while(t<=T){
#w=y1[which(sindx==i & time==Time[t])]
#res=density(na.omit(w),from=0,to=y1max)
#ft=res[[2]]
#lines(x_axis,ft,col=cr[t],lty=1,lwd=0.5)
#t=t+1}
#}
#i=i+1}
#par(mfrow=c(1,1))
#dev.off()}
#
##--------------------------------------------------#
## Plot: States' employement distribution evolution #
##--------------------------------------------------#
#y2max=300
#cr=rainbow(T,s=1,v=1,start=0,end=0.9) # Color: rainbow
#{pdf(paste("E:/Research/SyntheticControl/Empirical/TempPDFs","/distevol_emp.pdf",sep=""))
#par(mfrow=c(8,4),mar=c(2,2,2,2))
#i=1
#while(i<=N){
#if(i==1){Sys.sleep(1)}
#if(Ni[i]>minN){
#w=y2[which(sindx==i & time==Time[1])]
#mt=as.character(Sname[i])
#res=density(na.omit(w),from=0,to=y2max)
#x_axis=res[[1]]
#f1=res[[2]]
#plot(x_axis,f1,type="l",lty=1,lwd=0.5,col=cr[1],xlab=" ",ylab=" ",main=mt)
#t=2
#while(t<=T){
#w=y2[which(sindx==i & time==Time[t])]
#res=density(na.omit(w),from=0,to=y2max)
#ft=res[[2]]
#lines(x_axis,ft,col=cr[t],lty=1,lwd=0.5)
#t=t+1}
#}
#i=i+1}
#par(mfrow=c(1,1))
#dev.off()}
