---
title: "Forecasating the sales data"
output: 
  flexdashboard::flex_dashboard:
             social: menu
             source: embed
             vertical_layout: fill
             orientation: columns
             theme: cerulean
             
---


```{r,echo=FALSE,include=FALSE}
rm(list = ls())

#install.packages("tidyverse")
#install.packages("datatable")

library(dplyr)
library(tidyverse)
library(lubridate)
library(forecast) #boxcox
library(tseries)
library(vtable)
library(flexdashboard) #Presentation quality
library(data.table)



```


INTRODUCTION
=======================================================================

## ABOUT DATA {.tabset}





We are trying to forecast a time series data of a retail company.
The data is about sales of the thousands of products of different families sold at stores.The training data includes dates, store and product information, whether that item was being promoted, as well as the sales numbers.
The data is publically available on Kaggle website.
The dataset have sale details for various products on day to day basis. The data has various aspects which can be explored going ahead. For now, we will be using single time series of total daily sales across all the products and try to forecast the sales with different approaches.
I find this tak challenging as it involves various factors such as seasonality,holiday events along with forecasting sales.

### Data summary

#### File Descriptions and Data Field Information

sales_train.csv


The training data, comprising time series of features store_nbr, family, and on promotion as well as the target sales.
store_nbr identifies the store at which the products are sold.
family identifies the type of product sold.
on promotion gives the total number of items in a product family that were being promoted at a store at a given date.


```{r,echo=FALSE,include=FALSE}

  sales <- readr::read_csv('D:\\VAIBHAV\\HOMEWORK\\Time Series\\ASSIGNMENT\\sales_train.csv')
  


```


#### Structure of Data


```{r,echo=FALSE}

#str(sales)

sales$date<- dmy(sales$date)

str(sales)

```
```{r,echo=FALSE}

summary(sales)

```



### BASIC DATA QUALITY CHECKS



#### Null value check


```{r,echo=FALSE}
attach(sales)
print("Unique block numbers are the months in the data")
unique(date_block_num)

print("Number of NA values in columns are")
lapply(sales,function(x) { length(which(is.na(sales)))})

```

#### Summary

```{r,echo=FALSE}

summary(sales)
```




We could see in summry that there is some inconsistancy in price and item count per day column. We have some negative values which does not make sense. Lets see in detail.
```{r,echo=FALSE}

print('Number of negative values in item_cnt_day')
sum(item_cnt_day<0)

print("Number of negative values in item_price")
sum(item_price<0)
```

#### Boxplot for both columns
```{r,echo=FALSE}

boxplot(item_price,main = " Item Price")

boxplot(item_cnt_day, main = "Item count per day")

```

There is value in item price which is too far from other values, as it is only one value there is high chance of error, so we replace that by median and plot new boxplot.

```{r,echo=FALSE}


item_price[item_price < 0] <- 0

item_cnt_day[item_cnt_day < 0] <- 0

item_price[item_price > 307970] <- median(item_price) 
```



#### Updated box Plot
```{r,echo=FALSE}

boxplot(item_price,main="Item Price")

boxplot(item_cnt_day,main="Item count per day")

```



#### Summary statistics 
```{r,echo=FALSE}
vtable(sales,lush = TRUE)
```


Creating our time series to track total daily sales. Keeping only two columns hereafter.
```{r,echo=FALSE}

sales <-   sales %>% group_by(date) %>% 
           mutate(daily_Sales= sum(item_price*item_cnt_day)/100000) %>% 
           arrange(date) %>% 
           ungroup()


sales <- subset(sales,select = c(date,daily_Sales))

sales <- sales %>% group_by(date) %>% filter(!duplicated(date))

```

```{r,echo=FALSE}

summary(sales)
```


#### Daily total sales plot

```{r,echo=FALSE}

sales %>% 
ggplot()+
  geom_line(aes(date,daily_Sales))+
    theme_bw()+
    xlab("Date")+
    ylab("TOTAL DAILY SALES")

```

Looking at plot, the series looks quite stationary, but have seasonal factors involved. The trend and variance seems to be same overall series.




