arima_model.Rmd

---
params:
  series: '2hrs ph3'
  model: Arima(order=c(1, 0, 0), seasonal=c(1, 0, 0))
  #model: Arima(order=c(1, 0, 0), xreg=fourier(., K=2))
  #file:///D:/Personal/Dropbox/Facultate/ARIMA/v2/2hrs_ph3/Arima(order=c(1,%200,%200),%20seasonal=c(1,%200,%200),%20method=CSS,%20xreg=%7B%20%20%20%20cbind(dummies%20=%20getNthObsDummies(--d7ae590d5a60d061cd092b908e88f148.html
  transformation: identity() #log() %>% (function(d) {d[d==-Inf]<-0; d})
  diffs: identity() # diff() %>% diff()
  sdiffs: identity() # diff(lag=frequency(data$train))
  startday: -10
  traindays: 7
  testdays: 3
  xreg: NULL
  #xreg: fourier(., K=2)
  #
  obs: FALSE
title: '`r params$model` on `r params$series`'
author: 'Barbu Paul - Gheorghe'
output:
  html_document:
    toc: true
    toc_float: true
    number_sections: true
    self_contained: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r includes, warning=FALSE, include=FALSE}
library(forecast)
library(ggplot2)
library(scales)
library(tictoc)

source('datasets.R')
source('funcs.R')
```

# Parameters

Series: **`r params$series`**.

Model: **`r params$model`**.

Transformation: **`r params$transformation`**.

Differences: **`r params$diffs`**.

Seasonal differences: **`r params$sdiffs`**.

As observations: **`r params$obs`**.

Starting from: **`r params$startday`**.

Train days: **`r params$traindays`**.

Test days: **`r params$testdays`**.

xreg: **`r params$xreg`**.

# Original data

## Training & test data
The data has been previously cleaned, negative values were made 0.
Time series 80/20 rule goes in sequential order, not random.

```{r data, echo=FALSE, message=FALSE, warning=FALSE}
if(params$obs)
{
  data <- get_obs(datasets[[params$series]]$series,
                 params$startday,
                 params$traindays, params$testdays)
} else {
  data <- get_days(datasets[[params$series]]$series,
                   params$startday,
                   params$traindays, params$testdays)
}
# 
# autoplot(ts.union(data$test, data$train), facets = FALSE) +
#  aes(linetype = series, color= series) +
#  scale_linetype_manual(labels = c("Actual", "Forecasted"), values = c(2, 1)) #+
# 
# autoplot(ts.union(data$test, data$train)) +
#  aes(linetype = series, color= series) +
#  scale_linetype_manual(values = c(2, 1), name="Data series", labels=c("Test", "Train")) +
#   scale_color_manual(values = c(2, 1), name="Data series", labels=c("Test", "Train")) +
#  xlab("Observation [days]") +
#   ylab("Energy [Watts]") +
#   ggtitle("Original data") #+ 
#   #guides(colour=guide_legend(title="Data series")) 

lbls <- c("Test", "Train")
autoplot(ts.union(data$test, data$train)) +
  geom_point(aes(shape=series)) +
  scale_shape_discrete(labels=lbls) +
  scale_color_discrete(labels=lbls) +
  scale_x_continuous(breaks=pretty_breaks()) + 
  labs(shape="Data series", color="Data series") +
  xlab("Observation [days]") +
  ylab("Energy [Watts]") +
  ggtitle("Original data")
# 
# test_layer <- autolayer(data$test, series="Test") 
# test_values <- test_layer$layer_data()
# 
# autoplot(data$train, series="Train") +
#   test_layer +
#   geom_point() +
#   geom_point(data = test_values, aes(x = timeVal, y = seriesVal), shape=17) + 
#   xlab("Observation [days]") +
#   ylab("Energy [Watts]") +
#   
#   ggtitle("Original data") + 
#   #guides(colour=guide_legend(title="Data series"))
```

Total data points: **`r length(datasets[[params$series]]$series)`** representing
**`r round(length(datasets[[params$series]]$series)/frequency(data$train))`** days.

Number of data points per day: **`r frequency(data$train)`**
(gathered once every **`r (24*60)/frequency(data$train)`** minutes).

```{r subseries-plots, echo=FALSE, message=FALSE}

# TODO: try catch this
if(!params$obs)
{
season.labels = seq(frequency(data$train))
# ggseasonplot(data$train, season.labels = season.labels) +
#   xlab("Season") +
#   ylab("Energy [Watts]") +
#   guides(colour=guide_legend(title="Day")) +
#   ggtitle('Seasonal plot for training data') 

ggseasonplot(data$train, season.labels = season.labels) +
  xlab("Season") +
  ylab("Energy [Watts]") +
  ggtitle('Seasonal plot for training data') +
  geom_point(aes(shape=year)) +
  labs(color="Day", shape="Day")
}

