library(ncdf4)
library(lattice)
library(Metrics)
library(multDM)

mydatar.nc <- nc_open("T2_daily_1982_2years")
mydatar <- ncvar_get(mydatar.nc,"T2")
lon <- mydatar.nc$dim$x$vals
lat <- mydatar.nc$dim$y$vals
dim(mydatar)

mydata1.nc <- nc_open("T2_daily_1982_0hour")
mydata1 <- ncvar_get(mydata1.nc,"T2")

mydata2.nc <- nc_open("T2_daily_1982_6months")
mydata2 <- ncvar_get(mydata2.nc,"T2")

ncol <- length(lat)
nrow <- length(lon)

nreplications=1000
pval=0.05

rmse1 <- matrix(data=NA,nrow=nrow,ncol=ncol)
rmse2 <- matrix(data=NA,nrow=nrow,ncol=ncol)
rmsed <- matrix(data=NA,nrow=nrow,ncol=ncol)
rmsedtpc <- matrix(data=NA,nrow=nrow,ncol=ncol)
dm <- matrix(data=NA,nrow=nrow,ncol=ncol)
bootstrap <- matrix(data=NA,nrow=nrow,ncol=ncol)

q1=pval/2
q2=1-q1

for (i in 1:nrow){
for (j in 1:ncol){

  rmse1[i,j] = rmse(mydatar[i,j,],mydata1[i,j,])
  rmse2[i,j] = rmse(mydatar[i,j,],mydata2[i,j,])
  rmsed[i,j] = rmse2[i,j] - rmse1[i,j]  ### absoluite diff
  rmsedtpc[i,j] = ((rmse2[i,j] - rmse1[i,j]) / rmse1[i,j])*100  ### % diff
  
  aux <- DM.test(mydata1[i,j,], mydata2[i,j,], mydatar[i,j,], loss.type="SE",c=FALSE,H1="same")
  dm[i,j] = aux$p.value
  
  rmsedrandom=replicate(nreplications, NA)
  xy=c(mydata1[i,j,],mydata2[i,j,])
  for (k in 1:nreplications) {
    idx=sample(1:length(xy), replace=F)
    x_b=xy[idx[1:(length(xy)/2)]]
    y_b=xy[idx[(length(xy)/2+1):length(xy)]]
    rmsedrandom[k] = rmse(x_b,mydatar[i,j,])-rmse(y_b,mydatar[i,j,])
  }
  th=quantile(rmsedrandom,c(q1,q2))
  if (rmsed[i,j]<th[1] || rmsed[i,j]>th[2]) {
    bootstrap[i,j]=1
  } else {
    bootstrap[i,j]=0
  }  
  
}
}