#### Cummulative sales per day.
```{r,echo=FALSE}
sales %>% ggplot(aes(colour="Blue"))+
             
             geom_line(aes(sales$date,cumsum(sales$daily_Sales)),size=1.5,show.legend = FALSE)+
             theme_bw()+
             ggtitle("TOTAL CUMMULATIVE SALES BY DATE")+
             ylab("Cummulative Sales")+
             xlab("DATE")

#?GeomLine
             
```














APPROACHING TO MODELLING
=======================================================================



## Understanding the series {.tabset}

### UNDERSTANDING THE SERIES

```{r,echo=FALSE,include=FALSE}

sales %>% ggplot()+
             geom_line(aes(sales$date,sales$daily_Sales))+
              theme_bw()+
              #geom_smooth(aes(sales$date,sales$daily_Sales),method = "lm")+
              ggtitle("TOATAL DAILY SALES")+
                xlab("Date")+
                ylab("Total daily sales")
              

```

The series appeared to be bit stationary before in terms of both variance and mean. Plotting again with trend line.The trend is bit constant here.



#### Lets see if we have any auto-correlation in the series to move ahead
#### ACF & PACF GRAPHS
```{r,echo=FALSE,message=FALSE,warning=FALSE}

attach(sales)

par(mfrow=c(1,2))
acf(sales$daily_Sales,lag.max = 40)
Pacf(sales$daily_Sales, lag.max = 40)

```


From above plots we could see that there is strong autocorrelation in the data



#### Checking wheather the given series is stationary or not

```{r,echo=FALSE,warning=FALSE}

adf.test(sales$daily_Sales)
kpss.test(sales$daily_Sales)

```

The ADF test has Null hypothesis that given series is not stationary and KPSS test has Null Hypothesis as the given series is stationary. Thus both test suggest that given series is already stationary.


We will try to check how much difference it makes if we transformed the given series.

#### Lets check the histogram of our sales
```{r,echo=FALSE}

hist(sales$daily_Sales,main="Daily sales",xlab = "DAILY SALES", ylab = "Frequency")



```

The distribution looks right skewed.



#### We will log transform the series and plot the histogram.

```{r,echo=FALSE}

sales_1 <- sales

train_log = sales_1 %>%
  mutate(sale_log = log1p(daily_Sales))

train_log %>%
  ggplot()+
      geom_line(aes(date,sale_log))+
      theme_bw()+
      ggtitle("SALES over Time (Log)")+
      ylab("DAILY SALE (Log)")+
      xlab("Date")


```

Now the distribution of the sales is close to normal. But Remember, we dont look for the ditribution of Y variable to be normal, that is not our concern. So we just leave it here.






#### Lets try box-cox transformation and check
```{r,echo=FALSE}

train_box_cox <- sales_1 %>% 
                 mutate(sale_box_cox = forecast::BoxCox(daily_Sales,lambda="auto"))

train_box_cox %>% ggplot()+
                  geom_line(aes(date,sale_box_cox))+
                  theme_bw()+
                  ggtitle("SALES PER DAY-BOX-COX TRANSFORMED")+
                  ylab("TOTAL SALES")



```

There seem to be no difference in original and transformed series here, we have already checked that our series is stationary. Box-Cox performed well here by not doing anything.



#### Comparing histogram of all three
```{r,echo=FALSE}

hist(sales$daily_Sales,main="Daily sales",xlab = "DAILY SALES", ylab = "Frequency")

hist(train_log$sale_log,main="Daily sales Log transformed",xlab = "DAILY SALES", ylab = "Frequency")


hist(train_box_cox$sale_box_cox,main="Daily sales box-cox transformed",xlab = "DAILY SALES", ylab = "Frequency")
```


#### Performing stationarity test on all.
```{r,echo=FALSE}

adf.test(train_log$sale_log)
kpss.test(train_log$sale_log)
print("'
      ")

adf.test(sales$daily_Sales)
kpss.test(sales$daily_Sales)


adf.test(train_box_cox$sale_box_cox)
kpss.test(train_box_cox$sale_box_cox)
```

