accuracy - VAR with 6 variables

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:

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)

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)

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

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)

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(ir.1)

Calculate the IRF

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

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(ir.3)

Calculate impulse response

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

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(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(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")

Evaluating forecast accuracy

accuracy(var.aic.for, df_test)

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.