accuracy - VAR with 6 variables

[max.resende]

Dear R Comunity,

I am having difficulties to do accuracy analysis after a VAR prediction.
I keep getting this error message:
First argument should be a forecast object or a time series.

I think it is because I am estimating a VAR with 6 variables. The one I need is energy_se230. I have tried to upload my data file, but I wasn't able to.I am new at this forum.

Link to data file:
https://drive.google.com/file/d/1Uzq9VVK_vd37E0wIG2gqU6FPEXGUWFzT/view?usp=sharing

Can anyone help me, please?

Thank you all,

Max

library(readxl)
prod <- read_excel("C:/Users/Max/Desktop/Anglo-American/Consumo_Energia/Sprint_2/prod_sprint2.xlsx")
View(prod)
dim(prod)
head(prod)

Define date as date:

as.Date(prod$Date, tryformats = c("%Y/%m/%d","%Y-%m-%d"), origin ="2018-12-23")
str(prod)

library(xts)

subsample$variables to plot

Demanda de energia diária - definida como time seriews

a_energy <- xts(prod$energy_se230, prod$Date)
plot(a_energy, main = 'Electricity Demand (MWh)', xlab = 'Date', ylab = 'Demand (MWh)')

############# BRITAGEM

a_britag_h_p <- xts(prod$Britag_HorasParadas, prod$Date)
plot(a_britag_h_p, main = 'Britagem Horas Paradas', xlab = 'Date', ylab = 'Horas')

a_britag_h_t <- xts(prod$Britag_HrTrab, prod$Date)
plot(a_britag_h_t, main = 'Britagem Horas Trabalhadas', xlab = 'Date', ylab = 'Horas')

a_britag_prod <- xts(prod$Britag_Prod, prod$Date)
plot(a_britag_prod, main = 'Britagem Produção', xlab = 'Date', ylab = 'Ore (ton)')

############# USINA

a_usin_h_p <- xts(prod$Usina_HorasParadas, prod$Date)
plot(a_usin_h_p, main = 'Usina Horas Paradas', xlab = 'Date', ylab = 'Horas')

a_usin_h_t <- xts(prod$Usin_HrTrab, prod$Date)
plot(a_usin_h_t, main = 'Usina Horas Trabalhadas', xlab = 'Date', ylab = 'Horas')

a_usin_alim <- xts(prod$Usin_AlimMoag, prod$Date)
plot(a_usin_alim, main = 'Usina Alimentação', xlab = 'Date', ylab = 'ore (ton)')

a_usin_prod <- xts(prod$Usin_Prod, prod$Date)
plot(a_usin_prod, main = 'Usina Produção', xlab = 'Date', ylab = 'ore (ton)')

############# PIPELINE

a_miner_h_p <- xts(prod$Mineroduto Horas Paradas, prod$Date)
plot(a_miner_h_p, main = 'Mineroduto Horas Paradas', xlab = 'Date', ylab = 'Horas')

a_miner_h_t <- xts(prod$Miner_HrTrab, prod$Date)
plot(a_miner_h_t, main = 'Mineroduto Horas Trabalhadas', xlab = 'Date', ylab = 'Horas')

a_miner_bomb <- xts(prod$Miner_Bomb, prod$Date)
plot(a_miner_bomb, main = 'Mineroduto Bombeamento', xlab = 'Date', ylab = 'Ore (ton)')

Merger multiple timeseries array into one dataframe

df_ts_energy <- data.frame(a_energy, a_britag_h_p, a_britag_h_t, a_britag_prod, a_usin_h_p, a_usin_h_t, a_usin_alim, a_usin_prod, a_miner_h_p, a_miner_h_t, a_miner_bomb)

#####################################################################################

Time Series dataframe

#####################################################################################
library(xts)

First Difference energy_se230

energy_se230 <- xts(prod$energy_se230, prod$Date)
diff_energy_se230 <- diff(energy_se230)
plot(diff_energy_se230, main = 'Differenced Daily Electricity demand', xlab = 'Date', ylab = 'Log Demand (MWh)')

First difference Britag_Prod

britag_prod <- xts(prod$Britag_Prod, prod$Date)
diff_britag_prod <- diff(britag_prod)
plot(diff_britag_prod, main = 'FDifferenced Daily Primary Crushing', xlab = 'Date', ylab = 'Ore (ton)')

First difference Usin_AlimMoag

usin_alim <- xts(prod$Usin_AlimMoag, prod$Date)
diff_usin_alim <- diff(usin_alim)
plot(diff_usin_alim, main = 'Differenced Daily Griding feed', xlab = 'Date', ylab = 'Ore (ton)')

First difference Usin_Prod

