Forecasting ETS and ARIMA

#a. 

vic_pigs <- aus_livestock %>%
  filter(State == "Victoria", Animal == "Pigs")


vic_pigs %>%
  autoplot(Count) +
  theme_ipsum_es() + 
  labs(y = "Count"
         , x = "Date"
         , title = "Pigs Slaughtered, Victoria") +
  geom_smooth(formula = y ~ x, method = "loess")

model_pigs <- vic_pigs %>%
  model(ETS(Count~error("A") + trend("N") + season("N")))

report(model_pigs)

#Actual Model
pig_decomp <- vic_pigs %>%
  model(ETS(Count~error("A") + trend("N") + season("N"))) %>%
  components()

#Plot ETS decomposition
autoplot(pig_decomp) 

  # +
  # labs(title =
  #   "Decomposition of total pigs slghtrd using ETS")

#Forecast 4 months using "Model Pigs", using ETS model
fcst_pigs <- 
  model_pigs %>%
  forecast(h = 4)

#Plot Forecast
vic_pigs_plot <-
  vic_pigs %>%
  filter(yearmonth(Month) > yearmonth("2015 Dec")) %>%
  index_by(Date = as_date(Month))

fcst_pigs %>%
  autoplot(vic_pigs_plot) +
  labs(x = "Month",
       y = "Slaughter, Head",
       title = "Monthly Slaughter + Forecast",
       subtitle = "Pigs in Victoria") +
   theme_ipsum_es()

#b.

resid_pigs<-augment(model_pigs)

#standard deviation
sd <-sd(resid_pigs$.resid)
knitr::kable(paste("Standard Deviation = ",sd),"html")

# Confidence Interval

upper <- fcst_pigs$.mean[1] + 1.96*sd
mean <- fcst_pigs$.mean[1]
lower <- fcst_pigs$.mean[1] - 1.96*sd

knitr::kable(rbind(upper, mean, lower), "html")

# knitr::kable(paste("Upper Limit equals = ",upper),"html")
# knitr::kable(paste("Point Forecast = ", mean),"html")
# knitr::kable(paste("Lower Limit equals = ",lower),"html")

#hilo function - 

int <-hilo(fcst_pigs)

head(int)

int_hilo <- int$`95%`[1]

knitr::kable(paste("Limits equals = ", int_hilo),"html")

#borrowed from Jacob Stemmerich
tail(pigs)

tail(vic_pigs)

#used "pigs" from fpp2 as ts not aus_livestock as tibble

fc = ses(pigs, h=4)
fc_ses <-forecast(fc, h=4)

summary(fc)

summary(fc_ses) %>%
  kbl()%>%
  kable_material(c("striped","hover", full_width = T))

ses_manual = function(y, alpha, level){
  x = level
  for (i in 1:length(y)) {
    x = alpha*y[i] + (1 - alpha)*x
  }
  cat("Forecast of next observation by SES function: ",
      as.character(x),
      sep = "\n")
}

# Find Alpha and Intial Level of fc
a = fc$model$par[1]
l = fc$model$par[2]

# run ses_manual()

knitr::kable(paste(
  "Forecast of next observation is: ", as.character(fc$mean[1])
  ), "html")

#ses_manual(y = pigs, alpha = a, level = l)

fc %>%
  plot(pigs,
       xlab = "Date"
       , ylab = "Slaughter"
       , main = "Pig Slaughter Fcts"
       , col = brewer.pal(4, "Set2"))

################## DID NOT WORK

# ets <-ets(pigs, model = "ANN")
# fc_ets <-forecast(ets, h = 4)
# summary(fc_ets)
# 
# fc_ets %>%
#   plot(vic_pigs$Count)
# 
# ets_manual = function(y1, alpha1, level1){
#   x1 = level1
#   for (i in 1:length(y1)) {
#     x1 = alpha1*y1[i] + (1 - alpha1)*x1
#   }
#   print(x1)
# }
# 
# # Find Alpha and Intial Level of fc
# a1 = ets$model$par[1]
# l1 = ets$model$par[2]
# 
# # run ses_manual()
# ets_manual(y1 = pigs, alpha1 = a1, level1 = l1)

