VAR model_cointeg2.Rmd

---
title: "VAR model_cointeg"
author: "Hyelim Park"
#date: '2018 6 28 '
output: html_document
---

```{r, warning=FALSE,message=FALSE, echo=FALSE}
# R version 3.5.1 (2018-07-02) -- "Feather Spray"
# Copyright (C) 2018 The R Foundation for Statistical Computing
# Platform: x86_64-apple-darwin15.6.0 (64-bit)

library(MASS)
library(leaps)
library("ggplot2")
library(nlme)
library(vars)
library(aod)
library(zoo)
library(plyr)
library(forecast)
library(raster)
library(plyr)
require(rgeos)
require(RColorBrewer)
library(rgdal)
library(maptools)
library(SpatialTools)
library(geosphere)
library(plotly)



setwd("~/Dropbox/TESIS/R_code/Theme 1.8")

Capture <- read.csv(file="fishery_2002_2015.csv",sep = ",",header=TRUE)
Capture <- Capture[,c(seq(-3,-29,-2))]

for(i in c(2:ncol(Capture))) {
  Capture[,i] <- as.numeric(as.character(Capture[,i]), na.pass=T)}
rm(i)


marine_areas <- data.frame("zone" = c(18,21,27,31,34,37,41,47,48,51,57,58,61,67,71,77,81,87,88),
                           "Ocean.Area" = c("Arctic Sea","Atlantic, Northwest","Atlantic, Northeast","Atlantic, Western Central","Atlantic, Eastern Central ","Mediterranean and Black Sea","Atlantic, Southwest ","Atlantic, Southeast ","Atlantic, Antarctic","Indian Ocean, Western ","Indian Ocean, Eastern ","Indian Ocean, Antarctic and Southern","Pacific, Northwest ","Pacific, Northeast ","Pacific, Western Central ","Pacific, Eastern Central","Pacific, Southwest","Pacific, Southeast ","Pacific, Antarctic"))


marine_areas <- marine_areas[order(marine_areas$Ocean.Area),]

marine_areas[,3:17] <- Capture




CRSWG84 <- CRS(" +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
setwd("~/Dropbox/data/FAO_AREAS")
fao <- readShapeSpatial("FAO_AREAS.shp",
                        proj4string = CRSWG84)

#plot(fao.union)
fao.union <- (unionSpatialPolygons(fao, fao$F_AREA))
zone <- row.names(fao.union)
fao.union <- SpatialPolygonsDataFrame(fao.union, 
                                      data = data.frame("zone" = zone,
                                                        "area" = data.frame("area" = areaPolygon(fao.union))/10^6,
                                                        row.names=row.names(fao.union)))



FAO_cap <- merge(fao.union@data, marine_areas[,c(-2,-3)], "zone")


setwd("~/Dropbox/TESIS/R_code/Theme 1.8")

sst <- read.csv("sst4.csv", sep=" ")
chl <- read.csv("chlor_a.csv", sep=" ")
pic <- read.csv("pic.csv", sep=" ")
poc <- read.csv("poc.csv", sep=" ")


makedf <- function(fao_zone){
  #fao_zone = 18
  var<- cbind.data.frame(
      (colMeans(sst[sst$fao.zone %in% fao_zone,-1])),
      (colMeans(chl[chl$fao.zone %in% fao_zone,-1])),
      ((colMeans(pic[pic$fao.zone %in% fao_zone,-1])*100090)/(colMeans(poc[poc$fao.zone %in% fao_zone,-1]))))
    
    colnames(var) <- c("sst","chl","pic.poc")
    
    var2 <- cbind.data.frame((colSums(FAO_cap[FAO_cap$zone %in% fao_zone,-1:-2]))*19.78426/sum(FAO_cap$area[FAO_cap$zone %in% fao_zone]),
                               var[-15,] , 2002:2015)
  colnames(var2) <- c("tasa_FAO", "sst", "chl","pic.poc","year")
  
  return(var2)}


makedf_ts <- function(fao_zone){
  #fao_zone = 21
  var<- cbind.data.frame(
      (colMeans(sst[sst$fao.zone %in% fao_zone,-1])),
      (colMeans(chl[chl$fao.zone %in% fao_zone,-1])),
      ((colMeans(pic[pic$fao.zone %in% fao_zone,-1])*100090)/(colMeans(poc[poc$fao.zone %in% fao_zone,-1]))))
    
    colnames(var) <- c("sst","chl","pic.poc")
    
    var2 <- cbind.data.frame((colSums(FAO_cap[FAO_cap$zone %in% fao_zone,-1:-2]))*19.78426/sum(FAO_cap$area[FAO_cap$zone %in% fao_zone]),
                               var[-15,] , 2002:2015)
  colnames(var2) <- c("tasa_FAO", "sst", "chl","pic.poc","year")
  
  
  var2 <- ts(var2[,-5], start=c(2002), end=c(2015))
  var3 <- apply(var2[,-5], 2, function(c){
    c <- ts(c, start=c(2002), end=c(2015))
    #print(c)
    nd <-ndiffs(c)
   # cat(colnames(c), nd)
    if (nd > 0) {
      col.modi <-c(rep(NA, nd),diff(c,differences=nd))}
    else
      col.modi <- c
    return(col.modi)
    
  })
     var3 <-ts(as.data.frame(var3), start=c(2002), end=c(2015))
   var3 <- na.omit(var3)
  return(var3)}

makedf_ts_conNA <- function(fao_zone){
  #fao_zone = 27
 var<- cbind.data.frame(
      (colMeans(sst[sst$fao.zone %in% fao_zone,-1])),
      (colMeans(chl[chl$fao.zone %in% fao_zone,-1])),
      ((colMeans(pic[pic$fao.zone %in% fao_zone,-1])*100090)/(colMeans(poc[poc$fao.zone %in% fao_zone,-1]))))
    
    colnames(var) <- c("sst","chl","pic.poc")
    
    var2 <- cbind.data.frame((colSums(FAO_cap[FAO_cap$zone %in% fao_zone,-1:-2]))*19.78426/sum(FAO_cap$area[FAO_cap$zone %in% fao_zone]),
                               var[-15,] , 2002:2015)
  colnames(var2) <- c("tasa_FAO", "sst", "chl","pic.poc","year")
  var2 <- ts(var2[,-5], start=c(2002), end=c(2015))
  var3 <- apply(var2[,-5], 2, function(c){
    c <- ts(c, start=c(2002), end=c(2015))
    #print(c)
    nd <- ndiffs(c)
    #print(nd)
    if (nd > 0) {
      col.modi <-c(rep(NA, nd),diff(c,differences=nd))}
    else
      col.modi <- c
    return(col.modi)
    
  })
     var3 <-ts(as.data.frame(var3), start=c(2002), end=c(2015))
  # var3 <- na.omit(var3)
  return(var3)}




```




