Correlaciones.Rmd

---
title: "R Notebook"
output:
  html_document:
    df_print: paged
fig_width: 12
fig_height: 4
---

# Leemos los df
```{r}
library(ggplot2)
# install.packages('pacman')
pacman::p_load(forecast)
dev.new(width=12, height=5)
knitr::opts_chunk$set(fig.width=12, fig.height=5) 

# install.packages('tidyverse')
df_train = read.csv('data/train.csv')
df_train$Date = as.Date(df_train$Date)
df_test = read.csv('data/test.csv')
df_store = read.csv('data/store.csv')
head(df_store)
head(df_train)
head(df_test)
```

# Store con minima distancia
```{r}
min(df_store[!is.na(df_store$CompetitionDistance),'CompetitionDistance'])
```

```{r}
df_store[df_store$CompetitionDistance==20,]
```


# Obtenemos las stores de test
```{r}
stores_test = unique(df_test$Store)
```


#Filtramos del DF las stores de test
```{r}
df_train_filtered = df_train[df_train$Store %in% stores_test, ]
head(df_train_filtered)
```


# Filtramos la tienda 1
```{r}
df_train_tienda_1 = df_train_filtered[df_train_filtered$Store==1,]
head(df_train_tienda_1)
```

# Filtramos la tienda 3
```{r}
df_train_tienda_3 = df_train_filtered[df_train_filtered$Store==3,]
head(df_train_tienda_3)
```

# Chequeo de la cantidad de días con venta en el día domingo
```{r}
nrow(df_train_filtered[(df_train_filtered$DayOfWeek==7),])
nrow(df_train_filtered[(df_train_filtered$DayOfWeek==7)&(df_train_filtered$Sales==0),])
```

# Correlación cruzada con convolve de la ventas tienda 1, y clientes
```{r fig.width=6, fig.height=4}
df_train_tienda_1$ventas_clientes = convolve(df_train_tienda_1$Sales, df_train_tienda_1$Customers, conj = TRUE)

ggplot(df_train_tienda_1, aes(x = Date, y = ventas_clientes)) +
    geom_line(color="green", size = 1.2) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
    ggtitle("Cantidad de ventas a lo largo de los últimos años")
```


# Correlación cruzada Clientes Ventas
```{r}
ggCcf(df_train_tienda_1$Customers, df_train_tienda_1$Sales,lag.max = 300)
```
# Autocorrelación de la serie de ventas
```{r}
ggAcf(df_train_tienda_1$Sales, lag.max = length(df_train_tienda_1$Sales))
```
# Correlación cruzada tienda 1 y 3
```{r}

ggplot(df_train_tienda_1, aes(x = Date, y = Sales)) +
    geom_line(color="red", size = 1.2) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
    ggtitle("Ventas Tienda 1") +
    xlab("Fecha") +
    ylab("Ventas")
```


```{r}
ggCcf(df_train_tienda_1$Sales,df_train_tienda_3$Sales)
```


```{r}
library(tidyverse)
library(ggplot2)
library(tsibble)
library(feasts)

y <- tsibble(
  Date = df_train_tienda_1$Date,
  Sales = df_train_tienda_1$Sales,
  index = Date
)
autoplot(y,Sales) +
  labs(title = "Ventas diarias Tienda 1",
       y = "Ventas", x='Fecha')

```

#Filtrando cuando se encuentra cerrado
```{r}
melsyd_economy <- df_train_filtered %>%
  filter(Open == 1, Store == 1) %>%
  as_tsibble(index = Date)

autoplot(melsyd_economy, Sales) +
  labs(title = "Ventas - Domingos y feriados filtrados")
```

```{r}
melsyd_economy <- df_train_filtered %>%
  filter(Store == 1) %>%
  mutate(Sales = ifelse(Open == 1, Sales, NA),
         Customers = ifelse(Open == 1, Customers, NA),
         Promo = ifelse(Open == 1, Promo, NA),
         StateHoliday = ifelse(Open == 1, StateHoliday, NA),
         SchoolHoliday = ifelse(Open == 1, SchoolHoliday, NA),
         ) %>%
  as_tsibble(index = Date)
```

# Se puede observar cierta estacionalidad a lo largo del año
```{r}
melsyd_economy %>%
  gg_season(Sales, labels = "both") +
  labs(title = "Seasonal plot")
```

# Estacionalidad semanal
```{r}
melsyd_economy %>% gg_season(Sales, period = "week") +
  theme(legend.position = "none") +
  labs(title="Weekly Seasonality")
```

# Boxplots estacionalidad semanal
```{r}
ggplot(melsyd_economy %>%
  filter(Open==1) %>%
    mutate(DayOfWeek = as.factor(DayOfWeek)),  aes(x= DayOfWeek, y=Sales, color = DayOfWeek)) + 
  geom_boxplot()+
  labs(title = "Ventas diarias Tienda 1",
       y = "Ventas", x='Día de la semana')
```
# La mediana y los extramos de la caja se comportan de manera distinta dependiendo el día de la semana