##Borrowed from Miguel Angel Mendez with a few changes.

SES <- function(pars = c(alpha, l0), y){

  ## understood that this is the difference
  error <- 0
  SSE <- 0
  alpha <- pars[1]
  l0 <- pars[2]
  y_hat <- l0
  
  for(index in 1:length(y)){
    error <- y[index] - y_hat
    SSE <- SSE + error^2
    
    y_hat <- alpha*y[index] + (1 - alpha)*y_hat 
  }
  
  return(SSE)
}

opt_SES_pigs <- optim(par = c(0.5, pigs[1]), y = pigs, fn = SES)

knitr::kable(paste(
  "Optimal parameters for alpha and next focst obs using opt function: ",
  "\n",
  as.character(opt_SES_pigs$par[1]),
  ", ",
  as.character(opt_SES_pigs$par[2]),
  sep = ""
  ), "html")

knitr::kable(paste(
  "Optimal parameters for alpha and next focst obs manual ses function: ",
  "\n",
  as.character(a),
  ", ",
  as.character(l),
  sep = ""
  ), "html")

#a.
mex_exports <- global_economy %>%
  filter(Country == "Mexico")

mex_exports %>%
  autoplot(Exports) +
  labs(y = "Count", title = "Mexico exports"
       , x = "Date") +
  theme_ipsum_es() + 
  geom_smooth(formula = y ~ x, method = "loess")

#b.

ets_mex <- mex_exports %>%
  model(ETS(Exports ~ error("A") + trend("N") + season("N")))

fc_mex <- ets_mex %>%
  forecast(h = 10)

fc_mex %>%
  autoplot(mex_exports) + 
  labs(y = "Count", title = "Mexico exports with Forecast"
       , x = "Date") +
  theme_ipsum_es()

knitr::kable(
ets_mex %>%
  accuracy(), "html")

#c.

comp_mex <- mex_exports %>%
  model(
    ANN= ETS(Exports ~error("A")+trend("N")+season("N"))
    ,AAN= ETS(Exports ~error("A")+trend("A")+season("N"))
    #,ANA= ETS(Exports ~error("A")+trend("N")+season("A"))
    )

accuracy(comp_mex) %>%
  kbl()%>%
  kable_material(c("striped","hover", full_width = T))

#d. - #e.

comp_mex %>%
  forecast(h = 10) %>%
  autoplot(mex_exports, level=NULL) +
  labs(y = "Count", title = "Mexico exports with Forecast"
       , x = "Date") +
  theme_ipsum_es()

#f

knitr::kable(comp_mex %>%
  forecast(h=1) %>%
  mutate(interval = hilo(Exports, 95)) %>%
           pull(interval), "html")

comp2 <-comp_mex %>%
  dplyr::select(Country,AAN)%>%
  accuracy()%>%
  transmute(Country,s=RMSE)

knitr::kable(comp_mex %>%
  dplyr::select(Country,ANN)%>%
  forecast(h=1)%>%
  left_join(comp2,by="Country")%>%
  mutate(low= .mean -1.96*comp2$s,
         high= .mean +1.96*comp2$s)%>%
  dplyr::select(Country,Exports,low,high), "html")

#intuition
cna_gdp <- global_economy %>%
  filter(Country=="China") %>%
  mutate(gdp_bill = GDP/1000000000)

cna_gdp %>%
  autoplot(gdp_bill) + 
  labs(y = "GDP in billions", title = "China GDP in USD$"
       , x = "Date") +
  theme_ipsum_es() +
  geom_smooth(method = "loess", formula = y ~ x, color="red", size = 0.5) +
  scale_y_continuous(labels = scales::dollar)

  #+ geom_smooth(formula = y ~ x, method = "loess")