```{r,echo=FALSE, warning=FALSE, message=FALSE}


cointeg <- function(fao_zone){
 # cat("fao area is ",fao_zone, "\n")
  
  # fao_zone <-(81)
  dff <- makedf(fao_zone)
 # dff_ts <- makedf_ts(fao_zone)
  var_combi <- data.frame("id"=1:7)
   if(length(fao_zone) == 1) {
  var_combi$zone <- rep(fao_zone, nrow(var_combi))} else {
     var_combi$zone <- rep("Global", nrow(var_combi))
  }
  
  #var_combi <- var_combi[,-1]
  vv<- rbind(ldply(list(2,3,4,c(2,3),c(2,4),c(3,4),c(2,3,4)),function(l) {
  #cat(l, length(l),"\n")
    #l=2
  vv <- data.frame(matrix(c(colnames(dff)[c(l)],rep(NA,3-length(l))),ncol = 3),stringsAsFactors=FALSE)
  #print(vv)
  
   cajo <- ca.jo(dff[,c(1,l)], type="eigen",spec="transitory",ecdet="none",K=2)
    S <- summary(cajo)
    valo <- as.data.frame(cbind(S@teststat, S@cval))
  
  if(length(which(valo[,1]>valo[,4])) >0 ){
      r <- rownames(valo[which(valo[,1]>valo[,3]),])
      vv$rank <-max(as.numeric(substr(r, start=5, stop=6)))+1
     
     } else{
       vv$rank <- 0
    }
    return(vv)}))
  

  colnames(vv)[1:3]<-c("variables_1","variables_2","variables_3")  
  
  var_combi <- cbind(var_combi,vv)
    
  var_level <- mdply(var_combi[,c(-1,-6)], function(zone, variables_1, variables_2, variables_3){
  i_d <-as.data.frame(t( apply(dff[,-5],2,  function(c){
    c <- ts(c, start=c(2002), end=c(2015))
    nd <-ndiffs(c)
  })))
  max_id <- max(i_d)
  vari <- c(na.omit(c(variables_1, variables_2, variables_3)))
  #cat(zone, vari, "\n")
  
  vv <-  VAR(dff[,c("tasa_FAO",vari)], lag.max = 2, ic = c("AIC"))
  ff <- tryCatch(summary(vv)$varresult$tasa_FAO$fstatistic,
                 error = function(e) as.numeric(NA))
  #cat("ok...\n")
  pvalue <-tryCatch(pf(ff[1],ff[2],ff[3],lower.tail=F),
                    error = function(e) as.numeric(NA))
  rmse_value <-  tryCatch(accuracy(vv$varresult[[1]])[,2],
                          error = function(e) as.numeric(NA))
  normality <-  tryCatch(serial.test(vv, lags.bg=2, type="ES")$serial$p.value,
                         error = function(e) as.numeric(NA))
  jbtest <-  tryCatch(normality.test(vv)$jb.mul$JB$p.value,
                         error = function(e) as.numeric(NA))

  bic_value <- CV(vv$varresult$tasa_FAO)[[c(4)]]
  r2 <- CV(vv$varresult$tasa_FAO)[[c(5)]]
  causal <- tryCatch(causality(vv, cause = colnames(vv$y)[-1])$Granger$p.value,
                     error = function(e) NA) 
  
  lt <- data.frame("max_I(d)"=max_id,"lag"=vv$p,"p_value"=pvalue, "R2"=r2, "BIC"=bic_value,
                   "RMSE"=rmse_value, "normality"=normality, "JB-Test"=jbtest,
                   "causal" = causal)
 # cat("MIAU\n")
  
  lt <- cbind(i_d,lt)
  return(lt)
})  
    
   
  
 var_level <- cbind(var_level,var_combi[6])
  
  
  
  return(var_level)}

```