Both test indicates that the series is stationary
As our series is already stationary, we will move ahead with orginal series.





### FINALISING THE MODEL

#### As our series is now stationary lets try to figure out what type of series it is.

Checking for auto correlation.
```{r,echo=FALSE}
par(mfrow=c(1,2))

acf(sales$daily_Sales,main="Daily sales ACF")
pacf(sales$daily_Sales,main="Daily sales PACF")
```

The above ACF and PACF plot shows us that there is high auto correlation. The ACF has dampning effect and PACF shows no significant spike after 7. This confirms that this is AR process. The spikes in PACF are significant at level 3 and 7. so we will check for both levels.




```{r,echo=FALSE}

model_1 = arima(sales$daily_Sales,order = c(7,0,0))


summary(model_1)

```



```{r,echo=FALSE}

forecast::checkresiduals(model_1)
```

The residual are looking Good. But there seems to be litle seaonal factor left according to residual graphs. But we seems to be close to the best model.
Here we will try different models near to (3,0,0) and comapre there AIC's.

```{r,echo=FALSE}

model_2 = arima(sales$daily_Sales,order = c(3,0,0))


summary(model_2)

```


```{r,echo=FALSE}

forecast::checkresiduals(model_2)
```

If we compare the two, we see that there are still some significant spike left at lag 7 in ARIMA(3,0,0). There are no significant spikes for model ARIMA(7,0,0). The error distribution looks normal expect some large error at last. That might be the error becuase of some high seasonal value. Overall the ARIMA(7,0,0) looks good.




#### To further check if any other close model is better than ARIMA(7,0,0). We will compare the AIC's of them.
```{r,echo=FALSE}

AIC(
  arima(sales$daily_Sales,order=c(7,0,0)),
  
  arima(sales$daily_Sales,order=c(8,0,0)),
  
  arima(sales$daily_Sales,order=c(4,0,0)),
  
  arima(sales$daily_Sales,order=c(5,0,0)),
  
  arima(sales$daily_Sales,order=c(6,0,0)),
  
  arima(sales$daily_Sales,order=c(7,0,1)),
  
  arima(sales$daily_Sales,order=c(6,0,1)),
  
  arima(sales$daily_Sales,order=c(5,0,1)),
  
  arima(sales$daily_Sales,order=c(9,0,0))
  
  
  
)
```

The model looks better comapared to others. But we got some low values when we added MA component in the model. Also when we look at the residual plot, we could see there some spikes at regular intervals. This might be taken care if we add MA component. Lets try different values and compare the AIC's again.

```{r,echo=FALSE}
AIC(

  arima(sales$daily_Sales,order=c(7,0,2)),
  
  arima(sales$daily_Sales,order=c(7,0,3)),
  
  arima(sales$daily_Sales,order=c(7,0,4)),
  
  arima(sales$daily_Sales,order=c(6,0,2)),
  
  arima(sales$daily_Sales,order=c(6,0,3)),
  
  arima(sales$daily_Sales,order=c(6,0,4)),
  
  arima(sales$daily_Sales,order=c(5,0,2)),
  
  arima(sales$daily_Sales,order=c(5,0,3)),
  
  arima(sales$daily_Sales,order=c(8,0,4)),
  
  arima(sales$daily_Sales,order=c(8,0,2)),
  
  arima(sales$daily_Sales,order=c(8,0,3))
  
  
  
)
```

```{r,echo=FALSE}
AIC(

  arima(sales$daily_Sales,order=c(3,0,0)),
  
  arima(sales$daily_Sales,order=c(3,0,1)),
  
  arima(sales$daily_Sales,order=c(3,0,2)),
  
  arima(sales$daily_Sales,order=c(3,0,3)),
  
  arima(sales$daily_Sales,order=c(4,0,1)),
  
  arima(sales$daily_Sales,order=c(4,0,2)),
  
  arima(sales$daily_Sales,order=c(4,0,3)),
  
  arima(sales$daily_Sales,order=c(5,0,1)),
  
  arima(sales$daily_Sales,order=c(5,0,2)),
  
  arima(sales$daily_Sales,order=c(5,0,3)),
  
  arima(sales$daily_Sales,order=c(2,0,1))
  
)
  
```