################## modeling

#bring lambda from Box Cox transformation first

lamb_cna <- cna_gdp %>%
  features(gdp_bill, features = guerrero) %>%
  pull(lambda_guerrero)

cna_bc_plot <- cna_gdp %>%
  autoplot(box_cox(gdp_bill, lamb_cna)) + 
  labs(x = "Date"
       , y = "GDP Transformed from Billions of $"
       , title = "China GDP"
       , subtitle = "lambda = -0.03446284"
       ) +
  geom_smooth(method = "loess", formula = y ~ x, color="red", size = 0.5) + 
  theme_ipsum_es() +
  scale_y_continuous(labels = scales::number)

cna_bc_plot

# holy s... variance was removed considerably, smoothed from an exponential move.

# bunch of models

cna_models <- cna_gdp %>%
  model(
    ets1 = ETS(gdp_bill ~ error("A") + trend("A") + season("N")) # no seasonality
    , ets2 = ETS(gdp_bill ~ error("A") + trend("N") + season("N")) #remove the trend--mmm
    , ets_damp = ETS(gdp_bill ~ error("A")+trend("Ad")+season("N"))
    , ets_bc = ETS(box_cox(gdp_bill, lamb_cna))
    , ets_bc_damp = ETS(box_cox(gdp_bill, lamb_cna) ~ trend("Ad"))
    , ets_log = ETS(log(gdp_bill))
  )

cna_models %>%
  forecast(h = 15) %>%
  autoplot(cna_gdp, level = NULL) + 
  labs(x = "Date"
       , y = "GDP in Billions of $"
       , title = "China GDP"
       #, subtitle = "lambda = -0.03446284"
       ) +
  #geom_smooth(method = "loess", formula = y ~ x, color="red", size = 0.8) + 
  theme_ipsum_es() +
  scale_y_continuous(labels = scales::dollar)

cna_models %>% accuracy() %>%
  kbl()%>%
  kable_material(c("striped","hover", full_width = T))

gas_aus <- aus_production %>%
  filter(Quarter > yearquarter("1980 Q1")) %>%
  select(Quarter, Gas)
  
gas_aus %>% 
  autoplot(Gas) +
  labs(y = "Petajouls"
        , x = "Date"
        , title = "Gas Production in Australia") + 
  theme_ipsum_es() + 
  geom_smooth(formula = y ~ x, method = "loess", size = 0.5)

# Model

models_gas <- 
  gas_aus %>%
  model(
    gas_m = ETS(Gas ~ error("M") + trend("A") + season("M"))
    , gas_damp = ETS(Gas ~ error("A") + trend("Ad") + season("M"))
    , gas_damp_sm = ETS(Gas ~ error("A") + trend("Ad",phi=0.90) + season("M"))
    )
 
gas_fc <-
  models_gas %>%
  forecast(h = 10) 
  
gas_fc %>% 
  autoplot(gas_aus, level = NULL) +
  labs(y = "Petajouls", title = "Gas Production in Australia"
       , x = "Date") +
  facet_grid(.model~.)

lambda_turnover <- 
  aus_turnover %>%
  features(Turnover, features = guerrero) %>%
  pull(lambda_guerrero)

aus_turnover_trans <- aus_turnover_train %>%
  mutate(bc = box_cox(Turnover,lambda_turnover))

stl_ets_holt <- aus_turnover_trans %>% 
  model("STL(BoxCox)" = STL(bc ~ season(window="periodic"),robust=T)
        ,"ETS (BoxCox)" = ETS(bc)
        )

#different response variable- separate

holts <- aus_turnover_trans %>%
    model("Holt Winters Damp" =ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
      
rbind(accuracy(stl_ets_holt),accuracy(holts)) %>%
  kbl()%>%
  kable_material(c("striped","hover", full_width = T))