# Estacionalidad anual
```{r}
melsyd_economy %>% gg_season(Sales, period = "year") +
  theme(legend.position = "none") +
  labs(title="Year Seasonality")
```

#Subseries Plot, estacionalidad dentro de la semana
```{r}
melsyd_economy %>%
  gg_subseries(Sales,period = 'week') +
  labs(title = "Ventas por día de la semana en Tienda 1",
       y = "Ventas", x='Día de la semana')
```


# Venta a lo largo del día del mes
```{r}
melsyd_economy %>%
  gg_subseries(Sales,period = 'month') +
  labs(title = "Ventas por día del mes en Tienda 1",
       y = "Ventas", x='Día del mes')
```

# Boxplot Estacionalidad a lo largo del mes
```{r}
ggplot(melsyd_economy %>%
  filter(Open==1) %>%
    mutate(DayOfMonth = as.factor(lubridate::day(Date))),  aes(x= DayOfMonth, y=Sales, color = DayOfMonth)) + 
  geom_boxplot()+
  labs(title = "Ventas por día del mes para la Tienda 1",
       y = "Ventas", x='Día del mes')
```

# Días para que termine el mes
```{r}
ggplot(melsyd_economy %>%
  filter(Open==1) %>%
    mutate(DayToEndOfMonth = as.factor(lubridate::days_in_month(Date) - lubridate::day(Date))),  aes(x= DayToEndOfMonth, y=Sales, color = DayToEndOfMonth)) + 
  geom_boxplot()
```
# Cuando faltan pocos días para que termine el mes la venta es alta.

# Día posterior a día cerrado
```{r}
ggplot(melsyd_economy %>%
    mutate(Abierto_Dia_Anterior = as.factor(lag(Open))) %>%
    filter(Open==1) ,  aes(x= Abierto_Dia_Anterior, y=Sales, color = Abierto_Dia_Anterior)) + 
  geom_boxplot()+
  labs(title = "Ventas diarias Tienda 1",
       y = "Ventas", x='Se encontraba abierto el día anterior?')
```

# Día posterior a día cerrado, excluyendo lunes, dado que todos los domingos se encuentra cerrado
```{r}
ggplot(melsyd_economy %>%
    mutate(Abierto_Dia_Anterior = as.factor(lag(Open))) %>%
    filter(Open==1, DayOfWeek!=1) ,  aes(x= Abierto_Dia_Anterior, y=Sales, color = Abierto_Dia_Anterior)) + 
  geom_boxplot()
```
# Boxplots Efecto de vacaciones de escuela
```{r}
ggplot(melsyd_economy %>%
  filter(Open==1) %>%
    mutate(SchoolHoliday = as.factor(SchoolHoliday)),  aes(x= SchoolHoliday, y=Sales, color = SchoolHoliday)) + 
  geom_boxplot()+
  labs(title = "Ventas diarias Tienda 1",
       y = "Ventas", x='Ese día es vacaciones escolares?')
```
# La mediana parece ser más baja, no obstante hay mucho solapamiento

```{r}
melsyd_economy %>%
  autoplot(Customers) +
  labs(title = "Clientes diarios Tienda 1",
       y = "Clientes", x='Fecha')
```

# Scatter Plot Clientes y Ventas
```{r}
melsyd_economy %>%
  ggplot(aes(x = Customers, y = Sales)) +
  geom_point()
```

# Alta correlación entre clientes y ventas
```{r}
cor(df_train_tienda_1$Customers, df_train_tienda_1$Sales)
cor(df_train_tienda_1$Customers, df_train_tienda_1$Sales, method = 'spearman')
```
# Lag Plots
# Scatter plot contra diferentes lags, a medida que nos alejamos en el tiempo, se encuentran menos relacionadas
```{r}
melsyd_economy %>%
  gg_lag(Sales, geom = "point") +
  labs(x = "lag(Sales, k)")
```


# Scatter Plot para diferentes lags
```{r}
melsyd_economy %>%
  gg_lag(Sales, geom = "point", lags = c(1,2,3,7,14,30,365)   ) +
  labs(x = "lag(Sales, k)")
```

# Autocorrelacion 40 lags
```{r}
melsyd_economy %>%
  ACF(Sales, lag_max = 40, na.action = na.pass) %>%
  autoplot() + labs(title="Sales")
```


# Autocorrelacion 400 lags
```{r}
melsyd_economy %>%
  ACF(Sales, lag_max = 400, na.action = na.pass) %>%
  autoplot() + labs(title="Sales")
```
# Para lags pequeños tenemos una alta correlación, cada 7 días se ven picos, lo que indica que hay cierta estacionalidad
# También vemos un pico grande aproximadamente al año
# Cuando hay tendencia, suele haber una autocorrelación alta y positiva en los lags pequeños, valores cercanos en tiempo, cercanos también en valor
# Cuando la data es estacional, las autocorrelaciones son altas para los lags estacionales (en múltiplos de dichos lags)


