# https://github.com/pcm-dpc/COVID-19/tree/master/dati-regioni
# https://medium.com/@tomaspueyo/coronavirus-act-today-or-people-will-die-f4d3d9cd99ca


# install.packages("TSdist")

 library(forecast)
 library(waveslim)
 library(xts)
 library("meboot")
 library(RTransferEntropy)
 library("meboot")
library("TSdist")

	setwd("G:\\Lavori\\covid\\")

   y 	   		<-  read.table( "COVID_13_3_2020A.txt", sep="\t", header = T) 
	head(y,1) 

	names(y)
 
	xx 		<-  y[, -c( 1,3)]
	head(xx,1) 
      names(xx)
      CasiVsTamp <-  xx[,10]/xx[,11]  
	DecVsCasi  <-  xx[,9 ]/xx[,10]  * 100

	x <- cbind(xx, CasiVsTamp, DecVsCasi   )
head(x, 1)
attach(x)

 		dates 	<- seq(as.Date("2020-02-24"), length=18, by="day")

nobs <-  length( dates )

		newdata 	<- x[ which(codice_regione=="1"),]
 		PIEMONTE 	<- xts( newdata , order.by=dates )


		as.vector( PIEMONTE[,12][nobs ])


 
 		newdata 	<-  x[ which(codice_regione=="2"),]
     	 	ValleAosta  <-  xts( newdata , order.by=dates )

   		newdata 	<- x[ which(codice_regione=="3"),]
      	Lombardia  <-  xts( newdata , order.by=dates )

   		newdata 	<- x[ which(codice_regione=="4"), ]
      	Trento   <-  xts( newdata[seq(2, 36, by = 2) ,] , order.by=dates )

   		newdata 	<- x[ which(codice_regione=="4"), ]
		Bolzano  <-  xts( newdata[seq(1, 36, by = 2) ,] , order.by=dates )


  		newdata 	<- x[ which(codice_regione=="5"),]
      	Veneto       <-  xts( newdata , order.by=dates )

	  	newdata 	 <- x[ which(codice_regione=="6"),]
      	Friuli       <-  xts( newdata , order.by=dates )

	  	newdata 	<-  x[ which(codice_regione=="7"),]
      	Liguria     <-  xts( newdata , order.by=dates )

	  	newdata 	<- x[ which(codice_regione=="8"),]
      	Emilia       <-  xts( newdata , order.by=dates )

  		newdata 	<- x[ which(codice_regione=="9"),]
      	Toscana  	<-  xts( newdata , order.by=dates )

		newdata 	<- x[ which(codice_regione=="10"),]
      	Umbria  	<-  xts( newdata , order.by=dates )

		newdata 	<- x[ which(codice_regione=="11"),]
      	Marche  	<-  xts( newdata , order.by=dates )

		newdata 	<- x[ which(codice_regione=="12"),]
      	Lazio  	<-  xts( newdata , order.by=dates )

		newdata 	<- x[ which(codice_regione=="13"),]
      	Abbruzzo  	<-  xts( newdata , order.by=dates )

		newdata 	<- x[ which(codice_regione=="14"),]
      	Molise  	<-  xts( newdata , order.by=dates )

		newdata 	<- x[ which(codice_regione=="15"),]
      	Campania  	<-  xts( newdata , order.by=dates )

		newdata 	<- x[ which(codice_regione=="16"),]
      	Puglia  	<-  xts( newdata , order.by=dates )


		newdata 	<- x[ which(codice_regione=="17"),]
      	Basilicata 	<-  xts( newdata , order.by=dates )

		newdata 	<- x[ which(codice_regione=="18"),]
      	Calabria 	<-  xts( newdata , order.by=dates )

		newdata 	<- x[ which(codice_regione=="19"),]
      	Sicilia 	<-  xts( newdata , order.by=dates )

		newdata 	<- x[ which(codice_regione=="20"),]
      	Sardegna 	<-  xts( newdata , order.by=dates )

## Graph
aa <- cbind(PIEMONTE [,12], Lombardia[,12], Veneto  [,12],   Liguria [,12], Emilia [,12] )

  plot.xts(aa, col=2:6, ylim= c(0,.4)  ,main="")
  addLegend(legend.loc = "topleft", legend.names = c("PIEMONTE","Lombardia","Veneto","Liguria", "Emilia Romagna") ,
  col =  c(2,3,4,5,6) , lwd = 2)




## binomial

  probBIASED <-  PIEMONTE [,12][16]
BIN <- matrix (NA,30,1)

 for (j in 1:30) {
 BIN [j]  <-   dbinom(j, size= 100, prob=probBIASED   ) 
			}
 max (BIN)


 ts.plot (BIN )
 

BIN_I <- matrix (NA,20,1)

 for (j in 1:20) {
BIN_I [j]  <-   dbinom(j, size= 100, prob= .03) 
			}
 ts.plot (BIN_I )


######3  transfer entropy

 x <-    Lombardia [,12]
 y <-    Veneto    [,12]

transfer_entropy(x, y,lx = 2, ly = 2, nboot = 555,quantiles = c(1,99))


 x <-    Lombardia [,12]
 y <-    PIEMONTE  [,12]

transfer_entropy(x, y,lx = 1, ly = 1, nboot = 555,quantiles = c(5,95))

 x <-    Lombardia [,12]
 y <-    Campania [,12]

transfer_entropy(x, y,lx = 1, ly = 1, nboot = 555,quantiles = c(5,95))


 x <-    Lombardia [,12]
 y <-    Lazio     [,12]

transfer_entropy(x, y,lx = 1, ly = 1, nboot = 555,quantiles = c(5,95))



 x <-    Lombardia [,12]
 y <-    Toscana [,12]

transfer_entropy(x, y,lx = 1, ly = 1, nboot = 555,quantiles = c(5,95))

 x <-    Lombardia [,12]
 y <-    Emilia       [,12]

transfer_entropy(x, y,lx = 1, ly = 1, nboot = 555,quantiles = c(5,95))

###### 4  Pueyo
 
 names(Lombardia  )
 
 M    <-  as.vector (    Lombardia     [nobs, 9]  )
 MvC  <-  as.vector (    Lombardia     [nobs, 13]  )   # 8% (morti/casi)
 K    <- 8.1
 Lombardia_Contagiati = M* (100/MvC) *8