if(!params$obs)
{
  ggsubseriesplot(data$train, labels = season.labels) +
    xlab("Season") +
    ylab("Energy [Watts]") +
    ggtitle('Subseries plot for each season in the training data') +
    scale_x_continuous(breaks = 0.5 + (1:frequency(data$train)),
                       labels = 1:frequency(data$train))
}
```

## ACF & PACF of training data

```{r acf-pacf, echo=FALSE}
data$train %>% ggtsdisplay(main='Training data', xlab='Observation [days]', ylab='Energy [Watts]')
```

#Data transformations & differences

## Box-Cox transformation

* Original data points: $y_{1},\dots,y_{T}$
* Transformed data points: $w_{1},\dots,w_{T}$

$$
w_t  =
    \begin{cases}
      \log(y_t)               & \text{if $\lambda=0$};  \\
      (y_t^\lambda-1)/\lambda & \text{otherwise}.
    \end{cases}
$$  

```{r boxcox, echo=FALSE}
lambda <- BoxCox.lambda(data$train)
```

Suggested BoxCox transformation on the training data:
$\mathbf{\lambda=`r lambda`}$

```{r boxcox-applied, echo=FALSE}
autoplot(BoxCox(data$train, lambda)) +
  xlab("Observation [days]") +
  ylab("Transformed energy") +
  ggtitle(bquote("Transformed training data, " ~ lambda ~ " = " ~ .(lambda)))
```

```{r apply-transformation, echo=FALSE}
use_transformation <- !is.null(params$transformation) && params$transformation != 'identity()'

eval(parse(text=paste('data$train %>%',
                      params$transformation))) -> data$train.transformed
```

## Transformation

Transformation: **`r params$transformation`**

```{r plot-transformed-data, echo=FALSE, eval=use_transformation}
autoplot(data$train.transformed) +
  xlab("Observation [days]") +
  ylab("Transformed energy") +
  ggtitle(paste("Transformed data using:", params$transformation)) + 
  guides(colour=guide_legend(title="Data series"))
```

**ATTENTION:** continuing with **`r if(use_transformation) 'transformed' else 'original'`** data!

## ACF & PACF after transforming & differencing the data

Suggested number of differences to make data stationary:
**`r ndiffs(data$train.transformed)`**.

Suggested number of seasonal differences to make data stationary:
**`r nsdiffs(data$train.transformed)`**.

```{r apply-diffs, echo=FALSE}
use_diffs <- !is.null(params$diffs) && params$diffs != 'identity()'

eval(parse(text=paste('data$train.transformed %>%',
                      params$diffs))) -> data$train.transformed.diffed
```

```{r apply-sdiffs, echo=FALSE}
use_sdiffs <- !is.null(params$sdiffs) && params$sdiffs != 'identity()'

eval(parse(text=paste('data$train.transformed.diffed %>%',
                      params$sdiffs))) -> data$train.transformed.sdiffed
```

### Differences

Differences: **`r params$diffs`**.

```{r plot-diff, echo=FALSE, eval=use_diffs}
data$train.transformed.diffed %>%
  ggtsdisplay(main=paste('Transformed & differenced data using:', params$diffs),
              xlab='Observation [days]',
              ylab='Transformed & differenced energy')
```

### Seasonal differences

Seasonal differences: **`r params$sdiffs`**.

```{r plot-sdiff, echo=FALSE, eval=use_sdiffs}
data$train.transformed.sdiffed %>%
  ggtsdisplay(main=paste('Transformed & diffed & s. diffed data using:',
                         params$diffs, 'and', params$sdiffs), 
              xlab='Observation [days]', 
              ylab='Transformed & diffed & s. diffed energy')
```

# Model fitting & evaluation

```{r model-fitting, include=FALSE}
tic("Model fitting")

eval(parse(text=paste('data$train.transformed %>%', params$model))) -> fit

times <- toc()
```

The model **`r params$model`** was evaluated on **`r params$traindays`** days
(representing **`r frequency(data$train.transformed) * params$traindays`** data points)
in **`r times$toc - times$tic`** seconds
(**`r (times$toc - times$tic)/60`** minutes).

## Model summary

```{r summary-fit}
summary(fit)
```

## Residuals statistics

```{r include=FALSE}
fit %>% residuals() -> data$residuals
```


```{r residuals-stats, echo=FALSE}
checkresiduals(fit,
               xlab='Observation [days]',
               ylab='Energy [Watts]') 