We could see that we got the lowest AIC for ARIMA(6,0,4). Lets try to check its residual graphs.

```{r,echo=FALSE}
model_7 <- arima(sales$daily_Sales,order = c(6,0,3))

forecast::checkresiduals(model_7)
```

The AIC is reduced compare to ARIMA(7,0,0), but residual has still some spikes left. Keeping the lag at 7 in mind, we will try MA component for lag 7.

```{r,echo=FALSE}

model_3 <- arima(sales$daily_Sales,order = c(7,0,3))

forecast::checkresiduals(model_3)

```

We are still left with the spike. We will try the next best AIC  that is ARIMA (7,0,2).



```{r,echo=FALSE}

model_4 <- arima(sales$daily_Sales,order = c(7,0,2))

forecast::checkresiduals(model_4)

```

The model chosen at start ARIMA(7,0,0) is gave closely similar results to ARIMA(7,0,2). Lets try to compare them.

```{r,echo=FALSE}

summary(model_1)

summary(model_4)
```

There is not much of difference in the models, thus we will go ahead with a simpler model .i.e ARIMA(7,0,0)



We havent checked the ARIMA(3,0,0) counter parts,as we have selected ARIMA(7,0,0) now, we have options of ARIMA(3,0,0) counterparts with lower AIC's. Lets check there residuals



```{r,echo=FALSE}


model_5 <- arima(sales$daily_Sales,order = c(4,0,3))

forecast::checkresiduals(model_5)


model_6 <- arima(sales$daily_Sales,order = c(5,0,3))

forecast::checkresiduals(model_6)



```


The spike at 7 is significant in all the models except ARIMA(7,0,0). Thus we will select ARIMA(7,0,0) as our best model. This could be the result of weekly seasonality present in the data set.









#### From both AIC and BIC, the best model is ARIMA(7,0,0)



```{r,echo=FALSE}



best_fit <- arima(sales$daily_Sales,order = c(7,0,0))

forecast::checkresiduals(best_fit)

```










### AUTO ARIMA SUGGESTIONS.

#### Fitting the model
```{r,echo=FALSE}

auto_fit <- auto.arima(sales$daily_Sales,seasonal= TRUE,stepwise=FALSE,approximation=FALSE,max.p = 15,max.q = 15)


summary(auto_fit)

```




We are getting different model ARIMA(1,0,2) than our best estimation from auto.arima

#### Checking the residuals.

Performing Ljung Box test for the both best model and auto model at different lags, the Ljung-Box test has hypothesis that the residuals are independent. Thus the P value less than 0.05 indicates that there exist a correlation in the residuals.

```{r,echo=FALSE}

Box.test(best_fit$residuals,type = "Ljung-Box",lag=1)

Box.test(auto_fit$residuals,type = "Ljung-Box",lag=1)

```



```{r,echo=FALSE}

Box.test(best_fit$residuals,type = "Ljung-Box",lag=3)

Box.test(auto_fit$residuals,type = "Ljung-Box",lag=3)
```


```{r,echo=FALSE}

Box.test(best_fit$residuals,type = "Ljung-Box",lag=7)

Box.test(auto_fit$residuals,type = "Ljung-Box",lag=7)

```

```{r,echo=FALSE}

Box.test(best_fit$residuals,type = "Ljung-Box",lag=10)

Box.test(auto_fit$residuals,type = "Ljung-Box",lag=10)

```

#### The auto arima model fails the Ljung-Box test at lag seven,thus our best model will be ARIMA(7,0,0)


### IN SAMPLE RMSE

```{r,echo=FALSE}

resid_1 = best_fit$residuals

resid_2 = auto_fit$residuals

pred_1 = sales$daily_Sales-resid_1


pred_2 = sales$daily_Sales-resid_2

ggplot()+
  geom_line(aes(sales$daily_Sales,pred_1),color='red',alpha=0.7)+
  
  geom_line(aes(sales$daily_Sales,pred_2),color='blue',alpha=0.7)+
  ggtitle("BEST MODEL VS AUTO ARIMA")+
  xlab("Daily sales")+
  ylab("Prediction")
  
#  geom_line(aes(train_main_1$Week,train_main_1$week_bx))

```


