FEVD protocol

# version information: v0.2
# date of update: 23/Mar/2018

library(vars) # "vars" is the name of an R package used for vector autoregression (VAR) analysis

data_pre <- read.csv("filename.csv")  # "filename" is the name of the csv file
target   <- 3  # Number of target variables
term_num <- 12 # Number of terms in the data

# An example of auto-scaled data prepared for forecast error variance decomposition (FEVD) approach
# Names |Variable 1 |Variable 2 |Variable 3 | ...
# March |-0.4       | 0.5       |-0.5       | ...
# April |-0.3       | 0.6       | 2.8       | ...
# May   | 0.4       |-0.1       |-0.2       | ...
# ...   |...        |...        |...        | ...

# VAR (function)
var_v2 <- function(data_pre,target,term_num){

cur_dir <- getwd()
dir.create("var_v2")
setwd("var_v2")
data_pre <- as.matrix(data_pre)
target   <- as.integer(target)
term_num <- as.integer(term_num)

time_data <- data_pre[,1]##In this case, the first colum is month
time_data <- as.character(time_data)
select_ta <- data_pre[,2:(target+1)]##target microbe
fluc_data <- data_pre[,-1:-(target+1)]##taget fluctuation data

for(k in 1:length(select_ta[1,])){

result_1 <- matrix(NA,term_num * 2 - 1,length(fluc_data[1,]))
	rownames(result_1) <- c(seq(1:(term_num - 1)),"RMSE",seq(1:(term_num - 1)))
	result_2 <- matrix(NA,term_num + 3,length(fluc_data[1,]))
	rownames(result_2) <- c(time_data,"total","average","sd")

for(i in 1:length(fluc_data[1,])){
		var_1           <- cbind(select_ta[,k],fluc_data[,i])
		colnames(var_1) <- c(colnames(select_ta)[k],colnames(fluc_data)[i])
			
		var_ts<-as.matrix(ts(var_1))
		mode(var_ts) <- "numeric"
		var_selected <- VARselect(var_ts, lag.max=nrow(data_pre)-2,type="both")
		var_res      <- VAR(var_ts,p=var_selected$selection[1])
		
		var_res_part<-var_res$varresult
		predict_value  <- var_res_part[[1]]$fitted.values
		measured_value <- var_res_part[[1]]$model$y
		RMSE           <- sqrt((sum((measured_value-predict_value)^2))/length(measured_value))
		result_frag_1  <- c(predict_value,RMSE,measured_value)
		
		result_1[,i]   <- result_frag_1
		
		var_res_fevd  <- fevd(var_res,n.ahead = term_num)
		cont_val      <- var_res_fevd[[1]][,2]
		cont_val_sum  <- sum(cont_val)
		cont_val_mean <- mean(cont_val)
		cont_val_sd   <- sd(cont_val)
		result_frag_2 <- c(cont_val,cont_val_sum,cont_val_mean,cont_val_sd)
		
		result_2[,i] <- result_frag_2
		}
	colnames(result_1) <- colnames(fluc_data)
	colnames(result_2) <- colnames(fluc_data)
	
	write.table(result_1,paste(colnames(select_ta)[k],"_var_model.csv",sep = "")
                    ,col.names = colnames(result_1), row.names = rownames(result_1),sep = ",")
	write.table(result_2,paste(colnames(select_ta)[k],"_FEVD.csv",sep = "")
                    ,col.names = colnames(result_2), row.names = rownames(result_2),sep = ",")
	}
setwd(cur_dir)
}

# Get the results of VAR (useage)
var_v2(data_pre, target, term_num)