Instalando os pacotes
Projetando sem alterar o VECM
### preparing data
rm(list = ls(all.names = TRUE))
gc()
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2736605 146.2 5347722 285.6 3556154 190
## Vcells 4589068 35.1 10146329 77.5 8388608 64
## importando dados
path <- "d:/Users/thomas.daher/Desktop/cimento/dados_parte1/"
file_name <- "Obras hidreletricas.xlsx"
sheet_name <- "serie_curta"
file_path <- paste(path, file_name, sep = "")
cimento <- read_excel(file_path, sheet = sheet_name)
cimento$Data <- as.Date(cimento$Data, format = "%Y/%m/%d")
cimento <- filter(cimento, Data >= as.Date("1990-09-01"))
cimento_17 <- rev(diff(ts((cimento$cimento17), start = c(1990, 9), frequency = 12)))
cimento_51 <- rev(diff(ts((cimento$cimento51), start = c(1990, 9), frequency = 12)))
sheet_name <- "serie_longa"
cimento_long <- read_excel(file_path, sheet = sheet_name)
cimento_long$Data <- as.Date(cimento_long$Data, format = "%Y/%m/%d")
cimento_17_long <- rev(diff(ts((cimento_long$cimento17[-1]), start = c(1979, 5), frequency = 12)))
dset <- cbind(cimento_17,cimento_51)
vecm_model <- VECM(
dset,
2,
r = 1,
include = c("none"),
beta = NULL,
estim = c("2OLS"),
LRinclude = c("none"),
exogen = NULL
)
summary(vecm_model)
## #############
## ###Model VECM
## #############
## Full sample size: 102 End sample size: 99
## Number of variables: 2 Number of estimated slope parameters 10
## AIC 4424.352 BIC 4452.898 SSR 2.024838e+12
## Cointegrating vector (estimated by 2OLS):
## cimento_17 cimento_51
## r1 1 -1.079909
##
##
## ECT cimento_17 -1 cimento_51 -1
## Equation cimento_17 -1.2910(0.4056)** 0.6807(0.3229)* -0.0149(0.3957)
## Equation cimento_51 0.0069(0.1959) 0.4480(0.1559)** -0.0531(0.1911)
## cimento_17 -2 cimento_51 -2
## Equation cimento_17 -0.0985(0.2271) -0.1862(0.1805)
## Equation cimento_51 -0.0130(0.1097) -0.0378(0.0872)
forecast_one_period <- function(coef_matrix, lagged_values1, lagged_values2) {
# Extract the cointegrating vector coefficients
cointegrating_vector <- coef_matrix[1,]
# Extract the ECT coefficients
ect_coefficients <- coef_matrix[2,]
# Extract the lag coefficients
lag_coefficients <- coef_matrix[3:nrow(coef_matrix), ]
# Calculate the error correction term (ECT)
ect <- (lagged_values1[2] - lagged_values1[1])
# Calculate the forecast for cimento_51 using lag 1 and lag 2
forecast_values <- cointegrating_vector[2] +
lagged_values1[2] * lag_coefficients[2, 2] +
lagged_values1[1] * lag_coefficients[1, 2] +
lagged_values2[1] * lag_coefficients[4, 2] +
lagged_values2[2] * lag_coefficients[3, 2] +
ect
# Return the forecasts
return(forecast_values)
}
# Create a dataframe to store forecasted values
forecast_51 <- data.frame(cimento_51 = numeric(0))
### predicting
exogenous <- cimento_17_long
num_periods <- (length(cimento_17_long) - length(cimento_51))
exogenous <- tail(exogenous, n = num_periods)
# matrix de coeficientes
coef_matrix <- matrix(c(1, -1.079909, -1.291004652, 0.006915527,0.6806797, 0.4480329, -0.01492402, -0.05309733, -0.09847807, -0.01298321, -0.18624985, -0.03776402), nrow = 6, byrow = TRUE)
# Extract the last observed values for cimento_17 and cimento_51
lagged_values1 <- as.vector(c(tail(cimento_17, 1), tail(cimento_51, 1)))
lagged_values2 <- as.vector(c(tail(cimento_17, 2)[1], tail(cimento_51, 2)[1]))
# Iterate over the desired number of periods
for (i in 1:num_periods) {
# Calculate the one-period-ahead forecast using the provided coefficients and lagged values
forecast_values <- forecast_one_period(coef_matrix, lagged_values1, lagged_values2)
# Append the forecasted values to the dataframe
forecast_51 <- rbind(forecast_51, forecast_values)
# Update lagged values for the next iteration
lagged_values2 <- lagged_values1
lagged_values1 <- as.vector(c(head(exogenous[1]), forecast_values))
exogenous <- exogenous[-1]
}
# Print the forecast dataframe
print(forecast_51)
## X5.21310637124409
## 1 5.213106
## 2 -5.287579
## 3 -12.100444
## 4 -13.087842
## 5 -13.322636
## 6 -13.497248
## 7 -14.792996
## 8 -31.558055
## 9 -39.890621
## 10 -46.117407
## 11 -47.795768
## 12 -47.900921
## 13 -46.905984
## 14 -45.545733
## 15 -44.426641
## 16 -42.945132
## 17 -41.324726
## 18 -39.682452
## 19 -38.126168
## 20 -36.668587
## 21 -35.407806
## 22 -34.275129
## 23 -33.171562
## 24 -32.134578
## 25 -31.122134
## 26 -30.177328
## 27 -29.281580
## 28 -28.436341
## 29 -27.644445
## 30 -26.910437
## 31 -26.215992
## 32 -25.563199
## 33 -24.957779
## 34 -24.389769
## 35 -23.854517
## 36 -23.354756
## 37 -22.887212
## 38 -22.451041
## 39 -22.041786
## 40 -21.660273
## 41 -21.307817
## 42 -20.982888
## 43 -20.673544
## 44 -20.385019
## 45 -20.115752
## 46 -19.863473
## 47 -19.627592
## 48 -19.407439
## 49 -19.202036
## 50 -19.010397
## 51 -18.831601
## 52 -18.664786
## 53 -18.509150
## 54 -18.363944
## 55 -18.228454
## 56 -18.102860
## 57 -17.985316
## 58 -17.876166
## 59 -17.773932
## 60 -17.678535
## 61 -17.589307
## 62 -17.506217
## 63 -17.428287
## 64 -17.355841
## 65 -17.287998
## 66 -17.224630
## 67 -17.165849
## 68 -17.110966
## 69 -17.059488
## 70 -17.011652
## 71 -16.967051
## 72 -16.925488
## 73 -16.886444
## 74 -16.850109
## 75 -16.816182
## 76 -16.784589
## 77 -16.755087
## 78 -16.727349
## 79 -16.701592
## 80 -16.677577
## 81 -16.655141
## 82 -16.634229
## 83 -16.614699
## 84 -16.596487
## 85 -16.579473
## 86 -16.563620
## 87 -16.548809
## 88 -16.535007
## 89 -16.522128
## 90 -16.510095
## 91 -16.498866
## 92 -16.488381
## 93 -16.478604
## 94 -16.469480
## 95 -16.460972
## 96 -16.453039
## 97 -16.445621
## 98 -16.438699
## 99 -16.432241
## 100 -16.426226
## 101 -16.420605
## 102 -16.415373
## 103 -16.410479
## 104 -16.405920
## 105 -16.401660
## 106 -16.397688
## 107 -16.393986
## 108 -16.390523
## 109 -16.387293
## 110 -16.384285
## 111 -16.381474
## 112 -16.378849
## 113 -16.376404
## 114 -16.374125
## 115 -16.371999
## 116 -16.370010
## 117 -16.368158
## 118 -16.366428
## 119 -16.364818
## 120 -16.363314
## 121 -16.361907
## 122 -16.360594
## 123 -16.359372
## 124 -16.358233
## 125 -16.357169
## 126 -16.356176
## 127 -16.355248
## 128 -16.354381
## 129 -16.353572
## 130 -16.352818
## 131 -16.352114
## 132 -16.351459
## 133 -16.350846
## 134 -16.350275
## 135 -16.349743
## 136 -16.349245
Analise do teste
# desrevertendo e transformando em série novamente
cimento_51_barra <- rev(forecast_51$X5.21310637124409) %>% cumsum()
plot(cimento_51_barra)
exogenous <- cimento_17_long
exogenous <- tail(exogenous, n = num_periods)
cimento_17 <- rev(exogenous) %>% cumsum()
df <- as.data.frame(cbind(cimento_17,cimento_51_barra))
ggplot(df, aes(x = 1:(136))) +
geom_line(aes(y = cimento_17, color = "cimento_17"), size = 1) +
geom_line(aes(y = cimento_51_barra, color = "cimento_51_barra"), size = 1) +
labs(title = "Séries plotando somente o projetado versus realizado cimento 17", y = "valor") +
scale_color_manual(values = c("cimento_17" = "blue", "cimento_51_barra" = "red"))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
cimento_completo <- c(cimento_51_barra, cimento$cimento51)
# Create a data frame
df <- data.frame(
time = 1:length(cimento_completo),
cimento_17 = cimento_long$cimento17[-1],
cimento_51_barra = cimento_completo
)
ggplot(data = df, aes(x = time)) +
geom_line(aes(y = cimento_17, color = "cimento_17"), size = 1) +
geom_line(aes(y = cimento_51_barra, color = "cimento_51_barra"), size = 1) +
labs(title = "Séries plotando somente o projetado versus realizado cimento 17", y = "valor") +
scale_color_manual(values = c("cimento_17" = "blue", "cimento_51_barra" = "red"))