```{r,echo=FALSE}


                 
RMSE_best = sqrt(mean((pred_1 - sales$daily_Sales)^2))

                 
RMSE_best

RMSE_auto = sqrt(mean((pred_2 - sales$daily_Sales)^2))

RMSE_auto


```



#### We now will calculate MAE for fitted values 
```{r,echo=FALSE}

MAE_best = mean(abs(pred_1 - sales$daily_Sales))

MAE_best


MAE_auto = mean(abs(pred_2 - sales$daily_Sales))
MAE_auto

```
#### The MAE is less for the model we choose 


```{r,echo=FALSE}

best_fit %>% 
  forecast(h=30) %>% 
  autoplot()

```


```{r,echo=FALSE}

auto_fit %>% 
  forecast(h=30) %>% 
  autoplot()

```





```{r, echo=FALSE,warning=FALSE}

#length(sales$date)


prediction = predict(best_fit,n.ahead=30)

pred_data = data.frame(
  pred = prediction,
  date = seq(c(ISOdate(2015,10,31)), by = "day", length.out = 30) #as.Date(ymd(max(sales$date)):ymd(max(sales$date)+days(29)),origin="2013-01-01")
)

summary(pred_data)

pred_merge = sales %>%
  full_join(pred_data) %>%
  filter(date>ymd("2015-09-01")) %>%  #& date<=ymd("2015-03-01")) %>%
  mutate(
    pred_high = pred.pred+2*pred.se,
    pred_low =  pred.pred - 2*pred.se,
    pred.pred = pred.pred,
    error = pred.pred - daily_Sales
    )
pred_merge %>% 
  ggplot()+
  geom_line(aes(date,daily_Sales))+
  geom_line(aes(date,pred.pred),color='blue')+
  xlab("DATE - YEAR 2015")+
  ylab("SALES")+
  geom_ribbon(aes(x=date,ymin=pred_low,ymax = pred_high),fill='blue',alpha=0.2)



```

```{r,echo=FALSE}


```


USING FB PROPHET
=======================================================================

## Starting with PROPHET{.tabset}

### CHANGEPOINTS


#### Plotting Changepoint 
```{r,echo=FALSE,include=FALSE,warning=FALSE,message=FALSE}
dim(sales)
library(changepoint)

```

Identifying change points for trend in the series.The single changepoint detected is in January 2015, where we can see that the next trend is bit downward side.

```{r,echo=FALSE}

cp = cpt.mean(sales$daily_Sales)

sales %>% 
  ggplot()+
  geom_line(aes(date,daily_Sales))+
  geom_vline(aes(xintercept=sales$date[cpts(cp)]),color='red')+
  theme_bw()+
  xlab("DATE")+
  ylab("DAILY SALES")

cpts(cp)

```

Trying to check for multiple changing point with mean values.

The multiple change points can be seen in the plot.

```{r,echo=FALSE,warning=FALSE,message=FALSE}

mult_cp=cpt.mean(sales$daily_Sales,method='BinSeg',Q=10)

plot(mult_cp)


```

### OUTLIERS
  
#### Box Plot for detecting outliers
```{r,echo=FALSE}

box=boxplot(sales$daily_Sales,main="DAILY SALES")

```

The outliers according to the boxplot above are highlighted below.
One high value above 300 looks abnormal, but it lies in peak seasonal period, thus cannot be ruled out. Keepinng as it is.


#### Plotting all the outliers from the boxplot in the series.

```{r,echo=FALSE,warning=FALSE,message=FALSE}

sales %>% 
  mutate(
    anomaly = if_else(daily_Sales %in% box$out,TRUE,as.logical(NA))
  ) %>% 
  ggplot()+
  geom_line(aes(date,daily_Sales))+
  geom_point(aes(date,daily_Sales,color=anomaly))+
  theme_bw()+
  ylim(0,400)+
  ggtitle("Outliers According to Boxplot") + 
  ylab("DAILY SALES")+
  scale_colour_discrete(na.translate = F)

```