# FAO area scale
```{r,echo=FALSE, warning=FALSE, message=FALSE}

resume <- ldply(levels(FAO_cap$zone), function(z){
  df <- makedf(z)[,-5]
  dd <- data.frame("zone" = c(z,z),"test"=c("mean","sd"))
  dd[,2:5+1] <- rbind(apply(df, 2, mean),
  apply(df, 2, sd))
  colnames(dd)[2:5+1] <- c("tasa_FAO", "sst","chl","pic.poc")
  return(dd)})





lm_resume <- ldply(levels(FAO_cap$zone), function(z){
  #z=18
  df <- makedf(z)
  dd <- data.frame("zone" = z)
  dd$beta_tasaFAO<- c(summary(lm(tasa_FAO~year, df))$coeff[2,1])
  dd$p_tasaFAO <- summary(lm(tasa_FAO~year, df))$coeff[2,4]
  dd$beta_sst<- c(summary(lm(sst~year, df))$coeff[2,1])
  dd$p_sst <- summary(lm(sst~year, df))$coeff[2,4]
  dd$beta_chl<- c(summary(lm(chl~year, df))$coeff[2,1])
  dd$p_chl <- summary(lm(chl~year, df))$coeff[2,4]
  dd$beta_picpoc<- c(summary(lm(pic.poc~year, df))$coeff[2,1])
  dd$p_picpoc <- summary(lm(pic.poc~year, df))$coeff[2,4]
  return(dd)
})

lm_resume[lm_resume$p_tasaFAO < 0.05,1:3]
lm_resume[lm_resume$p_sst < 0.05,c(1,4:5)]
lm_resume[lm_resume$p_chl < 0.05,c(1,6:7)]
lm_resume[lm_resume$p_picpoc < 0.05,c(1,8:9)]




test <- ldply(levels(FAO_cap$zone), function(z){
  #cat("FAO area:", z, "\n");
  selec <- cointeg(z);
  #  print(selec)
})

#test


test_selec <- test[test$normality > 0.05 &  test$causal < 0.05,]

test <- test_selec[!(is.na(test_selec$zone)),]

test[order(test$BIC),]




```


