Data<- read.csv("D:/Salvador Ake/Trabajo/Artículos/Biocombustible/Mon_Data_R.csv") #View(Data) Data1<- Data[c(1,2,3,4,6,7,8)] #View(Data1) (Does not include RGDPW Real Gross Domestic Product World) #View(Data2) (Does not include any variable beyond 8/31/1998, is the full set of variables, already in time series mode) Data2<- ts(Data[c(80:295),c(2:8)] , start=c(1998, 8), end=c(2016, 7), frequency=12) pdf("D:/Salvador Ake/Trabajo/Artículos/Biocombustible/Original_data_Plot.pdf") plot(Data2) dev.off() #Begins function for unit root tests. If failed, it takes yields to the time series (except for rates, there it takes differences) .................................................................................................................. library(tseries) KPSS_for_sets<- function(Data_a, conditions_a, dir_a, name_table){ counter1 = 1 Data3<- data.frame(1:(nrow(Data_a)-1)) names_cols<-c("Obs") dir<- paste(dir_a, name_table,".csv") KPSS_tests<- c("H0","Level","Intercept","P$Value","Result") for(counter1 in 1:ncol(Data_a)){ variable <- Data_a[,counter1] conditions = c(conditions_a[counter1]) prueba<- kpss.test(variable, null = conditions,lshort = FALSE) if(prueba$p.value > .05) { result<- c(paste(colnames(Data_a)[counter1], "is stationary", sep = " "), conditions , prueba$p.value, "Do Not Reject H0" ) Data3<- cbind(Data3,variable[-1]) names_cols<- cbind(names_cols,colnames(Data_a)[counter1]) } else { result<- c(paste(colnames(Data_a)[counter1], "is stationary", sep = " "), conditions , prueba$p.value, "Reject H0" ) if(substr(colnames(Data_a)[counter1],1,4)=="RATE"){ assign(paste("dif",colnames(Data_a)[counter1],sep="_"),diff(as.numeric(Data_a[,counter1])) ) Data3<- cbind(Data3,get(paste("dif",colnames(Data_a)[counter1],sep="_"))) names_cols<- cbind(names_cols,(paste("dif",colnames(Data_a)[counter1],sep="_"))) }else{ assign(paste("Yield",colnames(Data_a)[counter1],sep="_"),((as.numeric(Data_a[,counter1][-1])/as.numeric(Data_a[,counter1])[-(length(Data_a[,counter1]))])-1)) Data3<- cbind(Data3,get(paste("Yield",colnames(Data_a)[counter1],sep="_"))) names_cols<- cbind(names_cols,(paste("Yield",colnames(Data_a)[counter1],sep="_"))) } } KPSS_tests<- cbind(KPSS_tests,result) } colnames(Data3)<- names_cols write.table(KPSS_tests, file = dir, sep = ",", col.names = NA, qmethod = "double") return(Data3) } #End function for unit root test.................................................................................................................. #View(Data3) Data_N<- KPSS_for_sets(Data2, c("Trend","Trend","Trend","Trend","Trend","Trend","Trend"), "D:/Salvador Ake/Trabajo/Artículos/Biocombustible/", "KPSS_test_Levels") Data_N<- KPSS_for_sets(Data_N,c("Trend","Level","Level","Level","Trend","Trend","Level"), "D:/Salvador Ake/Trabajo/Artículos/Biocombustible/", "KPSS_test_Yield") Data3<- ts(Data_N[,-1] , start=c(1998, 9), end=c(2016, 7), frequency=12) pdf("D:/Salvador Ake/Trabajo/Artículos/Biocombustible/Yield_data_Plot.pdf") plot(Data3) dev.off() #Initiates Vector Auto Regressive procedures.............................................................................................................................. library(vars) model1_data<- cbind(Data3[,1],Data3[,2],Data3[,3]) colnames(model1_data)<- c("Yield_FOI","Yield_FEI","Yield_FI") VAR_model1<- VAR(model1_data, p = 1, type = "none") VAR_model1_irf <- irf(VAR_model1, impulse = "Yield_FI", response = c("Yield_FOI", "Yield_FEI"), boot = FALSE) VAR_model1_stab <- stability(VAR_model1, type = "OLS-CUSUM") VAR_model1_norm_test<- normality.test(VAR_model1, multivariate.only = TRUE) VAR_model1_FVED <- fevd(VAR_model1, n.ahead = 5) VAR_model1_ARCH_test<- arch.test(VAR_model1, lags.single = 16, lags.multi = 5, multivariate.only = TRUE) #Ends Vector Auto Regressive procedures.............................................................................................................................. #Initiates de output record in a Word archive.......................................................................................................................... capture.output( summary(VAR_model1) ,causality(VAR_model1, cause = "Yield_FI", boot=TRUE, boot.runs=1000) ,print("Do not reject causality form Food Index to the VAR1 model") ,print("Reject instante causality form Food Index to the VAR1 model") ,arch.test(VAR_model1, lags.single = 16, lags.multi = 5, multivariate.only = TRUE) ,print("Do not reject ARCH effects in VAR1 model") ,normality.test(VAR_model1, multivariate.only = TRUE) ,print("Reject Multivariate normality for the residual of the VAR1 model") ,print(roots(VAR_model1, modulus = TRUE)) ,print("We show theeigenvalues of the VAR processdata(Canada) for the model") , file="D:/Salvador Ake/Trabajo/Artículos/Biocombustible/VAR_Model.doc") #library(ReporteRS) #Reporte_VAR <- docx(title = "VAR_model1") #Reporte_VAR <- addParagraph( Reporte_VAR, summary(VAR_model1), stylename = "rRawOutput" ) #writeDoc( Reporte_VAR, file = "D:/Salvador Ake/Trabajo/Artículos/Biocombustible/VAR_Report2.doc") #Ends de output record in a Word archive.......................................................................................................................... #Initiates the recording of the plots............................................................................................................................. pdf("D:/Salvador Ake/Trabajo/Artículos/Biocombustible/VAR_Report_Graphs.pdf") plot(VAR_model1_irf) plot(VAR_model1_stab) plot(VAR_model1_FVED) plot(VAR_model1) Hit Hit plot(VAR_model1_norm_test) Hit Hit plot(VAR_model1_ARCH_test) dev.off() #Ends the recording of the plots............................................................................................................................. library(rugarch) library(rmgarch) Data4<- Data3[,-c(3,6,7)] Data4<- ts(Data4 , start=c(1998, 9), end=c(2016, 7), frequency=12) # univariate normal GARCH(1,1) for each series garch1.spec <-ugarchspec(mean.model = list(armaOrder = c(1,0),include.mean=FALSE) , variance.model = list(garchOrder = c(1,1), model = "eGARCH"), distribution.model = "norm") # dcc specification - GARCH(1,1) for conditional correlations dcc.garch1.spec <- dccspec(uspec = multispec( replicate(ncol(Data4), garch1.spec) ), dccOrder = c(1,1), distribution = "mvt") dcc.fit_1 <- dccfit(dcc.garch1.spec, data = Data4) dcc.fit_1 # univariate normal GARCH(1,1) for each series garch1.spec2 <-ugarchspec(mean.model = list(armaOrder = c(1,0),include.mean=FALSE) , variance.model = list(garchOrder = c(1,1), model = "eGARCH"), distribution.model = "sged") # dcc specification - GARCH(1,1) for conditional correlations dcc.garch1.spec2 <- dccspec(uspec = multispec( replicate(ncol(Data4), garch1.spec2) ), dccOrder = c(1,1), distribution = "mvt") dcc.fit_2 <- dccfit(dcc.garch1.spec2, data = Data4) dcc.fit_2 Volatilities<-as.data.frame(dcc.fit_2@model$sigma) Volatilities2<- ts(Volatilidades, start=c(1998, 9), end=c(2016, 7), frequency=12) colnames(Volatilities2)<- colnames(Data4) pdf("D:/Salvador Ake/Trabajo/Artículos/Biocombustible/DCC_Graphs.pdf") plot(Volatilities2) plot(dcc.fit_2) dev.off() dev.off() slotNames(dcc.fit_2) names(dcc.fit_2@mfit) write.csv(dcc.fit_2@mfit$Qbar, "D:/Salvador Ake/Trabajo/Artículos/Biocombustible/Spillover_var_Biofuel.csv", row.names=TRUE) write.csv(dcc.fit_2@mfit$cvar, "D:/Salvador Ake/Trabajo/Artículos/Biocombustible/Var_covar_parameters_Biofuel.csv", row.names=TRUE) write.csv(dcc.fit_2@mfit$stdresid, "D:/Salvador Ake/Trabajo/Artículos/Biocombustible/residuals_Biofuel.csv", row.names=TRUE) write.csv(dcc.fit_2@model$sigma, "D:/Salvador Ake/Trabajo/Artículos/Biocombustible/volatilities_Biofuel.csv", row.names=TRUE) require(MVN) dcc.fit_2_norm_test<- hzTest((dcc.fit_2@mfit$stdresid),qqplot=FALSE) require(fNonlinear) dcc.fit_2_bds_test<- bdsTest((dcc.fit_2@mfit$stdresid), m = 3) Data_cor<- data.frame(1:(nrow(Data4)-0)) names_cols<-c("Obs") par(mfrow=c(1,1)) counter1 <- 1- for (counter1 in 1: ncol(Data4)){ counter2 <- 2 for(counter2 in 2: (ncol(Data4)-0) ){ if (counter2 <= counter1) { counter2 == (counter2 +1) } else{ counter2 == counter2 print(paste(counter1,"and", counter2, sep= "_")) print(counter2 <= counter1) assign(paste("Correlation",colnames(Data4)[counter1], "and", colnames(Data4)[counter2] ,sep="_"),as.vector((rcor(dcc.fit_2)[counter1,counter2,]))) Data_cor<- cbind(Data_cor, get(paste("Correlation",colnames(Data4)[counter1], "and", colnames(Data4)[counter2] ,sep="_"))) names_cols<- cbind(names_cols, paste(colnames(Data4)[counter1], "and", colnames(Data4)[counter2] ,sep="_")) } } } colnames(Data_cor)<- names_cols Data_cor<- Data_cor[,-1] require(xts) install.packages("xtsExtra", repos="http://R-Forge.R-project.org") require(PerformanceAnalytics) PerformanceAnalytics::chart.TimeSeries(as.xts(Data_cor)) Data_cor<- ts(Data_cor , start=c(1998, 8), end=c(2016, 7), frequency=12) pdf("D:/Salvador Ake/Trabajo/Artículos/Biocombustible/DCC_Correlations_Graphs.pdf") plot(Data_cor, cex = 0.8, cex.lab= 0.8, cex.main= 0.8, cex.sub= 0.8, cex.axis = 0.8, main = "Conditional Correlations from the DCC model") dev.off() #View(Data_cor)