# Descomponemos la serie en Tendencia, Estacionalidad Anual y Estacionalidad Semanal
```{r}
melsyd_economy_2 <- df_train_filtered %>%
  filter(Store == 1) %>%
  # mutate(Sales = ifelse(Open == 1, Sales, NA),
  #        Customers = ifelse(Open == 1, Customers, NA),
  #        Promo = ifelse(Open == 1, Promo, NA),
  #        StateHoliday = ifelse(Open == 1, StateHoliday, NA),
  #        SchoolHoliday = ifelse(Open == 1, SchoolHoliday, NA),
  #        ) %>%
  as_tsibble(index = Date)

dcmp <- melsyd_economy_2 %>%
  model(stl = STL(Sales))
components(dcmp)
melsyd_economy_2 %>% autoplot(Sales)
components(dcmp) %>% autoplot()

```


```{r}
components(dcmp) %>%
  as_tsibble() %>%
  autoplot(Sales, colour="gray") +
  geom_line(aes(y=trend), colour = "#D55E00") +
  labs(

    title = "Sales"
  )
```

# Datos con estacionalidad ajustada
```{r}
components(dcmp) %>%
  as_tsibble() %>%
  autoplot(Sales, colour = "gray") +
  geom_line(aes(y=season_adjust), colour = "#0072B2") +
  labs(title = "Ventas diarias Tienda 1 <estacionalidad ajustada>",
       y = "Ventas", x='Fecha')
```


# Moving Average 7
```{r}
# Considerando la serie completa
melsyd_economy_3 <- melsyd_economy_2 %>%
  mutate(
    `7-MA` = slider::slide_dbl(Sales, mean,
                .before = 3, .after = 3, .complete = TRUE)
  )

melsyd_economy_3 %>%
  autoplot(Sales) +
  geom_line(aes(y = `7-MA`), colour = "#D55E00") +
  labs(y = "Ventas",
       title = "Ventas diarias Tienda 1") +
  guides(colour = guide_legend(title = "series"))


# Filtrando cuando se encuentra cerrado
melsyd_economy_3 <- melsyd_economy %>%
  mutate(
    `7-MA` = slider::slide_dbl(Sales, mean,
                .before = 3, .after = 3, .complete = TRUE)
  )


melsyd_economy_3 %>%
  autoplot(Sales) +
  geom_line(aes(y = `7-MA`), colour = "#D55E00") +
  labs(y = "Ventas",
       title = "Ventas diarias Tienda 1") +
  guides(colour = guide_legend(title = "series"))

```


# Se puede calcular tendencia a partir de dos filtros moving average, en este caso, uno de 7 para la semana, y despues uno de 52 para el año  (*) habría que chequearlo...
```{r}
melsyd_economy_4 <- melsyd_economy_2 %>%
  mutate(
    `12-MA` = slider::slide_dbl(Sales, mean,
                .before = 3, .after = 3, .complete = TRUE),
    `2x12-MA` = slider::slide_dbl(`12-MA`, mean,
                .before = 26, .after = 26, .complete = TRUE)
  )
melsyd_economy_4 %>%
  autoplot(Sales, colour = "gray") +
  geom_line(aes(y = `2x12-MA`), colour = "#D55E00") +
  labs(y = "Ventas",
       title = "Ventas diarias Tienda 1 <tendencia como suma de dos filtros MA>")
```

# Descomposición clásica
```{r}
melsyd_economy_2 %>%
  model(
    classical_decomposition(Sales, type = "additive")
    # classical_decomposition(Sales~season(365), type = "additive")
  ) %>%
  components() %>%
  autoplot() +
  labs(title = "Descomposición lineal aditiva para las ventas Tienda 1")
```

# Descomposiciones clásicas
```{r}
# install.packages("seasonal")
library(seasonal)
# No funca
# x11_dcmp <- melsyd_economy_2[,c('Date','Sales')] %>%
#   model(x11 = X_13ARIMA_SEATS(Sales ~ x11(na.action = seasonal::na.x13))) %>%
#   components()
# autoplot(x11_dcmp) +
#   labs(title =
#     "Decomposition of total US retail employment using X-11.")
# 
# sum(is.na(melsyd_economy_2$Sales))
```

# STL Decomposition, varias ventajas, lo único, no considera estacionalidad intradiaria o dentro del calendario
```{r}
melsyd_economy_2 %>%
  model(
    STL(Sales ~ trend(window = 360) +
                   season(window = "periodic"),
    robust = TRUE)) %>%
  components() %>%
  autoplot()+
  labs(title = "Descomposición STL para las ventas Tienda 1")
```
# STL Features
```{r}
melsyd_economy_2 %>%
  features(Sales, feat_stl)
```

