# предсказания данных из
# https://petolau.github.io/Regression-trees-for-forecasting-time-series-in-R/
# только на сезонной arima без других пакетов

library(forecast)

DT <- read.csv("DT_load_17weeks.csv")

# Тренировка модели на 5712 строках, тест на послених 1000
trainIndexes <- 1:(nrow(DT)-1000)
testIndexes <- (nrow(DT)-1000+1):nrow(DT)
# На серьёзных данных не стоит обучать модель сразу на всех данных, лучше торговать каким-то скользящим окном,
# например обучать на 1000 точек, торговать следующих 100, потом обучить модель по новой, итд

# создание модели по ранее найденным параметрам, чтоб долго не ждать. Но с такими параметрами (q=24) всё равно будет долго.
arimaModel <- Arima(y = ts(DT$value[trainIndexes],frequency=48), order = c(1, 0, 24), seasonal = c(1, 0, 1))
# поиск параметров заново, если нужно
# arimaModel <- auto.arima(y = ts(DT$value[trainIndexes],frequency=48),
#                          seasonal.test = "ocsb",
#                          trace=TRUE,
#                          stepwise = FALSE,
#                          max.q = 48, 
#                          max.order = 48+5
#                         )


# можно взять и предсказать все 1000 будущих точек одним махом.
# prediction <- c(forecast(arimaModel,h=length(testIndexes))$mean)
# но мы так делать не будем. Модель в таком случае будет на каждой новой точке терять точность,
# т.к. не имеет досутпа к ошибкам своего предсказания

# поэтому предсказываем по одной точке, и после этого даём моделе информацию об ошибке, и предсказываем следующую точку.
# сами найденные ранее AR I MA коэффициенты заново искаться не будут, поэтому потеря производительности будет небольшой


# комбинируем два графика в одном окне. Для предсказания, и для ошибок
par(mfrow=c(2,1))
# истинные значения и предсказания

prediction <- rep(NA, length(testIndexes))
for(i in 1:length(testIndexes)){
  predictionIndex <- testIndexes[i]
  correctionIndexes <- (predictionIndex-48*7):(predictionIndex-1)
  arimaModelCorrected <- Arima(y = DT$value[correctionIndexes], model = arimaModel)
  prediction[i] <- forecast(arimaModelCorrected,h=1)$mean[1]
}

# комбинируем два графика в одном окне. Для предсказания, и для ошибок
par(mfrow=c(2,1))
# истинные значения и предсказания
plot(DT$value[testIndexes], type="l", col="blue", ylim=c(min(DT$value[testIndexes], prediction),max(DT$value[testIndexes], prediction)), main="Forecast & Value", ylab = "forecast")
lines(prediction, col="red")
grid()
legend(x = "bottomright", legend = c("value","forecast"), fill=c("blue","red"), bty="n")
# ошибки предсказания
plot(DT$value[testIndexes] - prediction, main="Forecast errors", ylab = "error")
grid()


summary(arimaModel)
cat("arima train r^2:", 1 - sum((residuals(arimaModel))^2)/sum((DT$value[trainIndexes]-mean(DT$value[trainIndexes]))^2), "\n")
cat("arima test r^2:", 1 - sum((DT$value[testIndexes] - prediction)^2)/sum((DT$value[testIndexes]-mean(DT$value[testIndexes]))^2), "\n")
cat("arima test diff r^2:", 1 - sum((diff(DT$value[testIndexes]) - diff(prediction))^2)/sum((diff(DT$value[testIndexes])-mean(diff(DT$value[testIndexes])))^2), "\n")
cat("arima test diff accuracy:", mean(diff(DT$value[testIndexes]) * diff(prediction) > 0), "\n")