## Area 21 (VAR): both I(1), no integ. ->VAR in dff : no sig.
```{r,echo=FALSE, warning=FALSE, message=FALSE}
cat("CHANGED!!!!\n")
level <- makedf(21)
dff <- makedf_ts(21)
cajo<- ca.jo(level[,c("tasa_FAO","chl")], type="eigen",spec="transitory",ecdet="none",K=2)
summary(cajo)

model_21 <-  VAR(dff[,c("tasa_FAO","chl")], p=1)
summary(model_21)

CV(model_21$varresult$tasa_FAO)
CV(model_21$varresult$chl)

accuracy(model_21$varresult$tasa_FAO)
accuracy(model_21$varresult$chl)


causality(model_21, cause = "chl")
causality(model_21, cause = "tasa_FAO")


serial.test(model_21, lags.bg=2, type="ES")
normality.test(model_21)


```

## Area 31 (VAR): different I(d), no cointeg.
```{r,echo=FALSE, warning=FALSE, message=FALSE}
level <- makedf(31)
model_31 <-  VAR(level[,c("tasa_FAO","sst")], lag.max = 2, ic = c("AIC"))
summary(model_31)

CV(model_31$varresult$tasa_FAO)
CV(model_31$varresult$sst)

accuracy(model_31$varresult$tasa_FAO)
accuracy(model_31$varresult$sst)


causality(model_31, cause = "sst")
causality(model_31, cause = "tasa_FAO")


serial.test(model_31, lags.bg=2, type="ES")
normality.test(model_31)

```


## Area 41 (VAR): different I(1), no cointeg. but no one-way causality -> VAR sst
```{r,echo=FALSE, warning=FALSE, message=FALSE}

level <- makedf(41)
dff <- makedf_ts(41)

cajo<- ca.jo(level[,c("tasa_FAO","sst")], type="eigen",spec="transitory",ecdet="none",K=2)
summary(cajo)

model_41 <-  VAR(level[,c("tasa_FAO","sst")], lag.max = 2, ic = c("AIC"))

summary(model_41)


CV(model_41$varresult$tasa_FAO)
CV(model_41$varresult$sst)

accuracy(model_41$varresult$tasa_FAO)
accuracy(model_41$varresult$sst)


causality(model_41, cause = "sst")
causality(model_41, cause = "tasa_FAO")


serial.test(model_41, lags.bg=2, type="ES")
normality.test(model_41)

```



## Area 57 (VAR): different order, no cointeg. 
```{r,echo=FALSE, warning=FALSE, message=FALSE}
level <- makedf(57)
dff <- makedf_ts(57)
model_57 <-  VAR(level[,c("tasa_FAO","chl")], lag.max = 2, ic = c("AIC"))

summary(model_57)

CV(model_57$varresult$tasa_FAO)
CV(model_57$varresult$chl)

accuracy(model_57$varresult$tasa_FAO)
accuracy(model_57$varresult$chl)


causality(model_57, cause = "chl")
causality(model_57, cause = "tasa_FAO")


serial.test(model_57, lags.bg=2, type="ES")
normality.test(model_57)

```