# BenchMark Models

# Filter Data for train
```{r}
df_train_tienda_1_train <- melsyd_economy_2 %>%
  filter(Date<'2015-01-01')
```


```{r}
# install.packages('fable')
library(fable)
```


# Mean
```{r}
df_train_tienda_1_train %>% model(MEAN(Sales))
```
# Naive last value
```{r}
df_train_tienda_1_train %>% model(NAIVE(Sales))
```

# Naive last cycle (seasonal) Year
```{r}
df_train_tienda_1_train %>% model(SNAIVE(Sales ~ lag("year")))
```

# Naive last cycle (seasonal) Week
```{r}
df_train_tienda_1_train %>% model(SNAIVE(Sales ~ lag("week")))
```

# Drift
```{r}
df_train_tienda_1_train %>% model(RW(Sales ~ drift()))
```

# Different Models including 0
```{r}
# Set training data from 1992 to 2006
# train <- aus_production %>%
#   filter_index("1992 Q1" ~ "2006 Q4")
# Fit the models
beer_fit <- df_train_tienda_1_train %>%
  model(
    Mean = MEAN(Sales),
    `Naïve` = NAIVE(Sales),
    `Seasonal naïve` = SNAIVE(Sales)
  )
# Generate forecasts for 14 quarters
beer_fc <- beer_fit %>% forecast(h = 42)
# Plot forecasts against actual values
beer_fc %>%
  autoplot(df_train_tienda_1_train, level = NULL) +
  autolayer( df_train_tienda_1_train %>%
    filter(Date>="2015-01-01") %>%
      select(Sales),
    colour = "black"
  ) +
  labs(
    title = "Forecasts para las ventas Tienda 1"
  ) +
  guides(colour = guide_legend(title = "Forecast"))
```

# Different Models Zeros as Null
```{r}
# Set training data from 1992 to 2006
# train <- aus_production %>%
#   filter_index("1992 Q1" ~ "2006 Q4")
# Fit the models
beer_fit <- melsyd_economy %>%
  model(
    # Actual_Sales = Sales,
    Mean = MEAN(Sales),
    `Naïve` = NAIVE(Sales),
    `Seasonal naïve (week)` = SNAIVE(Sales ~ lag("week")),
    `Seasonal naïve (year)` = SNAIVE(Sales ~ lag("year")),
    
                              )
# Generate forecasts for 14 quarters
beer_fc <- beer_fit %>% forecast(h = 42)
# Plot forecasts against actual values
beer_fc %>%
  autoplot(melsyd_economy 
            # %>% filter(Date>="2014-10-20")
             , level = NULL) +
  autolayer( melsyd_economy %>%
    filter(Date>="2015-01-01") %>%
      select(Sales),
    colour = "black"
  ) +
  labs(
    title = "Forecasts para las ventas Tienda 1"
  ) +
  guides(colour = guide_legend(title = "Forecast"))
```

# Simple Models
```{r}
# Set training data from 1992 to 2006
# train <- aus_production %>%
#   filter_index("1992 Q1" ~ "2006 Q4")
# Fit the models
beer_fit <- melsyd_economy %>%
  filter(Date<"2015-01-01") %>%
  model(
    # Actual_Sales = Sales,
    Mean = MEAN(Sales),
    `Naïve` = NAIVE(Sales),
    `Seasonal naïve (week)` = SNAIVE(Sales ~ lag("week")),
    # `Seasonal naïve (year)` = SNAIVE(Sales ~ lag("year")),
    
                              )
# Generate forecasts for 14 quarters
beer_fc <- beer_fit %>% forecast(h = 42)
# Plot forecasts against actual values
beer_fc %>%
  autoplot(melsyd_economy 
            %>% filter(Date>="2014-10-20")
             , level = NULL) +
  autolayer( melsyd_economy %>%
    filter(Date>="2015-01-01") %>%
      select(Sales),
    colour = "black"
  ) +
  labs(
    title = "Forecasts para las ventas Tienda 1"
  ) +
  guides(colour = guide_legend(title = "Forecast"))
```

# Ajustando el formato
```{r}
df_train_tienda_1 <- df_train_tienda_1 %>%
  as_tsibble(index = Date) 
```