```{#r,echo=FALSE,eval=FALSE,include=FALSE}

cleaned_data = tsclean(sales$daily_Sales)
clean_full_y = sales %>% 
  mutate(
    row_number = row_number(),
    anomaly = if_else(row_number %in% outliers$index,TRUE,as.logical(NA))
    #cleaned_data = cleaned_data
  ) %>% 
  ggplot()+
  geom_line(aes(date,cleaned_data))+
  geom_point(aes(date,cleaned_data,color=anomaly))+
  theme_bw()+
  ylim(0,400)+
  guides(color='none')+
  labs(title = 'Outliers Replaced') + 
  scale_colour_discrete(na.translate = F)

plot = clean_full_y

clean_full_y
```


```{r}
```


### FITTING THE PROPHET MODEL

#### The model is fitted on train data till date 2015-07-01 and forecast is done fo period of 120 days.
```{r,echo=FALSE,include=FALSE,warning=FALSE}

library(prophet)

prophet_sales = sales %>% 
    rename(ds = date, # Have to name our date variable "ds"
    y = daily_Sales)  # Have to name our time series "y"

train = prophet_sales %>% 
  filter(ds%
  filter(ds>=ymd("2015-07-01"))

model = prophet(train)

future = make_future_dataframe(model,periods = 120)

forecast = predict(model,future)
```


Model automatically disable daily seasonality as it detects there could be no daily seasonality present. 


```{r,echo=FALSE}

plot(model,forecast)+
ylab("TOTAL DAILY SALES")+xlab("Date")+theme_bw()
```

```{r}



```
We try check if there is any difference by turning on the daily seasonality or not.
#### Ckeck by Turning on the daily seasonality
```{r,echo=FALSE,include=FALSE,warning=FALSE}

model_1 = prophet(train,daily.seasonality = TRUE)

future_1 = make_future_dataframe(model,periods = 120)

forecast_1 = predict(model,future)


```
.


```{r,echo=FALSE,warning=FALSE,message=FALSE}

plot(model_1,forecast)+
ylab("TOTAL DAILY SALES")+xlab("Date")+theme_bw()

```

Not much of a difference here. Thus model correctly disabled the daily seasonality




#### Lets see our original series for comparison
```{r,echo=FALSE,warning=FALSE,message=FALSE}

ggplot()+
geom_point(aes(date,daily_Sales))+
  theme_bw()


```



```{r}
```

#### Interactive chart to check for perticular period. (Use scroller at bottom)

```{r,echo=FALSE,warning=FALSE,message=FALSE}

dyplot.prophet(model,forecast)

```

### CHECKING FOR COMPONENTS OF TIME SERIES

#### Plot shows the different components which includes, trend, weekly seasonality and yearly seasonality.

```{r,echo=FALSE}

prophet_plot_components(model,forecast)
```

#### Checking changepoint with phophet model

#### Plot shows that two three equidistant change points. The model includes one change point found using boxplot.
```{r,echo=FALSE,include=FALSE}

plot(model,forecast)+add_changepoints_to_plot(model)+theme_bw()+xlab("Date")+ylab("Daily Total Sales")

```


```{r,echo=FALSE,warning=FALSE,message=FALSE}

model = prophet(train,n.changepoints=5)

forecast_2 = predict(model,future)

plot(model,forecast_2)+add_changepoints_to_plot(model)+theme_bw()+xlab("Date")+ylab("Daily Sales")
```

#### COMPARING FORECAST

#### Comparing forecast with actual observed data. The model fits well but has greater variance at start.
```{r,echo=FALSE,warning=FALSE}

forecast_plot_data = forecast %>% 
  as_tibble() %>% 
  mutate(ds = as.Date(ds)) %>% 
  filter(ds>=ymd("2015-07-01"))

forecast_not_scaled = ggplot()+
geom_line(aes(test$ds,test$y))+
geom_line(aes(forecast_plot_data$ds,forecast_plot_data$yhat),color='red')+
  xlab("date")+
  ylab("Total sales")

forecast_scaled = forecast_not_scaled

forecast_scaled

```

