# Libraries ---------------------------------------------------------------

library('tidyverse')
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.8
## v tidyr   1.2.0     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library('lubridate')
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library('forecast')
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library('feasts')
## Loading required package: fabletools
## 
## Attaching package: 'fabletools'
## The following objects are masked from 'package:forecast':
## 
##     accuracy, forecast
library('fpp3') 
## -- Attaching packages -------------------------------------------- fpp3 0.4.0 --
## v tsibble     1.1.1     v fable       0.3.1
## v tsibbledata 0.4.0
## -- Conflicts ------------------------------------------------- fpp3_conflicts --
## x lubridate::date()      masks base::date()
## x dplyr::filter()        masks stats::filter()
## x fabletools::forecast() masks forecast::forecast()
## x tsibble::intersect()   masks base::intersect()
## x tsibble::interval()    masks lubridate::interval()
## x dplyr::lag()           masks stats::lag()
## x tsibble::setdiff()     masks base::setdiff()
## x tsibble::union()       masks base::union()
library('gridExtra') 
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

Loading the usual libraries for forecasting and modelling

# Data --------------------------------------------------------------------

aus_production[218,]
## # A tsibble: 1 x 7 [1Q]
##   Quarter  Beer Tobacco Bricks Cement Electricity   Gas
##     <qtr> <dbl>   <dbl>  <dbl>  <dbl>       <dbl> <dbl>
## 1 2010 Q2   374      NA     NA   2401       58041   236
beer <- ts(aus_production$Beer, frequency = 4, start = c(1956,1))
beer <- beer %>% as_tsibble()

Using Australian beer production because its fun and I was having an issue with finding a good data set.

# Exploratory graphing ----------------------------------------------------

beer %>% 
  autoplot() + 
  labs(title = "Australian Beer Production", 
       y = "Production", 
       x = "Date in Quarters")
## Plot variable not specified, automatically selected `.vars = value`

There is a lot of seasonality with this production.

# Modelling ---------------------------------------------------------------

#Building the model
fit <- beer %>% 
  model(NNETAR(value)) 
fit %>% report()
## Series: value 
## Model: NNAR(4,1,2)[4] 
## 
## Average of 20 networks, each of which is
## a 4-2-1 network with 13 weights
## options were - linear output units 
## 
## sigma^2 estimated as 260.7

I beleive the model is automatically built but not really specified in the text.

#Forecasting graph
fit %>% 
  forecast(h = 8) %>% 
autoplot(beer %>%
           filter(year(index) >= 1990)) +
  labs(title = "Australian Beer Production", 
       x = "Date in Quarters", 
       y = " Production")

I like how nueral nets allow asymmetric cyclical forecasting. This is unlike ARIMA models where everything has to be symmetrical.

fit %>% 
  forecast(h = 8) %>% 
  hilo()
## # A tsibble: 8 x 6 [1Q]
## # Key:       .model [1]
##   .model          index        value .mean                  `80%`
##   <chr>           <qtr>       <dist> <dbl>                 <hilo>
## 1 NNETAR(value) 2010 Q3 sample[5000]  422. [400.6312, 442.2002]80
## 2 NNETAR(value) 2010 Q4 sample[5000]  491. [470.0732, 511.4496]80
## 3 NNETAR(value) 2011 Q1 sample[5000]  418. [396.3948, 438.8591]80
## 4 NNETAR(value) 2011 Q2 sample[5000]  388. [367.1087, 409.4012]80
## 5 NNETAR(value) 2011 Q3 sample[5000]  426. [398.5424, 453.6001]80
## 6 NNETAR(value) 2011 Q4 sample[5000]  493. [464.7622, 520.0735]80
## 7 NNETAR(value) 2012 Q1 sample[5000]  422. [396.1115, 448.4765]80
## 8 NNETAR(value) 2012 Q2 sample[5000]  397. [371.6581, 422.5434]80
## # ... with 1 more variable: `95%` <hilo>

I used the hilo function to get the prediction intervals rather than the creating a simulation of the paths

fit %>%
  generate(times = 9, h = 30) %>%
  autoplot(.sim) +
  autolayer(beer, value)

I did also run the simulation of the paths. More information of how much value this brings would probably be beneficial