```

Box-Ljung test:

* H~0~ states that the data points are independently distributed.
* High p-value (> 0.05) doesn't reject H~0~.
* If Q\* is large: the autocorrelations do not come from a white noise series
(so Q\* should be small)

```{r residuals-pacf, echo=FALSE}
ggPacf(data$residuals) +
  ggtitle('PACF of residuals')
```

```{r qq-plot, echo=FALSE, message=FALSE}
ggplot(data.frame(y=data$residuals), aes(sample=y)) +
  stat_qq() +
  stat_qq_line() +
  ggtitle('Residuals QQ-plot')
```

## Forecasting
```{r forecast, echo=FALSE, include=FALSE}
tic('Forecast')

fourier.terms <- NULL

if(params$obs)
{
  h <- params$testdays
}else {
  h <- params$testdays * frequency(data$test)
}

if(!is.null(params$xreg))
{
  eval(parse(text=paste('data$test %>%', params$xreg))) -> fourier.terms
}

fit %>% forecast(h=h, 
                 xreg=fourier.terms) -> fcast

fcast.times <- toc()
```

Time elapapsed for forecasting **`r params$testdays`**
`r ifelse(params$obs, "data points", paste("days (representing", frequency(data$test) * params$testdays, "data points)"))`:
**`r fcast.times$toc - fcast.times$tic`** seconds
(**`r (fcast.times$toc - fcast.times$tic)/60`** minutes).

### Accuracy of the forecasts against the test data
```{r accuracy, echo=FALSE}
knitr::kable(accuracy(fcast, data$test))
#knitr::kable(accuracy(fcast, data$test), format='latex')
```

### Fitted values and forecasts plot against the real data
```{r fit-forecast-plot, echo=FALSE}
numtotaldays <- round(length(datasets[[params$series]]$series)/frequency(data$train))

startday <- params$startday
if(params$startday < 0) 
{
  startday <- end(datasets[[params$series]]$series)[1] + 1 + params$startday
}
# 
# clrs <- c("blueviolet", "blue", "darkgoldenrod4", "red")
# 
# plt <- autoplot(fcast) +
#   autolayer(fcast$mean, series="Forecast") +
#   autolayer(fitted(fit), series='Fitted') + 
#   autolayer(data$test, series='Test') +
#   xlab("Observation [days]") +
#   ylab("Energy [Watts]") +
#   guides(colour=guide_legend(title="Data series"), 
#          fill=guide_legend(title="Prediction interval")) +
#   scale_color_manual(values=clrs) +
#   scale_x_continuous(breaks=seq(from=startday, to=startday + params$traindays + params$testdays, by=2))
# 
# if(use_transformation)
# {
#   plt + autolayer(data$train.transformed, series = 'Train Transformed') 
# }else{
#   plt + autolayer(data$train, series = 'Train') 
# }


res_lbls <- c("Train", "Forecast", "Fitted", "Test")

# autoplot(ts.union(data$train, fcast$mean, fitted(fit), data$test)) +
#  aes(shape=series, linetype = series, color= series) +
#  scale_linetype_manual(values = c(12, 6, 2, 1), name="Data series", labels=lbls) +
#   scale_color_manual(values = clrs, name="Data series", labels=lbls) +
#  xlab("Observation [days]") +
#   ylab("Energy [Watts]") +
#   ggtitle(fcast$method)

autoplot(ts.union(data$train, fcast$mean, fitted(fit), data$test)) +
  geom_point(aes(shape=series)) +
  labs(color="Data series", shape="Data series") +
  scale_shape_discrete(labels=res_lbls) +
  scale_color_discrete(labels=res_lbls) +
  scale_x_continuous(breaks= pretty_breaks()) +
  xlab("Observation [days]") +
  ylab("Energy [Watts]") +
  ggtitle(fcast$method)


#https://stackoverflow.com/questions/49994580/autoplot-function-behaves-differently-when-i-load-ggfortify