#### Trying to forecast it for 250 days 
```{r,echo=FALSE,warning=FALSE}

two_yr_future = make_future_dataframe(model,periods = 250)
two_yr_forecast = predict(model,two_yr_future)
plot(model,two_yr_forecast)+theme_bw()+xlab("Date")+ylab("Total Sales")
```


```{r,eval=FALSE,include=FALSE,echo=FALSE}

train$floor = 0
train$cap = 25
future$floor = 0
future$cap = 25
# Set floor in forecsat data
two_yr_future$floor = 0
two_yr_future$cap = 25
#sat_model = prophet(train,growth='logistic',daily.seasonality = T)

sat_two_yr_forecast = predict(sat_model,two_yr_future)


```

.

```{r}
```


```{r,eval=FALSE,include=FALSE,echo=FALSE}

plot(sat_model,sat_two_yr_forecast)+ylim(0,25)+
  theme_bw()+xlab("Date")+ylab("Daily Sales")
```

```{r,echo=FALSE,warning=FALSE,eval=FALSE}

prophet_plot_components(model,forecast)

```

The prophet is considering seasonality correctly, We dont need to take account for daily seasonality here. We can check for additive and multiplicative seasonality for considration.

### SEASONALITY

#### Additive Seasonality 
```{r,echo=FALSE,warning=FALSE}

additive = prophet(train)

add_fcst = predict(additive,future)

plot(additive,add_fcst)

#summary(add_fcst)

#summary(train)

```


#### Plotting additive components.
```{r,echo=FALSE,warning=FALSE}

prophet_plot_components(additive,add_fcst)



```







#### Multiplicative seasonality
```{r,echo=FALSE}

multi = prophet(train,seasonality.mode = 'multiplicative')
multi_fcst = predict(multi,future)

plot(multi,multi_fcst)
```

#### Multiplicative components 
```{r,echo=FALSE}
prophet_plot_components(multi, multi_fcst)

```

#### There is not much of difference between both.

#### Checking Holiday Impact

#### Leveraging the holiday list available in the prophet model,We could see below which holidays has greater Impact.
```{r,echo=FALSE,warning=FALSE,message=FALSE}

#?prophet()

model = prophet(train,fit=FALSE)
model = add_country_holidays(model,country_name = 'US')
model = fit.prophet(model,train)
forecast = predict(model,future)
prophet_plot_components(model,forecast)

```


#### Checking the amount of impact of two main holidays. New year and Christmas 
```{r,echo=FALSE,warning=FALSE}

forecast %>% 
  filter(
    holidays!=0,
    ds == ymd("2014-12-25") | 
    ds == ymd("2015-01-01")) %>%
    select(ds,`Christmas Day`,`New Year's Day`)
```

#### Holiday Impact by each holiday in the US calender.
```{r,echo=FALSE,warning=FALSE}

forecast %>%
  filter(holidays != 0) %>%
  dplyr::select(-ds:-additive_terms_upper, -holidays:-holidays_upper, -weekly:-yhat, -contains("upper"), -contains("lower")) %>%
  mutate_all(~ if_else(. == 0, as.numeric(NA), .)) %>%
  summarize_if(is.numeric, ~ max(., na.rm = T)) %>%
  pivot_longer(
    cols = `Christmas Day`:`Washington's Birthday`,
    names_to = 'holiday', 
    values_to = 'effect'
  ) %>%
ggplot() +
  geom_col(aes(effect,holiday))+
  theme_bw()
```

### MODEL 'FIT' & CROSS VALIDATION

#### Checking In sample performance by Calculating the MSE,MAE and MAPE of the model.

#### The error values are not large, the fit is good. Generally, Prophet works good on daily data and with seasonality.
```{r,echo=FALSE,warning=FALSE}

forecast_metric_data = forecast %>% 
  as_tibble() %>% 
  mutate(ds = as.Date(ds)) %>% 
  filter(ds>=ymd("2015-07-01"))