## Area 58 (VAR): all I(0), no cointeg. 
```{r,echo=FALSE, warning=FALSE, message=FALSE}
level <- makedf(58)
dff <- makedf_ts(58)
model_58 <-  VAR(level[,c("tasa_FAO","chl")], lag.max = 2, ic = c("AIC"))

summary(model_58)

CV(model_58$varresult$tasa_FAO)
CV(model_58$varresult$chl)

accuracy(model_58$varresult$tasa_FAO)
accuracy(model_58$varresult$chl)


causality(model_58, cause = "chl")
causality(model_58, cause = "tasa_FAO")


serial.test(model_58, lags.bg=2, type="ES")
normality.test(model_58)

```


## Area 77 (VAR): same I(d) but d=0, 1 cointeg.
```{r,echo=FALSE, warning=FALSE, message=FALSE}
cat("CHANGED!!!!!\n")

level <- makedf(77)
model_77 <-  VAR(level[,c("tasa_FAO","sst")], lag.max = 2, ic = c("AIC"))
summary(model_77)

causality(model_77, cause = "sst")
causality(model_77, cause = "tasa_FAO")


CV(model_77$varresult$tasa_FAO)
CV(model_77$varresult$sst)

accuracy(model_77$varresult$tasa_FAO)
accuracy(model_77$varresult$sst)


causality(model_77, cause = "sst")
causality(model_77, cause = "tasa_FAO")


serial.test(model_77, lags.bg=2, type="ES")
normality.test(model_77)

```
## Area 88 (VAR): diff I(d), no cointeg. 
```{r,echo=FALSE, warning=FALSE, message=FALSE}
level <- makedf(88)
dff <- makedf_ts(88)
model_88 <-  VAR(level[,c("tasa_FAO","pic.poc")], lag.max = 2, ic = c("AIC"))

summary(model_88)

CV(model_88$varresult$tasa_FAO)
CV(model_88$varresult$pic.poc)

accuracy(model_88$varresult$tasa_FAO)
accuracy(model_88$varresult$pic.poc)


causality(model_88, cause = "pic.poc")
causality(model_88, cause = "tasa_FAO")


serial.test(model_88, lags.bg=2, type="ES")
normality.test(model_88)

```
## Forcasting tasa_FAO
```{r,echo=FALSE, warning=FALSE, message=FALSE}

single_model <- mdply(cbind(list(model_31, model_41,model_57,model_58, model_77,model_88),
                            c("31","41","57","58","77","88")),
                      function(m,c){
                        endo_v <- makedf(c)$tasa_FAO
                        endo <- cbind.data.frame(endo_v,endo_v,
                                 endo_v,endo_v)
                        
                        pred <- predict(m,n.ahead=5, ci=0.95)
                       
                        colnames(endo) <- c("fcst","lower","upper","CI")
                        df <- rbind(endo,as.data.frame(pred$fcst$tasa_FAO))

                        df$year <- 2002:2020
                        df$zone <- paste("FAO ",rep(c,19),sep="")
                        
                        return(df)
                      })

single_model <- single_model[,-1]
#colnames(single_model)

lm_test <- lapply(c("FAO 31","FAO 41","FAO 57","FAO 58","FAO 77","FAO 88"),function(z) { 
  dd <- single_model[single_model$zone %in% z,]
  ss <- summary(lm(fcst ~ year, dd))
 
  ff <- ss$fstatistic
  print(z)
  print(ss$coefficients)
 # print(pf(ff[1],ff[2],ff[3],lower.tail=F))
  })

single_model$lm_line <- unlist(lapply(c("FAO 31","FAO 41","FAO 57","FAO 58","FAO 77","FAO 88"),function(z) { 
 
  dd <- single_model[single_model$zone %in% z,]
  ss <- lm(fcst ~ year, dd)
 
  return(ss$fitted.values)
  }))


slope <- unlist(lapply(lm_test, function(f){
  if(f[2,4]>0.01&f[2,4]<0.05){p=paste(signif(f[2,1], digits =3),"*",sep = "")}
  if(f[2,4]>0.001&f[2,4]<0.01){p=paste(signif(f[2,1], digits =3),"**",sep = "")}
  if(f[2,4]<0.001){p=paste(signif(f[2,1], digits =3),"***",sep = "")}
  if(0.05<f[2,4]){p=paste(signif(f[2,1], digits =3),"",sep = "")}
  
  return(p)}))

annotate_slope <- list(
  x = c(2018,2016,2017,2015,2019:2020),
  y = single_model$lm_line[c(17,34,54,71,94,114)],
  text = slope,
  xref = "year",
  yref = "fcst",
  showarrow = TRUE,
  arrowhead = 7,
  ax = 20,
  ay = -40
)





gvar <- plot_ly(single_model,x = ~year, y = ~fcst, color = ~zone, mode = 'markers', type = 'scatter',showlegend=F)%>%
  layout(yaxis = list (title = "Capture rate (tonnes/km^2)"),
         xaxis = list (title = "Year"))%>%
  
  add_ribbons(showlegend = F,
              ymin = ~lower,
              ymax = ~upper,
              line = list(color = 'rgba(7, 164, 181, 0.05)'),
              fillcolor = 'rgba(7, 164, 181, 0.2)') %>%
  layout(annotations = annotate_slope)%>%
  add_trace(y = ~lm_line,mode = 'lines',showlegend=T)

gvar

#Sys.setenv("plotly_username"="pinguino.lim")
#Sys.setenv("plotly_api_key"="wMn1rceswv2Jabna8G3K")
#api_create(gvar, filename = "forecasting_FAO")

#fao.union@data[order(fao.union$X2014, decreasing = T),]


```
## Forcasting other variable
```{echo=FALSE, warning=FALSE, message=FALSE}
single_model <- mdply(cbind(list(model_21, model_31, model_41,model_57,model_58,model_67,model_67,model_77,model_88),
                            c("21","31","41","57","58","67_chl","67_pp","77","88"),
                            c("chl","sst","sst",
                      "chl","chl","chl","pic.poc","sst","pic.poc")),
                      function(m,c,p){
                        pred <- predict(m,n.ahead=5, ci=0.95)
                        endo <- cbind.data.frame(as.data.frame(pred$endog)[p],
                                               as.data.frame(pred$endog)[p],
                                                as.data.frame(pred$endog)[p],
                                               as.data.frame(pred$endog)[p])
                        colnames(endo) <- c("fcst","lower","upper","CI")
                        df <- as.data.frame(rbind(endo,eval(parse(text=sprintf("pred$fcst$%s", p)))))
                        df$year <- 2002:2020
                        df$zone <- rep(c,19)
                        df$variable <- p
                        
                        return(df)
                      })

single_model <- single_model[,-1]
#colnames(single_model)

lm_test <- lapply(c("21","31","41","57","58","67_chl","67_pp","77","88"),function(z) { 
  dd <- single_model[single_model$zone %in% z,]
  ss <- summary(lm(fcst ~ year, dd))
 
  ff <- ss$fstatistic
  print(z)
  print(ss$coefficients)
 # print(pf(ff[1],ff[2],ff[3],lower.tail=F))
  })




library(plotly)
#packageVersion('plotly')



gvar <- plot_ly(single_model,x = ~year, y = ~fcst, color = ~zone)%>%
  layout(yaxis = list (title = "Capture rate (tonnes/km^2)"),
         xaxis = list (title = "Year"))%>%
  add_lines()%>%
  add_ribbons(showlegend = F,
              ymin = ~lower,
              ymax = ~upper,
              line = list(color = 'rgba(7, 164, 181, 0.05)'),
              fillcolor = 'rgba(7, 164, 181, 0.2)')

gvar

#Sys.setenv("plotly_username"="pinguino.lim")
#Sys.setenv("plotly_api_key"="wMn1rceswv2Jabna8G3K")

#fao.union@data[order(fao.union$X2014, decreasing = T),]


```


