Ratio Root Mean Square Error

Import and Cleaning Dataset

PASW_all <- read_excel("PASW.xlsx",
                   col_types = c("text", "date", "numeric", "numeric", "numeric"))

## Exclude Dictrics of Coloumnia and Puerto Rico
PASW_all<-filter(PASW_all,STATE!="DC",STATE!="PR",DATE>="1998-01-01")
#ggplot(PASW)+geom_line(aes(y=PASW,x=DATE,col=STATE))


#COump Dataframe
dumy_date=seq.Date(from = as.Date("1998-01-01"),to = as.Date("2008-01-01"),by = "quarter")
Fix_df<-data_frame(STATE=rep(unique(countypop$abbr),length(dumy_date)))
Fix_df<-arrange(Fix_df,STATE)
Fix_df$DATE<-rep(dumy_date,length(unique(countypop$abbr)))

PASW<-left_join(Fix_df,PASW_all,by=c("STATE","DATE"))

##SELCET
PASW<-PASW%>%select(STATE,PASW,DATE)
PASW<-PASW %>% group_by(STATE) %>%  filter(!any(is.na(PASW)))
PASW$DATE<-as.Date(PASW$DATE)

ST<-unique(PASW$STATE)
n<-length(ST)
date<-seq.Date(from = as.Date("2005-01-01"),as.Date("2008-01-01"),by = "quarter") #Month after treatement
#date<-AER$DATE[69:97]
TD<-as.character(date)

RESULT<-data.frame(State=rep(ST,length(date)))
SE<-data.frame(SE=rep(0,length(date)))

Calculate All State MSE by looping a Sectors

#Looping
for(j in 1:length(date)){
  Clean_PASW<-filter(PASW,DATE <= TD[j])
  for(i in 1:n){
    Clean_PASW$Treat<-ifelse(Clean_PASW$DATE>="2005-01-01" & 
                              Clean_PASW$DATE<="2022-01-01" & 
                              Clean_PASW$STATE==ST[i],1,0)
    Clean_PASW<-as.data.frame(Clean_PASW)
    Clean_PASW$STATE<-factor(Clean_PASW$STATE)
    
    # SDID PAN
    matx<-panel.matrices(Clean_PASW,unit = 1,time = 3,outcome = 2,treatment = 4)
    sydid<-synthdid_estimate(matx$Y, matx$N0, matx$T0)
    #plot(sydid)
    RESULT$Estimate[i+((j-1)*n)]=summary(sydid)$estimate
    RESULT$Group[i+((j-1)*n)]=TD[j]
    #RESULT$SE[i+((j-1)*n)]=sqrt(vcov(sydid, method='placebo',replications = 5))
    if(ST[i]=="LA"){
      SE$Group[j]=TD[j]
      SE$SE[j]=sqrt(vcov(sydid, method='placebo',replications = 30))
      SE$estimate[j]=summary(sydid)$estimate
      #SE$SC[j]<-summary(sc_estimate(matx$Y, matx$N0, matx$T0))$estimate
    }
    
    weights = attr(sydid, "weights")
    Y = matx$Y
    T0 = matx$T0; T1 = ncol(Y) - T0
    N0 = matx$N0; N1 = nrow(Y) - N0
    omega.synth = c(weights$omega, rep(0, N1))
    
    omega.synth = c(weights$omega, rep(0, N1))
    omega.target = c(rep(0, N0), rep(1 / N1, N1))
    lambda.synth = c(weights$lambda, rep(0, T1))
    lambda.target = c(rep(0, T0), rep(1 / T1, T1))
    #over = attr(sydid, 'overlay')
    over=1
    
    intercept.offset = over * c((omega.target - omega.synth) %*% Y %*% lambda.synth)
    obs.trajectory = as.numeric(omega.target %*% Y)
    syn.trajectory = as.numeric(omega.synth %*% Y) + intercept.offset
    syn.trajectory_1 = as.numeric(omega.synth %*% Y)
    #plot(sydid)
    
    ## Residual
    Residu<-syn.trajectory-obs.trajectory
    Residu1<-syn.trajectory_1-obs.trajectory
    
    ## Root MSE Before Treatment
    residu_pre<-sqrt(mean(Residu[1:23]^2)) 
    residu_pre1<-sqrt(mean(Residu1[1:23]^2))
    leng=length(Residu)
    
    ## Root MSE After Treatment
    residu_past<-sqrt(mean(Residu[24:leng]^2))
    residu_past1<-sqrt(mean(Residu1[24:leng]^2))
    
    residu_past_point<-sqrt(mean(Residu[leng]^2))
    
    ## Ratio of Root MSE After and Before Treatment
    RESULT$Ratio_RMSE[i+((j-1)*n)]<-residu_past/residu_pre
    RESULT$Ratio_RMSE1[i+((j-1)*n)]<-residu_past1/residu_pre1
    RESULT$Ratio_RMSE_Point[i+((j-1)*n)]<-residu_past_point/residu_pre
    
    ## Save Output Value
    RESULT$MSE_pre[i+((j-1)*n)]<-residu_pre^2
    RESULT$TRUE_VAL[i+((j-1)*n)]<-obs.trajectory[leng]
    RESULT$CONTROL[i+((j-1)*n)]<-syn.trajectory[leng]
    RESULT$CONTROL_I[i+((j-1)*n)]<-syn.trajectory_1[leng]
    
  }
}

Calculate Ranking

## Ranking by looping for each period of time
for(j in 1:length(date)){
  RESULT$Rank[(((j-1)*n)+1):(j*n)]<-rank(-RESULT$Ratio_RMSE[(((j-1)*n)+1):(j*n)])
  RESULT$Rank_1[(((j-1)*n)+1):(j*n)]<-rank(-RESULT$Ratio_RMSE1[(((j-1)*n)+1):(j*n)])
  RESULT$Rank_Point[(((j-1)*n)+1):(j*n)]<-rank(-RESULT$Ratio_RMSE_Point[(((j-1)*n)+1):(j*n)])
}

RESULT$p_Value<-RESULT$Rank/n
RESULT$p_Value1<-RESULT$Rank_1/n
RESULT$Loss_Ratio<-RESULT$Estimate/RESULT$TRUE_VAL
RESULT$Predict<-RESULT$CONTROL+RESULT$Estimate
LA_PASW<-filter(RESULT,State=="LA")

for(j in 1:length(date)){
  bMSE<-filter(RESULT,MSE_pre<2*LA_PASW$MSE_pre[j],Group==TD[j])%>%nrow()
  print(bMSE)
  LA_PASW$Actual_p[j]<-LA_PASW$Rank[j]/bMSE
  LA_PASW$Actual_p_Poin[j]<-LA_PASW$Rank_Point[j]/bMSE
  LA_PASW$SDID_TE[j]<-LA_PASW$Estimate[j]*j-sum(LA_PASW$SDID_TE[1:j-1])
}

## [1] 44
## [1] 44
## [1] 45
## [1] 43
## [1] 43
## [1] 43
## [1] 43
## [1] 45
## [1] 44
## [1] 44
## [1] 44
## [1] 44
## [1] 44

Ratio Root Mean Square Error

2023-06-07

Import and Cleaning Dataset

Calculate All State MSE by looping a Sectors

Calculate Ranking