usin_prod <- xts(prod$Usin_Prod, prod$Date)
diff_usin_prod <- diff(usin_prod)
plot(diff_usin_prod, main = 'Differenced Daily Wet Route Prod', xlab = 'Date', ylab = 'Ore (ton)')

First difference Usin_HorasParadas

usin_h_p <- xts(prod$Usina_HorasParadas, prod$Date)
diff_usin_h_para <- diff(usin_h_p)
plot(diff_usin_h_para, main = 'Differenced Daily Work Downtime Wet Route', xlab = 'Date', ylab = 'Horas')

First difference Mineroduto horas paradas

mine_h_p <- xts(prod$Mineroduto Horas Paradas, prod$Date)
diff_mine_h_p <- diff(mine_h_p)
plot(diff_mine_h_p, main = 'Differenced Daily Work downtime Pipeline', xlab = 'Date', ylab = 'Horas')

Generate a time series_data frame

df_ts_prod <- data.frame(energy_se230, britag_prod, usin_alim, usin_prod, usin_h_p, mine_h_p)

Once selected variables ae non-stationary at level (ADf Test), we need todifferencied all ofthen.

diff_df_ts_prod <- diff(as.matrix(df_ts_prod))

Convert the data to axts object

df_forecast <- as.xts(diff_df_ts_prod)

Count the nu,ber of observations (days)

periodicity(df_forecast)
ndays(df_forecast)

###############################################################

SPLIT DATA frame for forecast: Proportin 80/20

Get all data until 2020 -01-01

df_train <- df_forecast["/2020-02"]

Get all data from 2020-01 until 2020 -04

df_test <- df_forecast["2020-03/2020-05"]

#############################################################

VAR MODEL

############################################################

Lag.lenght of VAR Model

library(vars)
varorder <- vars::VARselect(y = df_train, lag.max = 6, type = "const")
print(varorder$selection)

Estimate VAR model based on AIC Criteria

var.aic <- VAR(df_train , type = "none", lag.max = 6, ic = "AIC")
summary(var.aic)

Diagnostics of Fit for equation energy_se230

summary(var.aic, equation = "energy_se230")
plot(var.aic,names="energy_se230")

Diagnostic tests of VAR(p) specifications

ser11 <- serial.test(var.aic, lags.pt = 16, type = "PT.asymptotic")
ser11$serial

Kacque Beta Test

norm1 <- normality.test(var.aic)
norm1$jb.mul

ARCH test

arch1 <- arch.test(var.aic, lags.multi = 6)
arch1$arch.mul

plot(arch1, names = "energy_se230")
plot(stability(var.aic), nc = 2)

irf - CONTEMPORANEUS SHOCK

ARGUMENT ORTHO? Note that the ortho option is important, because it says something about the contemporaneous relationships between the variables.

If we know that such relationships do not exist, because the true variance-covariance matrix - or simply covariance matrix - is diagonal with zeros in the off-diagonal elements.

Calculate the IRF

ir.1 <- irf(var.aic, impulse = "britag_prod", response = c("energy_se230","usin_prod","usin_alim") , n.ahead = 20, ortho = FALSE)

Plot the IRF

plot(ir.1)

Calculate the IRF

ir.2 <- irf(var.aic, impulse = "britag_prod", response = "energy_se230", n.ahead = 20, ortho=TRUE, boot = TRUE)

Plot the IRF

plot(ir.2)

Calculate impulse response

ir.3 <- irf(var.aic, impulse = "usin_alim", response = "energy_se230", n.ahead = 20,ortho = TRUE, boot = TRUE)

Plot

plot(ir.3)

Calculate impulse response

ir.4 <- irf(var.aic, impulse = "usin_prod", response = "energy_se230", n.ahead = 20,rtho = TRUE, boot = TRUE)

Plot

plot(ir.4)

Calculate impulse response

ir.5 <- irf(var.aic, impulse = "usin_h_p", response = "energy_se230", n.ahead = 20,ortho = TRUE, boot = TRUE)

Plot

plot(ir.5)

Calculate impulse response

ir.6 <- irf(var.aic, impulse = "mine_h_p", response = "energy_se230", n.ahead = 20,ortho = TRUE, boot = TRUE)

Plot

plot(ir.6)

Estimate VAR model based on AIC Criteria

var.aic.for <- VAR(df_train , type = "none", lag.max = 6, ic = "AIC")
summary(var.aic.for)

prd <- predict(var.aic, equation = "energy_se230", n.ahead = 14, ci = 0.90, dumvar = NULL)
print(prd)

plot(prd, names="energy_se230")

Plot all equations at once

plot(prd, "single")

Serveral methods of forecasting will be implemented to measure which model would be the best for forecasting a horizon up to 14 days .

Evaluating forecast accuracy

accuracy(var.aic.for, df_test)