#Global scale

```{r,echo=FALSE, warning=FALSE, message=FALSE}
global_test <- cointeg(levels(FAO_cap$zone))

global_test[ global_test$p_value < 0.05 ,]

global_test_vecm <- global_test[global_test$rank != 0,]


```
##VECM_Global no exist -> VAR
```{r,echo=FALSE, warning=FALSE, message=FALSE}

global_level <- makedf(levels(FAO_cap$zone))
#global_diff <- makedf_ts(levels(FAO_cap$zone))

colMeans(global_level)
apply(global_level,2,FUN=sd)

summary(lm(tasa_FAO ~ year, as.data.frame(global_level)))
summary(lm(sst ~ year, as.data.frame(global_level)))

summary(lm(chl ~ year, as.data.frame(global_level)))
summary(lm(pic.poc ~ year, as.data.frame(global_level)))

cajo<- ca.jo(global_level[,c("tasa_FAO","sst")], type="eigen",spec="transitory",ecdet="none",K=2)
summary(cajo)



print("global_level")
global_model <-  VAR(global_level[,c("tasa_FAO","sst")], lag.max = 2, ic = c("AIC"))
summary(global_model)
CV(global_model$varresult$tasa_FAO)
CV(global_model$varresult$sst)

accuracy(global_model$varresult$tasa_FAO)
accuracy(global_model$varresult$sst)


causality(global_model, cause = "sst")
causality(global_model, cause = "tasa_FAO")


serial.test(global_model, lags.bg=2, type="ES")
normality.test(global_model)









#plot(predict(global_model,n.ahead=5, ci=0.95))

single_g_model <- mdply(cbind(list(global_model),
                            c("Global")),
                      function(m,c){
                        pred <- predict(m,n.ahead=5, ci=0.95)
                        endo <- cbind.data.frame(pred$endog[,1],pred$endog[,1],
                                                 pred$endog[,1],pred$endog[,1])
                        colnames(endo) <- c("fcst","lower","upper","CI")
                        df <- as.data.frame(rbind(endo,pred$fcst$tasa_FAO))
                        df$year <- 2002:2020
                        df$zone <- rep(c,19)
                        
                        return(df)
                      })

single_g_model <- single_g_model[,-1]
#colnames(single_model)

lm_test <- lapply(c("Global"),function(z) { 
  dd <- single_g_model[single_g_model$zone %in% z,]
  ss <- summary(lm(fcst ~ year, dd))
 
  ff <- ss$fstatistic
  print(z)
  print(ss$coefficients)
 # print(pf(ff[1],ff[2],ff[3],lower.tail=F))
  })

single_g_model$lm_line <- unlist(lapply(c("Global"),function(z) { 
  dd <- single_g_model[single_g_model$zone %in% z,]
  ss <- lm(fcst ~ year, dd)
 
  return(ss$fitted.values)
  }))

slope <- unlist(lapply(lm_test, function(f){
  if(f[2,4]>0.01&f[2,4]<0.05){p=paste(signif(f[2,1], digits =3),"*",sep = "")}
  if(f[2,4]>0.001&f[2,4]<0.01){p=paste(signif(f[2,1], digits =3),"**",sep = "")}
  if(f[2,4]<0.001){p=paste(signif(f[2,1], digits =3),"***",sep = "")}
  if(0.05<f[2,4]){p=paste(signif(f[2,1], digits =3),"",sep = "")}
  
  return(p)}))

annotate_slope <- list(
  x = 2017,
  y = single_g_model$lm_line[single_g_model$year %in% 2017],
  text = slope,
  xref = "year",
  yref = "fcst",
  showarrow = TRUE,
  arrowhead = 7,
  ax = 20,
  ay = -40
)

library(plotly)
#packageVersion('plotly')



gvar <- plot_ly(single_g_model,x = ~year, y = ~fcst, color = ~zone, 
                mode = 'markers', type = 'scatter',showlegend=F)%>%
  layout(yaxis = list (title = "Capture rate (tonnes/km^2)"),
         xaxis = list (title = "Year"))%>%
  add_ribbons(showlegend = F,
              ymin = ~lower,
              ymax = ~upper,
              line = list(color = 'rgba(7, 164, 181, 0.05)'),
              fillcolor = 'rgba(7, 164, 181, 0.2)') %>%
  layout(annotations = annotate_slope)%>%
  add_trace(y = ~lm_line,mode = 'lines',showlegend=T)

gvar

#api_create(gvar, filename = "forecasting_Global")

```