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$serialKacque Beta Test
norm1 <- normality.test(var.aic)
norm1$jb.mulARCH test
arch1 <- arch.test(var.aic, lags.multi = 6)
arch1$arch.mulplot(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)