## Plot QR Evaluation Results #### load("v9_forecasts/Evaluation.Rda") load("v9_forecasts/GAM-Grid.Rda") ## Plots for paper: #### # Deterministic: MAE plot_out <- ggplot(data=Evaluation[kfold=="All_cv" & Score %in% c("MAE_gam"),],aes(y=Value,x=Node,group=MODEL))+ geom_point(aes(col=MODEL,shape=MODEL)) + xlab("GSP Group") + ylab("Mean Absolute Error [MWh]") + theme_bw() + theme(legend.position="top",legend.title = element_blank(),text=element_text(family="serif")) plot_out ggsave("paper_plots/MAE_allCV.pdf", plot_out, width = 6, height = 4,device = "pdf") # Probabilistic: Pinball plot_out <- ggplot(data=Evaluation[kfold=="Test" & Tails=="QR" & Score %in% c("Pinball") & Quantile>=0.1 & Quantile<=0.9, .(Value=mean(Value)),by=c("Node","MODEL")], aes(y=Value,x=Node,group=MODEL))+ geom_point(aes(col=MODEL,shape=MODEL)) + xlab("GSP Group") + ylab("Pinball Score") + theme_bw() + theme(legend.position="top",legend.title = element_blank(),text=element_text(family="serif")) plot_out ggsave("paper_plots/Pinball_q10_q90_allCV.pdf", plot_out, width = 6, height = 4,device = "pdf") # Table PB_score <- Evaluation[kfold%in%c("All_cv","Test") & Tails=="QR" & Score %in% c("Pinball") & Quantile>=0.1 & Quantile<=0.9, .(`Pinball Score`=mean(Value)),by=c("kfold","MODEL")] PB_score <- dcast(PB_score,formula = MODEL~kfold,value.var = "Pinball Score") PB_score_2 <- Evaluation[kfold%in%c("Test") & Node %in% c("C","H","P") & Tails=="QR" & Score %in% c("Pinball") & Quantile>=0.1 & Quantile<=0.9, .(`Pinball Score`=mean(Value)),by=c("kfold","MODEL","Node")] PB_score_2 <- dcast(PB_score_2,formula = MODEL~kfold+Node,value.var = "Pinball Score") print(xtable::xtable(merge(PB_score,PB_score_2,by="MODEL")[5:1,],digits=3), include.rownames=FALSE) 100*(1-PB_score$Test/PB_score$Test[3]) 100*(1-PB_score_2[2,-1]/PB_score_2[3,-1]) rm(PB_score,PB_score_2) # Reliability plot_out <- ggplot(data=Evaluation[kfold=="Test" & Score %in% c("Reliability") & MODEL %in% c("GAM-Point","Vanilla-Point"),],aes(y=Value,x=Quantile,group=Node))+ geom_line(aes(col=Node)) + geom_point(aes(col=Node)) + facet_wrap( ~ MODEL,scales = "fixed") + theme_bw() + theme(legend.position="none",text=element_text(family="serif"), strip.background =element_rect(fill="white")) + ylim(0,1) + ylab("Observed frequency") + coord_fixed() + xlab("Forecast probability") plot_out ggsave("paper_plots/ReliabilityDiag_Test.pdf", plot_out, width = 6, height = 4,device = "pdf") plot_out <- ggplot(data=Evaluation[kfold=="All_cv" & Score %in% c("Reliability") & MODEL %in% c("GAM-Point","Vanilla-Point"),],aes(y=Value,x=Quantile,group=Node))+ geom_line(aes(col=Node)) + geom_point(aes(col=Node)) + facet_wrap( ~ MODEL,scales = "fixed") + theme_bw() + theme(legend.position="none",text=element_text(family="serif"), strip.background =element_rect(fill="white")) + ylim(0,1) + ylab("Observed frequency") + coord_fixed() + xlab("Forecast probability") plot_out ggsave("paper_plots/ReliabilityDiag_CV.pdf", plot_out, width = 6, height = 4,device = "pdf") # Interval width temp_eval <- Evaluation[kfold=="All_cv" & Score %in% c("Interval Width") & Quantile>=0.99 & MODEL =="GAM-Point" & Node %in%c("C","H","P"),] temp_eval[Node=="C",Node:=paste0("GSP Group ",Node," (Greater London)")] temp_eval[Node=="H",Node:=paste0("GSP Group ",Node," (South, High Embedded Solar)")] temp_eval[Node=="P",Node:=paste0("GSP Group ",Node," (North, High Embedded Wind)")] plot_out <- ggplot(data=temp_eval, aes(y=Value,x=100*Quantile,col=Tails,shape=Tails))+ geom_line() + geom_point() + facet_wrap( ~ Node,scales = "free_y",ncol = 1) + scale_x_continuous(breaks=seq(99,100,by=.2))+ theme_bw() + theme(legend.position="right",text=element_text(family="serif"), strip.background =element_rect(fill="white")) + xlab("Coverage [%]") + ylab("Average Interval Width") plot_out ggsave("paper_plots/SharpnessDiag_CV.pdf", plot_out, width = 6, height = 4,device = "pdf") rm(temp_eval) ## Worm Plots par(family = "serif",mar=c(3.1,3.1,4.1,.1),mgp=c(2,1,0)) new<-T; for(N in allNs){ WormPlot(U_data = NodeData_temp[[paste0(N,"_u")]][which(!NodeData_temp[[N]]$BadData)], cov_max_lag = 1,add=!new); new=F} dev.print(png, 'paper_plots/Worm_cGPD.png',width=6,height=4,units="in",res=600) new<-T;for(N in allNs){ WormPlot(U_data = NodeData_temp[[paste0(N,"_u_sGPD")]][which(!NodeData_temp[[N]]$BadData)], cov_max_lag = 1,add=!new,col=3); new=F} dev.print(png, 'paper_plots/Worm_sGPD.png',width=6,height=4,units="in",res=600) new<-T;for(N in allNs){ WormPlot(U_data = NodeData_temp[[paste0(N,"_u_QR")]][which(!NodeData_temp[[N]]$BadData)], cov_max_lag = 1,add=!new,col=4); new=F} dev.print(png, 'paper_plots/Worm_QR.png',width=6,height=4,units="in",res=600) ## Plots for Sup Material: #### # Det (GAM) plot_out <- ggplot(data=Evaluation[kfold!="All_cv" & Score %in% c("MAE_gam"),],aes(y=Value,x=Node,group=MODEL))+ geom_point(aes(col=MODEL)) + facet_wrap( ~ kfold,scales = "free") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white")) + ylab("MAE [MWh]") + xlab("GSP Group") plot_out ggsave("paper_sup_plots/MAE_CV.pdf", plot_out, width = 6, height = 4,device = "pdf") plot_out <- ggplot(data=Evaluation[kfold!="All_cv" & Score %in% c("sMAE_gam"),],aes(y=Value,x=Node,group=MODEL))+ geom_point(aes(col=MODEL)) + facet_wrap( ~ kfold,scales = "free") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white")) + ylab("MAE [standardised]") + xlab("GSP Group") plot_out ggsave("paper_sup_plots/MAE_Stand_CV.pdf", plot_out, width = 6, height = 4,device = "pdf") # Pinball plot_out <- ggplot(data=Evaluation[kfold=="All_cv" & Score %in% c("Pinball") ,],aes(y=Value,x=Quantile,group=MODEL))+ geom_line(aes(col=MODEL)) + facet_wrap( ~ Node,scales = "fixed") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white")) + ylab("Pinbal Loss [standardised]") + xlab("Quantile") plot_out ggsave("paper_sup_plots/Pinball_CV.pdf", plot_out, width = 6, height = 8,device = "pdf") plot_out <- ggplot(data=Evaluation[kfold=="Test" & Score %in% c("Pinball"),],aes(y=Value,x=Quantile,group=MODEL))+ geom_line(aes(col=MODEL)) + facet_wrap( ~ Node,scales = "fixed") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white")) + ylab("Pinbal Loss [standardised]") + xlab("Quantile") plot_out ggsave("paper_sup_plots/Pinball_Test.pdf", plot_out, width = 6, height = 8,device = "pdf") ## Pinball skill-score Eval_Wide <- dcast.data.table(Evaluation[Score=="Pinball" & Quantile>=0.1 & Quantile<=0.9 & Tails=="QR",], Quantile+Node+kfold~MODEL,value.var = "Value") # Grid vs Point Eval_Wide[,Skill:=100*(1-`GAM-Grid`/`GAM-Point`)] plot_out <- ggplot(data=Eval_Wide[kfold%in%c(1:3,"Test"),],aes(y=Skill,x=Quantile,group=Node,col=Node))+ geom_line() + geom_point() + facet_wrap( ~ kfold,scales = "fixed") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white")) + ylab("Pinbal Loss Skill Score [%]") + xlab("Quantile") + labs(color="GSP Group") plot_out ggsave("paper_sup_plots/Pinball_Skill_GridvsPo.pdf", plot_out, width = 6, height = 6,device = "pdf") # GAM vs Vanilla Eval_Wide[,Skill:=100*(1-`GAM-Point`/`Vanilla-Point`)] plot_out <- ggplot(data=Eval_Wide[kfold%in%c(1:3,"Test"),],aes(y=Skill,x=Quantile,group=Node,col=Node))+ geom_line() + geom_point() + facet_wrap( ~ kfold,scales = "fixed") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white")) + ylab("Pinbal Loss Skill Score [%]") + xlab("Quantile") + labs(color="GSP Group") plot_out ggsave("paper_sup_plots/Pinball_Skill_GAMvsV.pdf", plot_out, width = 6, height = 8,device = "pdf") # Reliability plot_out <- ggplot(data=Evaluation[kfold=="All_cv" & Score %in% c("Reliability"),],aes(y=Quantile-Value,x=Quantile,group=Node))+ geom_line(aes(col=Node)) + geom_point(aes(col=Node)) + facet_wrap( ~ MODEL,scales = "fixed") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white")) + ylim(-0.01,0.01) + ylab("Quantile Bias") + labs(color="GSP Group") plot_out ggsave("paper_sup_plots/Rel_qbias_cv.pdf", plot_out, width = 6, height = 8,device = "pdf") plot_out <- ggplot(data=Evaluation[kfold=="Test" & Score %in% c("Reliability"),],aes(y=Quantile-Value,x=Quantile,group=Node))+ geom_line(aes(col=Node)) + geom_point(aes(col=Node)) + facet_wrap( ~ MODEL,scales = "fixed") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white")) + ylim(-0.3,0.3) + ylab("Quantile Bias")+ labs(color="GSP Group") plot_out ggsave("paper_sup_plots/Rel_qbias_Test.pdf", plot_out, width = 6, height = 8,device = "pdf") plot_out <- ggplot(data=Evaluation[kfold=="All_cv" & Score %in% c("Reliability"),],aes(y=Value,x=Quantile,group=Node))+ geom_line(aes(col=Node)) + geom_point(aes(col=Node)) + facet_wrap( ~ MODEL,scales = "fixed") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white")) + ylim(0,1) + ylab("Empirical") + labs(color="GSP Group") plot_out ggsave("paper_sup_plots/Rel_cv.pdf", plot_out, width = 6, height = 8,device = "pdf") plot_out <- ggplot(data=Evaluation[kfold=="Test" & Score %in% c("Reliability"),],aes(y=Value,x=Quantile,group=Node))+ geom_line(aes(col=Node)) + geom_point(aes(col=Node)) + facet_wrap( ~ MODEL,scales = "fixed") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white"))+ ylim(0,1) + ylab("Empirical")+ labs(color="GSP Group") plot_out ggsave("paper_sup_plots/Rel_Test.pdf", plot_out, width = 6, height = 8,device = "pdf") ## Skill score - Point vs Grid by Quantile ## Sharpness plot_out <- ggplot(data=Evaluation[kfold=="All_cv" & Score %in% c("Interval Width") & Quantile>0.9,], aes(y=Value,x=Quantile,col=Tails,shape=MODEL))+ geom_line() + geom_point() + facet_wrap( ~ Node,scales = "fixed") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white"), legend.box="vertical") + xlab("Coverage") + ylab("Average Interval Width") plot_out ggsave("paper_sup_plots/Sharpness_90_cv.pdf", plot_out, width = 6, height = 8,device = "pdf") plot_out <- ggplot(data=Evaluation[kfold=="All_cv" & Score %in% c("Interval Width") & Quantile>=0.99 & MODEL =="GAM-Grid",], aes(y=Value,x=Quantile,col=Tails,shape=MODEL))+ geom_line() + geom_point() + facet_wrap( ~ Node,scales = "fixed") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white"), legend.box="vertical") + xlab("Coverage") + ylab("Average Interval Width") plot_out ggsave("paper_sup_plots/Sharpness_99_cv.pdf", plot_out, width = 6, height = 8,device = "pdf") plot_out <- ggplot(data=Evaluation[kfold=="Test" & Score %in% c("Interval Width") & Quantile>0.9,], aes(y=Value,x=Quantile,col=Tails,shape=MODEL))+ geom_line() + geom_point() + facet_wrap( ~ Node,scales = "fixed") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white"), legend.box="vertical") + xlab("Coverage") + ylab("Average Interval Width") plot_out ggsave("paper_sup_plots/Sharpness_90_Test.pdf", plot_out, width = 6, height = 8,device = "pdf") plot_out <- ggplot(data=Evaluation[kfold=="Test" & Score %in% c("Interval Width") & Quantile>=0.99 & MODEL =="GAM-Grid",], aes(y=Value,x=Quantile,col=Tails,shape=MODEL))+ geom_line() + geom_point() + facet_wrap( ~ Node,scales = "fixed") + theme_bw() + theme(legend.position="bottom",text=element_text(family="serif"), strip.background =element_rect(fill="white"), legend.box="vertical") + xlab("Coverage") + ylab("Average Interval Width") plot_out ggsave("paper_sup_plots/Sharpness_99_Test.pdf", plot_out, width = 6, height = 8,device = "pdf") ## Reserve Volume Change #### Reserve_Diff <- data.table(); for(N in allNs){ temp_MQR <- NodeData_temp[[paste0(N,"_qr")]][,NodeData_temp[[paste0(N,"_evtails")]]$mqr_is,with=F] temp_MQR_s <- copy(temp_MQR) for(prob in c(0.98,0.985,0.99,0.995,0.9975,0.999,0.9995,0.9999)){ tail_from <-NodeData_temp[[paste0(N,"_evtails")]]$tails_from temp_MQR[[paste0("q",prob*100)]] <- temp_MQR[[paste0("q",100-tail_from)]]+ qgpd(p=1-(prob-1+tail_from/100)/(tail_from/100), shape=NodeData_temp[[N]]$gpd_shape_r, scale=NodeData_temp[[N]]$gpd_scale_r) temp_MQR[[paste0("q",signif(100-prob*100,3))]] <- temp_MQR[[paste0("q",tail_from)]]- qgpd(p=1-(prob-1+tail_from/100)/(tail_from/100), shape=NodeData_temp[[N]]$gpd_shape_l, scale=NodeData_temp[[N]]$gpd_scale_l) temp_MQR_s[[paste0("q",prob*100)]] <- temp_MQR_s[[paste0("q",100-tail_from)]]+ qgpd(p=1-(prob-1+tail_from/100)/(tail_from/100), shape=NodeData_temp[[N]]$gpd_shape_r0, scale=NodeData_temp[[N]]$gpd_scale_r0) temp_MQR_s[[paste0("q",signif(100-prob*100,3))]] <- temp_MQR_s[[paste0("q",tail_from)]]- qgpd(p=1-(prob-1+tail_from/100)/(tail_from/100), shape=NodeData_temp[[N]]$gpd_shape_l0, scale=NodeData_temp[[N]]$gpd_scale_l0) } temp_MQR <- temp_MQR[,order(as.numeric(gsub("q","",names(temp_MQR)))),with=F] temp_MQR_s <- temp_MQR_s[,order(as.numeric(gsub("q","",names(temp_MQR_s)))),with=F] NaiveQs <- quantile((NodeData_temp[[N]]$node_n-NodeData_temp[[paste0(N,"_qr")]]$q50)[which(!NodeData_temp[[N]]$BadData)], as.numeric(gsub("q","",names(temp_MQR_s)))/100,na.rm = T) NaiveQs_long <- t(matrix(rep(temp_MQR[,q50],ncol(temp_MQR)),nrow=ncol(temp_MQR),byrow = T) + NaiveQs) dt <- data.table(Diff=colMeans((temp_MQR-temp_MQR_s)[which(!NodeData_temp[[N]]$BadData),],na.rm = T)*NodeData_temp[[paste0(N,"_norm")]][1], Diff_pc=colMeans((temp_MQR-temp_MQR_s)[which(!NodeData_temp[[N]]$BadData),],na.rm = T)/colMeans((temp_MQR_s-temp_MQR_s$q50)[which(!NodeData_temp[[N]]$BadData),],na.rm = T), Diff_c_gt_s_pc=colMeans((temp_MQR>temp_MQR_s)[which(!NodeData_temp[[N]]$BadData),],na.rm = T), # Diff_naive=colMeans((temp_MQR-NaiveQs_long)[which(!NodeData_temp[[N]]$BadData),],na.rm = T)*NodeData_temp[[paste0(N,"_norm")]][1], Diff_naive_pc=colMeans((temp_MQR-NaiveQs_long)[which(!NodeData_temp[[N]]$BadData),],na.rm = T)/NaiveQs, Diff_c_gt_naive_pc=colMeans((temp_MQR>NaiveQs_long)[which(!NodeData_temp[[N]]$BadData),],na.rm = T), # Quantile=as.numeric(gsub("q","",names(temp_MQR))), `GSP Group`=N) dt[Quantile>0.5,Diff_c_gt_naive_pc:=1-Diff_c_gt_naive_pc] dt[Quantile>0.5,Diff_c_gt_s_pc:=1-Diff_c_gt_s_pc] Reserve_Diff <- rbind(Reserve_Diff,dt[Diff!=0,]) if(N=="H"){ ## Plot ## Volume Plot issue <- NodeData_temp[[N]][,.(BadData=any(BadData)), by=issueTime][BadData==F,issueTime][1703] index <- which(NodeData_temp[[N]]$issueTime==issue) # Fan par(family = "serif",mar=c(3.1,3.1,.1,.1),mgp=c(2,1,0)) xx <- seq(22,47.5,by=0.5) {plot(temp_MQR[index,2:(ncol(temp_MQR)-1),with=F],targetTimes = xx, xlab="Lead-time [h]", ylab="Net-load [standardised]") legend("topleft",c("cGPD 99.9% Prediction Inverval", "Reserve Volume Avoided","Additional Reserve Volume"), pch = 15,col=c("cyan","green","red"),bty = "n") # Volume Saved # Naive <- temp_MQR[index,q50]+NaiveQs[1] Naive <- temp_MQR_s[index,q0.05] Vol_saved <- temp_MQR[index,q0.05] Vol_saved[Vol_saved=Naive] <- Naive[Vol_needed>=Naive] polygon(x=c(xx,rev(seq(22,47.5,by=0.5))), y = c(Vol_needed,rev(Naive)),col = rgb(1,0,0,alpha = 1))} dev.print(png, paste0("paper_plots/Reserve_s_vs_cGPD_",N,".png"),width=6,height=4,units="in",res=600) } rm(dt,NaiveQs) } ## Static vs Conditional ggplot(Reserve_Diff[Quantile<1,],aes(x=Quantile,y=Diff_pc*100,group=`GSP Group`,col=`GSP Group`))+ geom_line() + geom_point() +theme_bw() + ylab("Change in Reserve [%]") ggplot(Reserve_Diff[Quantile>99,],aes(x=Quantile,y=Diff_pc*100,group=`GSP Group`,col=`GSP Group`))+ geom_line() + geom_point() +theme_bw() + ylab("Change in Reserve [%]") ## Naive vs Conditional ggplot(Reserve_Diff[Quantile<1,],aes(x=Quantile,y=Diff_naive_pc*100,group=`GSP Group`,col=`GSP Group`))+ geom_line() + geom_point() +theme_bw() + ylab("Change in Reserve [%]") ggplot(Reserve_Diff[Quantile>99,],aes(x=Quantile,y=Diff_naive_pc*100,group=`GSP Group`,col=`GSP Group`))+ geom_line() + geom_point() +theme_bw() + ylab("Change in Reserve [%]") ## Table for paper require(xtable) reserve_xtable <- Reserve_Diff[Quantile%in%c(0.01,0.05,0.1,0.25,99.75,99.9,99.95,99.99),.(`Volume`=100*mean(Diff_naive_pc), `Periods`=100*mean(Diff_c_gt_naive_pc), `Volume:sGPD`=100*mean(Diff_pc), `Periods:sGPD`=100*mean(Diff_c_gt_s_pc)),by=Quantile] reserve_xtable[,Quantile:= paste0(Quantile,"%")] print(xtable(reserve_xtable,digits = c(0,2,1,1,1,1)),include.rownames = F) ## Pinball TS and significance test #### PB_ts <- list() for(MODEL in names(Model_list)[!1:length(Model_list) %in% grep("_qr",names(Model_list))]){ rm(NodeData_temp) load(paste0("v9_forecasts/",MODEL,".Rda")) PB_ts[[MODEL]] <- copy(NodeData_temp[["A"]][BadData==F,.(issueTime,targetTime)]) for(Ni in allNs){ print(paste("Evaluation - Model:",MODEL,"::: GSP-group:",Ni)) forecast_temp <- copy(NodeData_temp[[paste0(Ni,"_qr")]][NodeData_temp[[Ni]]$BadData==F,]) actuals_temp <- copy(NodeData_temp[[Ni]][BadData==F,node_n]) pb_temp <- copy(NodeData_temp[[Ni]][BadData==F,.(issueTime,targetTime)]) for(q in 1:9/10){ pb_temp[,(paste0("q",100 * q)):=(actuals_temp - forecast_temp[[paste0("q",100 * q)]]) * q * (actuals_temp >= forecast_temp[[paste0("q",100 * q)]]) + (actuals_temp - forecast_temp[[paste0("q",100 * q)]]) * (q - 1) * (actuals_temp < forecast_temp[[paste0("q", 100 * q)]])] } pb_temp[,(Ni):=rowMeans(pb_temp[,-c(1:2)])] PB_ts[[MODEL]] <- merge(PB_ts[[MODEL]],pb_temp[,.(issueTime,targetTime,temp_name = get(Ni))], by=c("issueTime","targetTime"),all.x=T) setnames(PB_ts[[MODEL]],"temp_name",Ni) } PB_ts[[MODEL]][,All:=rowMeans(PB_ts[[MODEL]][,-c(1:2)])] } # save(PB_ts,file=paste0("v9_forecasts/Evaluation_sig.Rda")) load("v9_forecasts/Evaluation_sig.Rda") ## Bootstrap nboot <- 250 NAMES <- names(PB_ts) REF <- "Vanilla-Point" NAMES <- NAMES[which(!NAMES %in% REF)] test_start <- "2018-01-01" GSPG <- "All" boot_ts <- copy(PB_ts[[1]][issueTime>=test_start,.(issueTime,targetTime)]) for(i in 1:length(names(PB_ts))){ boot_ts <- merge(boot_ts,PB_ts[[i]][,.(issueTime,targetTime,temp_name=get(GSPG))],by=c("issueTime","targetTime"),all.x=T) setnames(boot_ts,"temp_name",names(PB_ts)[i]) } ## Block bootstrap if >1 Block <- 48 bootdata <- data.table(Sample=1:nboot) for(i in 1:nboot){ bootdata[i,c(NAMES,REF):=as.list(colMeans(abs(boot_ts[issueTime>=test_start,][ rep(sample(1:(.N-Block+1),.N-Block+1,replace = T),each=Block)+rep(0:(Block-1),.N-Block+1), mget(c(NAMES,REF))]),na.rm = T))] } plotdata <- melt(bootdata[,100*(get(REF)-.SD)/get(REF),.SDcols=NAMES], measure.vars = 1:length(NAMES),variable.name = "Method", value.name = "Skill Score") ggplot(plotdata[`Skill Score`>=0,], aes(x=reorder(Method, -`Skill Score`), y=`Skill Score`)) + ylab("Pinball Skill Score [%]") + geom_boxplot() + theme_classic() + #ylim(c(10,15)) + ggtitle(paste("Pinball Skill Score Relative to",REF)) # ggsave("") ## DM Test require(forecast) NAMES <- names(PB_ts) test_start <- "2018-01-01" for(GSPG in c("All",allNs)){ DM <- Skill <- matrix(NA,length(NAMES),length(NAMES)) for(i in 1:length(NAMES)){ for(j in 1:length(NAMES)){ Es <- merge(PB_ts[[i]][issueTime>=test_start,.(issueTime,targetTime,e1=get(GSPG))], PB_ts[[j]][issueTime>=test_start,.(issueTime,targetTime,e2=get(GSPG))], by=c("issueTime","targetTime")) # Es <- Es[,.(e1=mean(e1),e2=mean(e2)),by="issueTime"] # acf(Es[,e1-e2],na.action = na.omit,lag.max = 52*7) # acf(Es[,e1-e2],na.action = na.omit,lag.max = 14) Skill[i,j] <- Es[!is.na(e1) & !is.na(e2),100*(1- mean(e2)/mean(e1))] if(i!=j){ DM[i,j] <- Es[!is.na(e1) & !is.na(e2),dm.test(e1 = e1, e2 = e2, h = 52,power = 1,alternative = "t")$p.value] } } } colnames(Skill) <- rownames(Skill) <- colnames(DM) <- rownames(DM) <- NAMES DM_melt <- reshape2::melt(DM,measure.vars = 1:length(NAMES), value.name = "DM test p-value") Skill_melt <- data.table(reshape2::melt(Skill,measure.vars = 1:length(NAMES), value.name = "Skill")) Skill_melt[,Skill:=format(round(Skill,1),nsmall = 1)] DM_melt <- data.table(merge(DM_melt,Skill_melt,by=c("Var1","Var2"))) # require(RColorBrewer) # ggplot(DM_melt, aes(x=Var2, y=Var1,fill=`DM test p-value`))+ # geom_tile(color = "white") + # geom_text(aes(label=Skill),size=3) + # scale_fill_gradientn(colours = rev(brewer.pal(n = 11, name = "RdYlGn")), # limits=c(0,0.1),oob=scales::squish) + # ggtitle(if(GSPG=="All"){NULL}else{ # paste0("Pinball Skill Scores [%] and Significance, GSP Group ",GSPG)}) + # theme_minimal()+ xlab("Method") + ylab("Reference") + # theme(axis.text.x = element_text(angle = 20))+ # scale_x_discrete(position = "top",limits=rev) # DM_melt[,DM:=round(`DM test p-value`,digits = 3)] DM_melt[,DMCOL:="bold"] DM_melt[`DM test p-value`>=0.05,DM:="p>0.05"] DM_melt[`DM test p-value`>=0.05,DMCOL:="plain"] DM_melt[`DM test p-value`<0.05,DM:="p<0.05"] DM_melt[`DM test p-value`<0.01,DM:="p<0.01"] DM_melt[`DM test p-value`<0.001,DM:="p<0.001"] DM_melt[`DM test p-value`<1e-6,DM:="p<1e-6"] DM_melt[Var1==Var2,DM:=""] DM_melt[Var1!=Var2,Skill:=paste0(Skill,"%")] DM_melt[Var1==Var2,Skill:=""] ggplot(DM_melt, aes(x=Var2, y=Var1))+ geom_tile(color = "grey",fill="white") + geom_text(aes(label=paste0(Skill,"\n",DM),fontface=DMCOL), size=3) + # paste("(",Counts,",",Duration,")"),hjust=-1 scale_x_discrete(position = "top",limits=rev)+ ggtitle(if(GSPG=="All"){NULL}else{ paste0("Pinball Skill Scores [%] and Significance, GSP Group ",GSPG)}) + theme_minimal()+ xlab("Method") + ylab("Reference") + theme(text=element_text(family="serif"),axis.text.x = element_text(angle = 20)) ggsave(paste0("paper_sup_plots/Sig/Skill_DM_",GSPG,".pdf")) }