# Simple Models
```{r}
# Set training data from 1992 to 2006
# train <- aus_production %>%
#   filter_index("1992 Q1" ~ "2006 Q4")
# Fit the models
beer_fit <- df_train_tienda_1 %>%
  filter(Date<"2015-01-01") %>%
  model(
    # Actual_Sales = Sales,
    Mean = MEAN(Sales),
    `Naïve` = NAIVE(Sales),
    `Seasonal naïve (week)` = SNAIVE(Sales ~ lag("week")),
    # `Seasonal naïve (year)` = SNAIVE(Sales ~ lag("year")),
    
                              )
# Generate forecasts for 14 quarters
beer_fc <- beer_fit %>% forecast(h = 42)
# Plot forecasts against actual values
beer_fc %>%
  autoplot(df_train_tienda_1 
            %>% filter(Date>="2014-10-20")
             , level = NULL) +
  autolayer( df_train_tienda_1 %>%
    filter(Date>="2015-01-01") %>%
      select(Sales),
    colour = "black"
  ) +
  labs(
    title = "Forecasts para las ventas Tienda 1"
  ) +
  guides(colour = guide_legend(title = "Forecast"))
```

# Para ver los valoes del fit, residuos y residuos innovados (cambio de escala)
```{r}
augment(beer_fit)
```

<!-- 5.4 Residual diagnostics -->
<!-- A good forecasting method will yield innovation residuals with the following properties: -->

<!-- The innovation residuals are uncorrelated. If there are correlations between innovation residuals, then there is information left in the residuals which should be used in computing forecasts. -->
<!-- The innovation residuals have zero mean. If they have a mean other than zero, then the forecasts are biased. -->

<!-- In addition to these essential properties, it is useful (but not necessary) for the residuals to also have the following two properties. -->

<!-- The innovation residuals have constant variance. This is known as “homoscedasticity”. -->
<!-- The innovation residuals are normally distributed. -->

# Residuos
```{r}
aug <- df_train_tienda_1 %>%
  model(SNAIVE(Sales ~ lag("week"))) %>%
  augment()
autoplot(aug, .innov) +
  labs(title = "Residuos del método naïve")
```

# Histograma Residuos
```{r}
aug %>%
  ggplot(aes(x = .innov)) +
  geom_histogram() +
  labs(title = "Histograma de los residuos")
```

# Autocorrelacion de los residuos
```{r}
aug %>%
  ACF(.innov) %>%
  autoplot() +
  labs(title = "Residuos del método naïve")
```
# Se pueden ver algunas autocorrelaciones altas.

# Tres gráficos en 1, Residuos a lo largo del tiempo, autocorrelación e histograma
```{r}
df_train_tienda_1 %>%
  model(SNAIVE(Sales ~ lag("week"))) %>%
  gg_tsresiduals()
```

# Intervalos de predicción
```{r}
df_train_tienda_1 %>%
  model(SNAIVE(Sales ~ lag("week"))) %>%
  forecast(h = 42) %>%
  autoplot(df_train_tienda_1) +
  labs(title="Forecast de Ventas <intervalos de predicción SNAIVE>", y="$US" )
```

# Intervalos de predicción
```{r}
melsyd_economy %>%
  model(SNAIVE(Sales ~ lag("week"))) %>%
  forecast(h = 42) %>%
  autoplot(melsyd_economy) +
  labs(title="Forecast de ventas <Intervalos de predicción SNAIVE>", y="Ventas" )
```

# Intervalos de predicción
```{r}
melsyd_economy %>%
  model(NAIVE(Sales)) %>%
  forecast(h = 42) %>%
  autoplot(melsyd_economy) +
  labs(title="Forecast de ventas <Intervalos de predicción SNAIVE>", y="Ventas" )
```

# Forecast for the last 42 days (6 weeks)
```{r}
recent_production <- melsyd_economy %>%
  slice(n()-42:0)
beer_train <- melsyd_economy %>%
  slice(1:(n()-43))

beer_fit <- beer_train %>%
  model(
    Mean = MEAN(Sales),
    `Naïve` = NAIVE(Sales),
    `Seasonal naïve` = SNAIVE(Sales),
    Drift = RW(Sales ~ drift())
  )

beer_fc <- beer_fit %>%
  forecast(h = 42)

beer_fc %>%
  autoplot(
    melsyd_economy %>% slice(n()-134:0),
    level = NULL
  ) +
  labs(title="Forecast de ventas", y="Ventas" ) +
  guides(colour = guide_legend(title = "Forecast"))
```

# Accuracy metrics
```{r}
accuracy(beer_fc, recent_production)
```


# Forecast for the last 42 days (6 weeks)
```{r}
recent_production <- melsyd_economy_2 %>%
  slice(n()-42:0)
beer_train <- melsyd_economy_2 %>%
  slice(1:(n()-43))

beer_fit <- beer_train %>%
  model(
    Mean = MEAN(Sales),
    `Naïve` = NAIVE(Sales),
    `Seasonal naïve` = SNAIVE(Sales),
    Drift = RW(Sales ~ drift())
  )

beer_fc <- beer_fit %>%
  forecast(h = 42)

beer_fc %>%
  autoplot(
    melsyd_economy_2 %>% slice(n()-134:0),
    level = NULL
  ) +
  labs(title="Forecast de ventas", y="Ventas" ) +
  guides(colour = guide_legend(title = "Forecast"))
```

