#BIG TROUT LAKE
dat1=read.csv("mm6010735.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat5=cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(53.83,length(dat1$Year)),rep(-89.87,length(dat1$Year)),rep(224,length(dat1$Year)))

dat5=dat5[dat5[,1]>=1999 & dat5[,1]<=2008,]
dat5=dat5[-c(8,10),]
#mean prediction
newdata=data.frame(x1=c(1999,2000,2001,2002,2003,2004,2005,2007),x2=dat5[,14],x3=dat5[,15],x4=dat5[,16])

newdata1=data.frame(latit=rep(dat5[,14],12),longi=rep(dat5[,15],12),eleva=rep(dat5[,16],12),years=rep(newdata$x1,12),mon=c(rep(1,length(dat5[,1])),rep(2,length(dat5[,1])),rep(3,length(dat5[,1])),rep(4,length(dat5[,1])),rep(5,length(dat5[,1])),rep(6,length(dat5[,1])),rep(7,length(dat5[,1])),rep(8,length(dat5[,1])),rep(9,length(dat5[,1])),rep(10,length(dat5[,1])),rep(11,length(dat5[,1])),rep(12,length(dat5[,1]))),tindex=c((newdata$x1-1978)*12+1,(newdata$x1-1978)*12+2,(newdata$x1-1978)*12+3,(newdata$x1-1978)*12+4,(newdata$x1-1978)*12+5,(newdata$x1-1978)*12+6,(newdata$x1-1978)*12+7,(newdata$x1-1978)*12+8,(newdata$x1-1978)*12+9,(newdata$x1-1978)*12+10,(newdata$x1-1978)*12+11,(newdata$x1-1978)*12+12))


mod3prediction=predict(fit1,newdata1)
mod4prediction=predict(fit2,newdata1)

mod3prediction=matrix(mod3prediction,length(dat5[,1]),length(mod3prediction)/length(dat5[,1]))
mod4prediction=matrix(mod4prediction,length(dat5[,1]),length(mod4prediction)/length(dat5[,1]))

mean1=predict(model1,newdata)
mean2=predict(model2,newdata)
mean3=predict(model3,newdata)
mean4=predict(model4,newdata)
mean5=predict(model5,newdata)
mean6=predict(model6,newdata)
mean7=predict(model7,newdata)
mean8=predict(model8,newdata)
mean9=predict(model9,newdata)
mean10=predict(model10,newdata)
mean11=predict(model11,newdata)
mean12=predict(model12,newdata)

meanmat=cbind(mean1,mean2,mean3,mean4,mean5,mean6,mean7,mean8,mean9,mean10,mean11,mean12)
marginalvar=c(op1$par[6],op2$par[6],op3$par[7],op4$par[7],op5$par[6],op6$par[7],op7$par[7],op8$par[7],op9$par[7],op10$par[6],op11$par[6],op12$par[6])

resipre1<-dat5[,2]-mean1
resipre2<-dat5[,3]-mean2
resipre3<-dat5[,4]-mean3
resipre4<-dat5[,5]-mean4
resipre5<-dat5[,6]-mean5
resipre6<-dat5[,7]-mean6
resipre7<-dat5[,8]-mean7
resipre8<-dat5[,9]-mean8
resipre9<-dat5[,10]-mean9
resipre10<-dat5[,11]-mean10
resipre11<-dat5[,12]-mean11
resipre12<-dat5[,13]-mean12

resipreM=cbind(resipre1,resipre2,resipre3,resipre4,resipre5,resipre6,resipre7,resipre8,resipre9,resipre10,resipre11,resipre12)
month1<-pnorm(resipreM[,1],mean=0,sd =op1$par[6])
month2<-pnorm(resipreM[,2],mean=0,sd =op2$par[6])
month3<-pnorm(resipreM[,3],mean=0,sd =op3$par[7])
month4<-pnorm(resipreM[,4],mean=0,sd =op4$par[7])
month5<-pnorm(resipreM[,5],mean=0,sd =op5$par[6])
month6<-pnorm(resipreM[,6],mean=0,sd =op6$par[7])
month7<-pnorm(resipreM[,7],mean=0,sd =op7$par[7])
month8<-pnorm(resipreM[,8],mean=0,sd =op8$par[7])
month9<-pnorm(resipreM[,9],mean=0,sd =op9$par[7])
month10<-pnorm(resipreM[,10],mean=0,sd =op10$par[6])
month11<-pnorm(resipreM[,11],mean=0,sd =op11$par[6])
month12<-pnorm(resipreM[,12],mean=0,sd =op12$par[6])

pri<-cbind(month1,month2,month3,month4,month5,month6,month7,month8,month9,month10,month11,month12)
priT<-t(pri)
priv<-as.vector(priT)

mod1=matrix(0,length(dat5[,1]),12)
mod2=matrix(0,length(dat5[,1]),12)
mod3=matrix(0,length(dat5[,1]),12)
mod4=matrix(0,length(dat5[,1]),12)

percentage1=rep(0,12)
percentage2=rep(0,12)
percentage3=rep(0,12)

for(g in 1:12){
  if(g==1){
    
    
    mod1[,g]<-abs(mean1-dat5[,g+1])
    mod2[,g]<-abs(mean1-dat5[,g+1])
    mod3[,g]<-abs(mod3prediction[,g]-dat5[,g+1])
    mod4[,g]<-abs(mod4prediction[,g]-dat5[,g+1])
    
    
    count=0
    for(l in 1:length(dat5[,g+1])){
      if(abs(mean1[l]-dat5[l,g+1])<abs(mean1[l]-dat5[l,g+1])){
        count=count+1
      }else{
        count=count
      }
    }
    percentage1[g]<-count/length(dat5[,g+1])
    
    count1=0
    for(i in 1:length(mean1)){
      if(abs(mean1[l]-dat5[l,g+1])<abs(mod3prediction[l,1]-dat5[l,g+1])){
        count1=count1+1
      }else{
        count1=count1
      }
    }
    percentage2[g]<-count1/length(fitted)
    
    count2=0
    for(i in 1:length(mean1)){
      if(abs(mean1[l]-dat5[l,g+1])<abs(mod4prediction[l,1]-dat5[l,g+1])){
        count2=count2+1
      }else{
        count2=count2
      }
    }
    percentage3[g]<-count2/length(fitted)
    
    
  }else if(g==2){
    RVM5<-RVineMatrix(RVM3$Matrix[(13-g):12,(13-g):12], family = RVM3$family[(13-g):12,(13-g):12], par = RVM3$par[(13-g):12,(13-g):12], par2 = RVM3$par2[(13-g):12,(13-g):12],check.pars = TRUE)
    dat<-c(rep(0,g))
    
    fitted<-rep(0,length(dat5[,g+1]))
    fitted1<-rep(0,length(dat5[,g+1]))
    fitted2<-rep(0,length(dat5[,g+1]))
    for(j in 1:length(dat5[,g+1])){
      a<-seq(-30,30,0.01)
      b<-rep(0,6001)
      c<-rep(0,6001)
      
      
      
      c2<-function(k){
        cop<-pnorm(k,mean=0,sd =op2$par[6])
        dat<-c(pri[j,1:(g-1)],cop)
        RVinePDF(dat,RVM5,verbose=TRUE)
      }
      
      f1<-function(k){
        dnorm(k,mean=0,sd =op2$par[6])
      }
      
      f<-function(k){
        k*f1(k)
      }
      
      
      
      
      #fitted
      for(i in 1:6001){
        b[i]<-c2(a[i])
        c[i]<-f1(a[i])
        d<-a*b*c*0.01
      }
      
      fitted[j]<-sum(d)+mean2[j]
      
      #fitted 1
      
      fitted1[j]<-integrate(f,-Inf,Inf)$value+mean2[j]
      
      
    }
    
    
    mod1[,g]<-abs(fitted-dat5[,g+1])
    mod2[,g]<-abs(fitted1-dat5[,g+1])
    mod3[,g]<-abs(mod3prediction[,g]-dat5[,g+1])
    mod4[,g]<-abs(mod4prediction[,g]-dat5[,g+1])
    
    
    count=0
    for(l in 1:length(fitted)){
      if(abs(fitted[l]-dat5[l,g+1])<abs(fitted1[l]-dat5[l,g+1])){
        count=count+1
      }else{
        count=count
      }
    }
    percentage1[g]<-count/length(fitted)
    
    count1=0
    for(i in 1:length(fitted)){
      if(abs(fitted[i]-dat5[l,g+1])<abs(mod3prediction[l,g]-dat5[l,g+1])){
        count1=count1+1
      }else{
        count1=count1
      }
    }
    percentage2[g]<-count1/length(fitted)
    
    count2=0
    for(i in 1:length(fitted)){
      if(abs(fitted[i]-dat5[l,g+1])<abs(mod4prediction[l,g]-dat5[l,g+1])){
        count2=count2+1
      }else{
        count2=count2
      }
    }
    percentage3[g]<-count2/length(fitted)
    
    
  }else{
    RVM5<-RVineMatrix(RVM3$Matrix[(13-g):12,(13-g):12], family = RVM3$family[(13-g):12,(13-g):12], par = RVM3$par[(13-g):12,(13-g):12], par2 = RVM3$par2[(13-g):12,(13-g):12],check.pars = TRUE)
    
    
    RVM4<-RVineMatrix(RVM3$Matrix[(14-g):12,(14-g):12], family = RVM3$family[(14-g):12,(14-g):12], par = RVM3$par[(14-g):12,(14-g):12], par2 = RVM3$par2[(14-g):12,(14-g):12],check.pars = TRUE)
    
    
    dat<-c(rep(0,g))
    
    fitted<-rep(0,length(dat5[,g+1]))
    fitted1<-rep(0,length(dat5[,g+1]))
    fitted2<-rep(0,length(dat5[,g+1]))
    for(j in 1:length(dat5[,g+1])){
      a<-seq(-30,30,0.01)
      b<-rep(0,6001)
      c<-rep(0,6001)
      
      c1<-rep(0,length(dat5[,g+1]))
      c1[j]<-RVinePDF(pri[j,1:(g-1)], RVM4, verbose = TRUE)
      
      
      c2<-function(k){
        cop<-pnorm(k,mean=0,sd =marginalvar[g])
        dat<-c(pri[j,1:(g-1)],cop)
        RVinePDF(dat,RVM5,verbose=TRUE)
      }
      
      f1<-function(k){
        dnorm(k,mean=0,sd =marginalvar[g])
      }
      
      f<-function(k){
        k*f1(k)
      }
      
      
      
      
      #fitted
      for(i in 1:6001){
        b[i]<-c2(a[i])
        c[i]<-f1(a[i])
        d<-a*b*c*0.01/c1[j]
      }
      
      fitted[j]<-sum(d)+meanmat[j,g]
      
      #fitted 1
      
      fitted1[j]<-integrate(f,-Inf,Inf)$value+meanmat[j,g]
      
      
      
      
      
      
      
      
      
      
      
    }
    
    mod1[,g]<-abs(fitted-dat5[,g+1])
    mod2[,g]<-abs(fitted1-dat5[,g+1])
    mod3[,g]<-abs(mod3prediction[,g]-dat5[,g+1])
    mod4[,g]<-abs(mod4prediction[,g]-dat5[,g+1])
    
    
    count=0
    for(l in 1:length(fitted)){
      if(abs(fitted[l]-dat5[l,g+1])<abs(fitted1[l]-dat5[l,g+1])){
        count=count+1
      }else{
        count=count
      }
    }
    percentage1[g]<-count/length(fitted)
    
    count1=0
    for(i in 1:length(fitted)){
      if(abs(fitted[i]-dat5[i,g+1])<abs(mod3prediction[i,g]-dat5[i,g+1])){
        count1=count1+1
      }else{
        count1=count1
      }
    }
    percentage2[g]<-count1/length(fitted)
    
    count2=0
    for(i in 1:length(fitted)){
      if(abs(fitted[i]-dat5[i,g+1])<abs(mod4prediction[i,g]-dat5[i,g+1])){
        count2=count2+1
      }else{
        count2=count2
      }
    }
    percentage3[g]<-count2/length(fitted)
  }
  
  
  
}

mod1sta=rep(0,length(dat5[,1]))
mod2sta=rep(0,length(dat5[,1]))
mod3sta=rep(0,length(dat5[,1]))
mod4sta=rep(0,length(dat5[,1]))

percentage1=rep(0,length(dat5[,1]))
percentage2=rep(0,length(dat5[,1]))
percentage3=rep(0,length(dat5[,1]))
for(i in 1:length(dat5[,1])){
  mod1sta[i]=mean(mod1[i,])
  mod2sta[i]=mean(mod2[i,])
  mod3sta[i]=mean(mod3[i,])
  mod4sta[i]=mean(mod4[i,])
  
  count=0
  count1=0
  count2=0
  for(j in 1:12){
    if(mod1[i,j]<=mod2[i,j]){
      count=count+1
    }
    if(mod1[i,j]<=mod3[i,j]){
      count1=count1+1
    }
    if(mod1[i,j]<=mod4[i,j]){
      count2=count2+1
    }
  }
  percentage1[i]=count/12
  percentage2[i]=count1/12
  percentage3[i]=count2/12
}