RMSE = sqrt(mean((test$y - forecast_metric_data$yhat)^2))
MAE = mean(abs(test$y - forecast_metric_data$yhat))
MAPE = mean(abs((test$y - forecast_metric_data$yhat)/test$y))


print(paste("RMSE:",round(RMSE,2)))
print(paste("MAE:",round(MAE,2)))
print(paste("MAPE:",round(MAPE,2)))

```

#### Cross Validation

Cross- Validation for the window of 30 days.
```{r,echo=FALSE,warning=FALSE,message=FALSE,include=FALSE}

df.cv <- cross_validation(model, initial = 700, period = 30, horizon = 30, units = 'days')
head(df.cv)

```

#### Cut-off points for the cross validation
```{r,echo=FALSE}
unique(df.cv$cutoff)

```

#### Plotting the comparison of cross validation and actual for different cut-off points.
```{r,echo=FALSE}

df.cv %>% 
  ggplot()+
  geom_point(aes(ds,y)) +
  geom_point(aes(ds,yhat,color=factor(cutoff)))+
  theme_bw()+
  xlab("Date")+
  ylab("Total sales")+
  scale_color_discrete(name = 'Cutoff')

```

### Conclusion


#### This is our Final part. Lets conclude our best model after comparing in sample and out of sample performance. Here we have our best model as ARIMA(7,0,0)  which we have chosen by trial and error. The auto model is model suggested by auto.arima and last is model by using prophet library.

#### Cross validation plot for PROPHET

```{r,echo=FALSE}
plot_cross_validation_metric(df.cv, metric = 'rmse')
```

```{r,echo=FALSE}

#dim(sales)



out_list = lapply(828:1004,function(x){
  train = sales$daily_Sales[1:x]
  test = sales$daily_Sales[(x+1):(x+29)]
  train_mod1 = arima(train,order=c(7,0,0))
  train_mod2 = arima(train,order=c(1,0,2))
  
  out_data = data.frame(
    train_thru = x,
    time_ahead = 1:29,
    time_point = x + 1:29,
    actual = test,
    naive_fcst = train[x],
    test_pred1 = predict(train_mod1,29)$pred,
    test_pred2 = predict(train_mod2,29)$pred,
    error1 = test - predict(train_mod1,29)$pred,
    error2 = test - predict(train_mod2,29)$pred
  )
})

out_data = rbindlist(out_list)
```

```{r,echo=FALSE,include=FALSE}
out_data %>%
  as_tibble()
```

```{r,echo=FALSE,message=FALSE,warning=FALSE}

error_metrics = out_data %>% 
  mutate(naive_error = actual - naive_fcst) %>%
  pivot_longer(
    cols = error1:naive_error,
    values_to = 'error'
  ) %>%
  group_by(time_ahead,name) %>% 
  summarize(
    RMSE = sqrt(mean(error^2)),
    MAE = mean(abs(error)),
    MAPE = mean(abs(error / actual))
    ) %>%
    ungroup()
```
#### CROSS VALIDATION GRAPH FOR ARIMA MODELS

```{r,echo=FALSE}

error_metrics %>% 
  ggplot()+
  geom_line(aes(time_ahead,RMSE,color=name))+
  xlab("Horizon (days)")
```



#### COMPARING IN SAMPLE RMSE AND MAE
```{r,echo=FALSE}

xyz = data.frame(
  Model = c("ARIMA","AUTO ARIMA","PROPHET"),
  `In Sample RMSE` = round(c(RMSE_best,RMSE_auto,RMSE),3),
  `In Sample MAE` = round(c(MAE_best,MAE_auto,MAE),3)
  #`Out of Sample RMSE` = round(c(out_rmse1,out_rmse2,out_rmse3,out_rmse9),3)
) 
#%>%
#datatable()

xyz
  
  

```



#### LOOKING ALL THE PLOTS AND VALUES ABOVE, WE CONCLUDE THAT PROPHET HAS OUTPERFORMED ALL THE MODELS. THE PROPHET HAS PROVED ITS PREDICTIVE POWER FOR SEASONAL DATA SETS.
#### THUS OUR WINNER IS 'PROPHET' MODEL.