# Accuracy metrics evaluación del modelo
```{r}
accuracy(beer_fc, recent_production)
```

# Walk forward validation (CV cross validation)

# Create Training Test Partitions
```{r}
# Time series cross-validation accuracy
google_2015_tr <- melsyd_economy_2 %>%
  stretch_tsibble(.init = 731, .step = 42) %>% # Toma al menos los primeros dos años y va tomando de a 6 semanas para testear
  relocate(Date, .id)
google_2015_tr
```

# Cross validation
```{r}
rbind(
# TSCV accuracy
google_2015_tr %>%
  model(RW(Sales ~ drift())) %>%
  forecast(h = 1) %>%
  accuracy(melsyd_economy_2)
,
# Training set accuracy
melsyd_economy_2 %>%
  model(RW(Sales ~ drift())) %>%
  accuracy()
)
```

# Cross validation
```{r}

# TSCV accuracy
rbind(
google_2015_tr %>%
  model(SNAIVE(Sales)) %>%
  forecast(h = 42) %>%
  accuracy(melsyd_economy_2)
,
# Training set accuracy
melsyd_economy_2 %>%
  model(SNAIVE(Sales)) %>%
  accuracy()
)
```

# Regresion Lineal
```{r}
fit_beer <- recent_production %>%
  model(TSLM(Sales ~ trend() + season()))
report(fit_beer)
```

# Regresion Lineal
```{r}
melsyd_economy_rl = melsyd_economy_2 %>%
  filter(Date<"2015-06-19")
melsyd_economy_rl$DayOfWeek = as.factor(melsyd_economy_rl$DayOfWeek)

fit_beer <- melsyd_economy_rl %>%
  model(TSLM(Sales ~ Open * (DayOfWeek + Date) - Open - DayOfWeek - Date))
  # model(TSLM(Sales ~ Date))
report(fit_beer)
```

```{r}
recent_production <- melsyd_economy_2 %>%
  slice(n()-42:0)
recent_production$DayOfWeek = as.factor(recent_production$DayOfWeek)

fc_beer <- forecast(fit_beer, new_data = recent_production)

fc_beer %>%
  autoplot(recent_production) +
  labs(
    title = "Forecast de ventas usando regresión"
  )
```

# Métricas del modelo
```{r}
accuracy(fc_beer, recent_production)

```


# Exponential Smoothing
```{r}
# # install.packages('forecast')
# library(Rcpp)
# library(forecast)
# # Estimate parameters
fit <- melsyd_economy_2 %>%
  filter(Date<"2015-06-19") %>%
  model(fable::ETS(Sales ~ error("A") + trend("N") + season("N")))
fc <- fit %>%
  forecast(h = 42)
```

# Plot
```{r}
fc %>%
  autoplot(melsyd_economy_2 %>%
             slice(n()-126:0)) +
  labs(title="Exponential Smoothing") +
  guides(colour = "none")
```

# Métricas del modelo
```{r}
accuracy(fc, recent_production )

```

# Modelos con tendencia y estacionalidad aditiva y multiplicativa
```{r}
aus_holidays <- melsyd_economy_2 %>%
  filter(Date<"2015-06-19") %>%
  summarise(Sales = sum(Sales)/1e3)

fit <- aus_holidays %>%
  model(
    additive = ETS(Sales ~ error("A") + trend("A") +
                                                season("A")),
    multiplicative = ETS(Sales ~ error("M") + trend("A") +
                                                season("M"))
  )
fc <- fit %>% forecast(h = 42)
fc %>%
  autoplot(aus_holidays %>%
           slice(n()-126:0), level = NULL) +
  labs(title="Sales <Estacionalidad Holts Winters> - Additive and Multiplicative") +
  guides(colour = guide_legend(title = "Forecast"))
```

# Métricas del modelo
```{r}
accuracy(fc, recent_production %>%
           mutate(Sales = Sales)  )
```

# Efecto Multiplicativo
```{r}
sth_cross_ped <- melsyd_economy_2 %>%
  summarise(Sales = sum(Sales)/1e3)

sth_cross_ped %>%
  filter(Date<"2015-06-19") %>%
  model(
    hw = ETS(Sales ~ error("M") + trend("Ad") + season("M"))
  ) %>%
  forecast(h = 42) %>%
  autoplot(sth_cross_ped %>% slice(n()-126:0)) +
  labs(title = "Estacionalidad Holts - Winters Multiplicativo")
```

# Descomposición del modelo Holts-Winters en Nivel, Tendencia, Estacionalidad
```{r}
dcmp <- melsyd_economy_2 %>%
  summarise(Sales = sum(Sales)/1e3) %>%
    model(
    hw = ETS(Sales ~ error("M") + trend("Ad") + season("M"))
  ) 
components(dcmp)
melsyd_economy_2 %>% autoplot(Sales)
components(dcmp) %>% autoplot()
```