#for_plot <- ggfortify:::fortify.forecast(fcast, ts.connect = TRUE)#, is.date=FALSE)
#for_plot$Index <- as.numeric(format(for_plot$Index,'%Y%m%d'))
#for_plot <- dplyr::rename(for_plot, Data="Train")
#for_plot <- ggfortify:::fortify.forecast(fcast, ts.connect = TRUE)
#
# ggplot(data = for_plot) +
#   geom_line(aes(x= Index, y = Data, color = "raw")) +
#   geom_point(aes(x= Index, y = Data, shape="raw", color="raw")) +
#   geom_line(aes(x= Index, y = Fitted, color = "fitted")) +
#   geom_line(aes(x= Index, y = `Point Forecast`, color = "point forecast")) +
#   geom_ribbon(aes(x= Index, ymin = `Lo 80`, ymax = `Hi 80`,  fill = "80"),  alpha = 0.2) +
#   geom_ribbon(aes(x= Index, ymin = `Lo 95`, ymax = `Hi 95`,  fill = "95"),  alpha = 0.2) +
#   #scale_x_continuous(breaks= pretty_breaks()) +
#   #scale_fill_manual("what", values = c("blue", "dodgerblue"))+
#   #scale_color_manual("why", values = c("blue", "red", "green"))
#   scale_fill_discrete("what") +
#   scale_color_discrete("why") +
#   scale_shape_discrete("why")
#
# for_plot_lines <- for_plot %>%
#   gather(key, value, 2:4) %>%
#   select(key, value, Index)



# for_plot %>%
#   filter(!is.na(`Point Forecast`)) %>%
#   gather(Lo, ymin, `Lo 95`, `Lo 80`) %>%
#   gather(Hi, ymax, `Hi 95`, `Hi 80`) -> for_plot_ribbon
# 
# for_plot %>%
#   filter(!is.na(`Point Forecast`)) %>%
#   gather(Lo, ymin, c( "Lo 95")) %>%
#   gather(Hi, ymax, c( "Hi 95")) -> for_plot_ribbon95
# 
# for_plot %>%
#   filter(!is.na(`Point Forecast`)) %>%
#   gather(Lo, ymin, c("Lo 80")) %>%
#   gather(Hi, ymax, c("Hi 80")) -> for_plot_ribbon80
# 
# ggplot(data = for_plot_lines) +
#   geom_line(aes(x= Index, y = value, color = key)) +
#   geom_point(aes(x= Index, y = value, shape = key, color=key)) +
#   geom_ribbon(data = for_plot_ribbon80, aes(x= Index, ymin = ymin, ymax = ymax, fill = Hi), alpha = 0.2)+
#   geom_ribbon(data = for_plot_ribbon95, aes(x= Index, ymin = ymin, ymax = ymax, fill = Hi), alpha = 0.2)+
#   scale_shape_discrete(labels=res_lbls) +
#   scale_color_discrete(labels=res_lbls) +
#   scale_fill_discrete(labels=c(80, 95)) +
#   labs(color="Data series", shape="Data series", fill="Prediction interval") +
#   scale_x_continuous(breaks=c(142, 150)) +
#   xlab("Observation [days]") +
#   ylab("Energy [Watts]") +
#   ggtitle(fcast$method)

```

Please note that the negative prediction intervals are due to the way these are
computed: the point forecast plus/minus the standard deviation multiplied 
by the 95% or 80% critical value of the two tailed Z distribution,
this usually means the point forecast plus/minus twice the standard error. 
([ref. page 10](https://people.duke.edu/~rnau/Notes_on_linear_regression_analysis--Robert_Nau.pdf)).

The lower part of the interval should be limited at zero during the
interpretation of the
results, but it is not on the plot due to symmetry and graphical reasons.



# Future work

* Add other statistical tests for stationarity, differencing, etc.
* Anderson-Darling statistical test for normality of the residuals reported along the QQ-plot
* Maybe it would be helpful to generate the (standardised) Residuals vs fitted plot, although it wouldn't be of much help since the residuals are plotted against time and their patterns are extracted and accounted for using the ACF & PACF plots, also the Ljung-Box test for independence is computed

# References
* [Forecasting: Principles and Practice, Rob J Hyndman and George Athanasopoulos, Monash University, Australia](https://otexts.org/fpp2/)
* [Statistical forecasting: notes on regression and time series analysis, Robert Nau, Fuqua School of Business, Duke University](http://people.duke.edu/~rnau/411home.htm)
* [Time Series Forecasting: Creating a seasonal ARIMA model using Python and Statsmodel, Sean Abu](http://www.seanabu.com/2016/03/22/time-series-seasonal-ARIMA-model-in-python/)
* [How to Create an ARIMA Model for Time Series Forecasting in Python, Jason Brownlee](https://machinelearningmastery.com/arima-for-time-series-forecasting-with-python/)
* [ Introduction to Forecasting with ARIMA in R, Ruslana Dalinina](https://www.datascience.com/blog/introduction-to-forecasting-with-arima-in-r-learn-data-science-tutorials)
* [Seasonal periods, Rob J Hyndman](https://robjhyndman.com/hyndsight/seasonal-periods/)
* [Math Notation for R Plot Titles: expression and bquote, Tyler Rinker](https://www.r-bloggers.com/math-notation-for-r-plot-titles-expression-and-bquote/)