# Descomposición del modelo Tendancia y Estacionalidad aditivas en Nivel, Tendencia, Estacionalidad
```{r}
dcmp <- melsyd_economy_2 %>%
  summarise(Sales = sum(Sales)/1e3) %>%
    model(
    hw = ETS(Sales ~ error("A") + trend("Ad") + season("A"))
  ) 
components(dcmp)
melsyd_economy_2 %>% autoplot(Sales)
components(dcmp) %>% autoplot()
```

# ARIMA
```{r}
melsyd_economy_2 %>%
  mutate(diff_sales = difference(Sales)) %>%
  features(diff_sales, ljung_box, lag = 10)
```

```{r}
melsyd_economy_2 %>%
  summarise(Sales = sum(Sales)/1e6) %>%
  transmute(
    `Sales ($million)` = Sales,
    `Log sales` = log(Sales),
    `Annual change in log sales` = difference(log(Sales), 365),
    `Doubly differenced log sales` =
                     difference(difference(log(Sales), 365), 1)
  ) %>%
  pivot_longer(-Date, names_to="Type", values_to="Sales") %>%
  mutate(
    Type = factor(Type, levels = c(
      "Sales ($million)",
      "Log sales",
      "Annual change in log sales",
      "Doubly differenced log sales"))
  ) %>%
  ggplot(aes(x = Date, y = Sales)) +
  geom_line() +
  facet_grid(vars(Type), scales = "free_y") +
  labs(title = "Corticosteroid drug sales", y = NULL)
```

# Test de estacionariedad
```{r}
melsyd_economy_2 %>%
  features(Sales, unitroot_kpss)
```

```{r}
melsyd_economy_2 %>%
mutate(diff_close = difference(Sales, 365)) %>%
features(diff_close, unitroot_kpss)
```

```{r}
melsyd_economy_2 %>%
mutate(diff_close = difference(Sales, 7)) %>%
features(diff_close, unitroot_kpss)
```

```{r}
melsyd_economy_2 %>%
  features(Sales, unitroot_ndiffs)
```

# ARIMA con estacionalidad

# Estacionalidad Anual
```{r}
melsyd_economy_2 %>%
  gg_tsdisplay(difference(Sales, 365),
               plot_type='partial', lag=36) +
  labs(title="Seasonally differenced", y="")
```

# Estacionalidad semanal
```{r}
melsyd_economy_2 %>%
  gg_tsdisplay(difference(Sales, 7),
               plot_type='partial', lag=36) +
  labs(title="Seasonally differenced", y="")
```

# Restando diferencia de orden 1 adicionalmente
```{r}
melsyd_economy_2 %>%
  gg_tsdisplay(difference(Sales, 28) %>% difference(),
               plot_type='partial', lag=36) +
  labs(title = "Double differenced", y="")
```

# Removiendo Nulos

# Estacionalidad Anual
```{r}
melsyd_economy %>%
  gg_tsdisplay(difference(Sales, 365),
               plot_type='partial', lag=36) +
  labs(title="Seasonally differenced", y="")
```

# Estacionalidad semanal
```{r}
melsyd_economy %>%
  gg_tsdisplay(difference(Sales, 7),
               plot_type='partial', lag=36) +
  labs(title="Seasonally differenced", y="")
```

# Restando diferencia de orden 1 adicionalmente
```{r}
melsyd_economy_2 %>%
  gg_tsdisplay(difference(Sales, 28) %>% difference(),
               plot_type='partial', lag=36) +
  labs(title = "Double differenced", y="")
```


```{r}
# install.packages('fable.prophet')
library(fable.prophet)

train <- melsyd_economy_2 %>%
  filter(Date<"2015-06-19")
fit <- train %>%
  model(
    arima = ARIMA(Sales),
    ets = ETS(Sales),
    prophet = prophet(Sales ~ season(type = "multiplicative"))
  )
```

```{r}
fc <- fit %>% forecast(h = 42)
fc %>% autoplot(melsyd_economy_2 %>%
                  slice(n()-134:0)
                  )
```

# Neural Nets (Redes Neuronales)
```{r}
melsyd_economy_2 %>%
  filter(Date<"2015-06-19") %>%
  model(NNETAR(Sales)) %>%
  forecast(h = 42, times = 100) %>% # El parámetro times permite correr distintas simulaciones, cambian los intervalos de confianza
  autoplot(melsyd_economy_2 %>%
                  slice(n()-134:0)) +
  labs(x = "Year", y = "Counts",
       title = "Yearly sunspots")
```

# Seleccion del mejor modelo

# Walk forward validation (CV cross validation)

# Create Training - Val Partitions para WFV
```{r}
# Time series cross-validation accuracy
df_cv <- melsyd_economy_2 %>%
  slice(0:(n()-85)) %>%
  stretch_tsibble(.init = 731, .step = 42) %>% # Toma al menos los primeros dos años y va tomando de a 6 semanas para testear
  relocate(Date, .id)


print(paste("Cantidad de folds: ", max(df_cv$.id)))
df_cv
```

# Cross validation
```{r}
rbind(
# TSCV accuracy
df_cv %>%
    model(
    Mean = MEAN(Sales),
    `Naïve` = NAIVE(Sales),
    `Seasonal naïve` = SNAIVE(Sales),
    Drift = RW(Sales ~ drift()),
    ETS = ETS(Sales),
    # Prophet = prophet(Sales ~ season(type = "multiplicative")),
    Prophet = prophet(Sales),
    # Redes_Neuronales = model(NNETAR(Sales))
  )
 %>%
  forecast(h =42) %>%
  accuracy(melsyd_economy_2) # Filtro sólo los días que estuvo abierto 
# ,
# # Training set accuracy
# melsyd_economy_2 %>%
#   model(RW(Sales ~ drift())) %>%
#   accuracy()
)
```

# Moving Average 7 días antes

# Create Training - Val Partitions para WFV
```{r}
# Time series cross-validation accuracy
df_cv <- melsyd_economy_2 %>%
    mutate(
    Sales = slider::slide_dbl(Sales, mean,
                .before = 6, .after = 0, .complete = FALSE)) %>%
  slice(0:(n()-85)) %>%
  stretch_tsibble(.init = 731, .step = 42) %>% # Toma al menos los primeros dos años y va tomando de a 6 semanas para testear
  relocate(Date, .id)


print(paste("Cantidad de folds: ", max(df_cv$.id)))
df_cv
```

# Cross validation
```{r}
rbind(
# TSCV accuracy
df_cv %>%
    model(
    Mean = MEAN(Sales),
    `Naïve` = NAIVE(Sales),
    `Seasonal naïve` = SNAIVE(Sales),
    Drift = RW(Sales ~ drift()),
    ETS = ETS(Sales),
    # Prophet = prophet(Sales ~ season(type = "multiplicative")),
    Prophet = prophet(Sales),
    # Redes_Neuronales = model(NNETAR(Sales))
  )
 %>%
  forecast(h =42) %>%
  accuracy(melsyd_economy_2 ) # Filtro sólo los días que estuvo abierto 
# ,
# # Training set accuracy
# melsyd_economy_2 %>%
#   model(RW(Sales ~ drift())) %>%
#   accuracy()
)
```

# Moving Average 28 días antes

# Create Training - Val Partitions para WFV
```{r}
# Time series cross-validation accuracy
df_cv <- melsyd_economy_2 %>%
    mutate(
    Sales = slider::slide_dbl(Sales, mean,
                .before = 27, .after = 0, .complete = FALSE)) %>%
  slice(0:(n()-85)) %>%
  stretch_tsibble(.init = 731, .step = 11) %>% # Toma al menos los primeros dos años y va tomando de a 6 semanas para testear
  relocate(Date, .id)


print(paste("Cantidad de folds: ", max(df_cv$.id)))
df_cv
```

# Cross validation
```{r}
rbind(
# TSCV accuracy
df_cv %>%
    model(
    Mean = MEAN(Sales),
    `Naïve` = NAIVE(Sales),
    `Seasonal naïve` = SNAIVE(Sales),
    Drift = RW(Sales ~ drift()),
    ETS = ETS(Sales),
    # Prophet = prophet(Sales ~ season(type = "multiplicative")),
    Prophet = prophet(Sales),
    # Redes_Neuronales = model(NNETAR(Sales))
  )
 %>%
  forecast(h = 42) %>%
  accuracy(melsyd_economy_2 ) # Filtro sólo los días que estuvo abierto 
# ,
# # Training set accuracy
# melsyd_economy_2 %>%
#   model(RW(Sales ~ drift())) %>%
#   accuracy()
)
```

# Resultados en TEST (los últimos 42 días)


```{r}
melsyd_economy_2 %>%
  slice(0:(n()-42)) %>%
    model(
      Prophet = prophet(Sales)
    )%>%
  forecast(h = 42) %>%
  accuracy(melsyd_economy_2)
```

# Gráfico en Test
```{r}
melsyd_economy_2 %>%
  slice(0:(n()-42)) %>%
  # filter(Open==1) %>%
    model(
      Prophet = prophet(Sales)
    )%>%
  forecast(h = 42) %>%
  autoplot(melsyd_economy_2 %>%
             # filter(Open==1) %>%
                  slice(n()-134:0)) +
  labs(title = "Ventas últimas 18 semanas y predicciones con Prophet")
```

# Descomposición Prophet
```{r}
fit <- melsyd_economy_2 %>%
  slice(0:(n()-42)) %>%
  model(
    prophet(Sales) ) 
fit %>%
  components() %>%
  autoplot()
```

# Análisis de residuos
```{r}
fit %>% gg_tsresiduals()
```