{"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# EDA+Prophet+ MLP Neural Network Forecasting\n\n**Arindam Dutta**\n\n22-09-2018 (Version-8)","metadata":{"_uuid":"7f592ad23986d4ebba5449c9c6486da1394a5556"}},{"cell_type":"markdown","source":"# Before we begin: ","metadata":{"_uuid":"0cb6bdedeb82326a4721cb690c15bb5e6cd8c204"}},{"cell_type":"markdown","source":"I would like to take the opportunity to appreciate the help of  Corey Levinson,Kirill Tsyganov, Alexey Kozionov, Alexander Andreev, Jaroslav Bologov, Nurlan Shagadatov & Makarychev Sergey**.\n\nIf you liked my work, then please don't forget to  upvote the Tutorial since it will keep me motivating to perform more in-depth reserach Forecasting Method.I hope you will enjoy our deep exploration into this dataset. Suggestions on improvements are welcome.\n\nLet the Forecasting begin. :)","metadata":{"_uuid":"35c4fab2fb3af76e33ffd833b279c8d047e313f8"}},{"cell_type":"markdown","source":"**What is Prophet in ML?**\n\nProphet is a procedure for forecasting time series data based on an additive model where non-linear trends are fit with yearly, weekly, and daily seasonality, plus holiday effects. It works best with time series that have strong seasonal effects and several seasons of historical data.\n\n**What type of model is Prophet?**\n\nAt its core, the Prophet procedure is an additive regression model with four main components: A piecewise linear or logistic growth curve trend. Prophet automatically detects changes in trends by selecting changepoints from the data. A yearly seasonal component modeled using Fourier series.\n","metadata":{}},{"cell_type":"markdown","source":"# Executive Summary:\n\nThis tutorial consists of comprehensive Exploratory Data Analysis with Prophet and MLP Neural Network Forecast Modeling for [Store Item Demand Forecasting Challenge](https://www.kaggle.com/c/demand-forecasting-kernels-only) competition.\n    The objective of the problem is to forecast 3 months of sales  for 50 different items at 10 different stores using the 5 years history of sales.\n    \nThe [data](https://www.kaggle.com/c/demand-forecasting-kernels-only/data) contains the following files:\n\n1. **train.csv :** the training data witht the history of sales.\n\n2. **test.csv:** the test data to predict the sales.\n\n3. **Sample submission.csv:** Sample submission file\n\n\n**Data fields**\n\ndate - Date of the sale data. There are no holiday effects or store closures.\n\nstore - Store ID\n\nitem - Item ID\n\nsales - Number of items sold at a particular store on a particular date.","metadata":{"_uuid":"8a00e7d2c1e3be485e07412c457ba310aa4faaff"}},{"cell_type":"markdown","source":" # Introduction:\nTime series forecasting is a skill that few people claim to know. Machine learning is cool. And there are a lot of people interested in becoming a machine learning expert. But forecasting is something that is a little domain specific. \n\nRetailers like **Walmart, Target** use forecasting systems and tools to replenish their products in the stores. An excellent forecast system helps in winning the other pipelines of the supply chain. If you are good at predicting the sale of items in the store, you can plan your inventory count well. You can plan your assortment well.\n\nA good forecast leads to a series of wins in the other pipelines in the supply chain. \n\n# What is time series?\nA time series is a sequence of observations collected at some time intervals. Time plays an important role here. The observations collected are dependent on the time at which it is collected.\n\nThe sale of an item say **bread** in a retail store like Walmart will be a time series. The sale could be at daily level or hourly level. The number of people flying from NewYork to Spain on daily basis is a time series. Time is important here. During Christmas holidays, this number would be humongous compared to the other days. This is know as **seasonality**. \n\n# What is the difference between a time series and a normal series?\nTime component is important here. The time series is dependent on the time. However a normal series say 1, 2, 3...100 has no time component to it. When the value that a series will take depends on the time it was recorded, it is a time series. \n\n# How to define a time series object in R\nts() function is used for equally spaced time series data, it can be at any level. Daily, weekly, monthly, quarterly, yearly or even at minutes level. If you wish to use unequally spaced observations then you will have to use other packages. \n\nts() is used for numerical observations and you can set frequency of the data. ts() takes a single frequency argument. There are times when there will be multiple frequencies in a time series. We use msts() multiple seasonality time series in such cases. I will talk about msts() in later part of the post. For now, let us define what is frequency. \n\n# Frequency\nWhen setting the frequency, many people are confused what should be the correct value. This is the simple definition of frequency. Frequency is the number of observations per cycle. Now, how you define what a cycle is for a time series? \n\nSay, you have electricity consumption at Bangalore at hourly level. The cycle could be a day, a week or even annual. I will cover what frequency would be for all different type of time series. \n\nBefore we proceed I will reiterate this.\n\n**Frequency is the number of observations per cycle.**\n\nWe will see what values frequency takes for different interval time series. \n\n**Daily data **  There could be a weekly cycle or annual cycle. So the frequency could be 7 or 365.25.\n\nSome of the years have 366 days (leap years). So if your time series data has longer periods, it is better to use frequency = 365.25. This takes care of the leap year as well which may come in your data. \n\n**Weekly data** There could be an annual cycle. frequency = 52 and if you want to take care of leap years then use frequency = 365.25/7 \n\n**Monthly data** Cycle is of one year. So frequency = 12 \n\n**Quarterly data** Again cycle is of one year. So frequency = 4 \n\n**Yearly data** Frequency = 1 \n\nThere could be smaller interval freuqncy as well like **Hourly, Half-Hourly, Mintues, Seconds** time series data as well\n","metadata":{"_uuid":"8329dabaf921b3c459abb4f5af8aba2fc2000fc8"}},{"cell_type":"markdown","source":"# Table of Contents\n\n## 1. [Load Libraries](#section1)\n## 2. [Load Data  & Structure of the Data](#section2)\n## 3. [Missing Values Detection](#section3)\n## 4. [Individual feature visualisations](#section4)\n>## 4. 1.[Histogram of Sale Price](#section4.1) \n>## 4.2.[Growth by date](#section4.2) \n>## 4.3.[Growth by Month of Different Year](#section4.3) \n>## 4.4.[Growth by Year](#section4.4) \n>## 4.5.[Growth By Store](#section4.5) \n>## 4.6.[ Yearly Growth By Item](#section4.6) \n  \n \n  \n## 5. [PART 1-What is the Prophet model?](#section5)\n>## 5. 1.[Background](#section5.1) \n>## 5.2.[Advantages over other Time Series Models](#section5.2) \n>## 5.3.[Modeling Framework & Theory](#section5.3) \n> ## 5.4.[Model building using PROPHET: store= 1, Product_ID=1](#section5.4) \n>## 5.5.[Customizing holidays and events](#section5.5) \n> ## 5.6.[IncludeFourier Order of Seasonality &  Additional Regressors](#section5.6) \n> ## 5.7.[ Including Additional Regressors e.g. NFL Sundays](#section5.7) \n> ## 5.8.[SMAPE Calculation](#section5.8) \n >## 5.9.[Automation for Prophet: Splitting data by Store and Item](#section5.9) \n > ## 5.10.[Final Prediction using Prophet](#section5.10) \n > ## 5.11.[How to Improve Propehet Model Further](#section5.11) \n  \n\n\n## 6.[PART 2: Forecasting Time Series with Neural Networks](#section6)\n>  ## 6. 1.[Background](#section6.1) \n>  ## 6.2.[Baseline Predcition with 1 Hidden Layer](#section6.2) \n>  ## 6.3.[Plotting the Forecast & Calculate the MSE (Mean Squared Error)](#section6.3) \n>  ## 6.4.[Shifting Input](#section6.4) \n>  ## 6.5.[Outplot & Model Fitting](#section6.5) \n\n\n\n## 7. [Extensions](#section7)\n>Choosing Oprimal No of Hidden Layer based on min(MSE).\n\n>Fitting ELM Feed Forward Network for forecasting.\n\n>Temporal Hierarchies forecasting\n\n>Automation based on MLP/ ELM NN Forecasting.\n\n\n\n## 8. [References](#section8)","metadata":{"_uuid":"ceee6dd7610969ac0a3ea55c08e84dd7ea7091fc"}},{"cell_type":"markdown","source":"<a id=\"section1\"></a>\n# 1.LOAD LIBRARIES:","metadata":{"_uuid":"eeacbc04b090cc226874be09d5dbe7b3d478c8cf"}},{"cell_type":"code","source":"rm(list=ls())\n\nsuppressMessages(library(data.table))\nsuppressMessages(library(DT))\nsuppressMessages(library(timeSeries))\nsuppressMessages(library(tidyverse))\nsuppressMessages(library(reshape))\nsuppressMessages(library(stringr))\nsuppressMessages(library(doBy))\nsuppressMessages(library(formattable))\nsuppressMessages(library(gridExtra))\n\nsuppressMessages(library(ggplot2))\nsuppressMessages(library(plotly))\nsuppressMessages(library(corrplot))\nsuppressMessages(library(wesanderson))\nsuppressMessages(library(RColorBrewer))\nsuppressMessages(library(gridExtra))\nsuppressMessages(library(zoo))\n\nsuppressMessages(library(forecast))\nsuppressMessages(library(prophet)) ### For Prophet Forecasting\nsuppressMessages(library(nnfor))    ### For Neural Network Forecasting\n\nset.seed(2018)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":false,"_kg_hide-input":false,"execution":{"iopub.status.busy":"2022-07-09T18:50:41.738888Z","iopub.execute_input":"2022-07-09T18:50:41.741478Z","iopub.status.idle":"2022-07-09T18:50:47.566406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" <a id=\"section2\"></a>\n # 2. Load Data , Structure & Summary of the Data","metadata":{"_uuid":"cd0aadfa9dcf68cd389c1fb06ac4be5a9c88a670","trusted":true}},{"cell_type":"code","source":"train=fread(\"../input/train.csv\")\nsprintf(\"The train data set has %d rows and %d columns\", nrow(train), ncol(train) )\nstr(train)\n\ntest  <- fread(\"../input/test.csv\")\nsprintf(\"The test data set has %d rows and %d columns\", nrow(test), ncol(test) )\nstr(test)\n\nprint(\"the summary of train sales is:\")\nsummary(train$sales)\n\n\n# Extraction of Year and Month of Year :\ntrain$Year=year(train$date)        #returns the year from date i.e. 2013, 2014 etc.\ntrain$Month=as.yearmon(train$date) #this yearmon() function is coming from zoo package returns the month of an year i.e Jan 2013, Feb 2015 etc","metadata":{"_uuid":"5742399a757bb704e128c6a3983b00a2d45b4dc2","execution":{"iopub.status.busy":"2022-07-09T18:51:11.225213Z","iopub.execute_input":"2022-07-09T18:51:11.258545Z","iopub.status.idle":"2022-07-09T18:51:30.781745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"section3\"></a>\n# 3. Missing Values Detection:","metadata":{"_uuid":"86bd882e1e6872572dbddacdb59b4227ccceec18","trusted":true}},{"cell_type":"code","source":"\ncolSums(is.na(train))\n# Function 1 : For ploting missing value\nplot_missing <- function(data, title = NULL, ggtheme = theme_gray(), theme_config = list(\"legend.position\" = c(\"bottom\"))) {\n  ## Declare variable first to pass R CMD check\n  feature <- num_missing <- pct_missing <- group <- NULL\n  ## Check if input is data.table\n  is_data_table <- is.data.table(data)\n  ## Detect input data class\n  data_class <- class(data)\n  ## Set data to data.table\n  if (!is_data_table) data <- data.table(data)\n  ## Extract missing value distribution\n  missing_value <- data.table(\n    \"feature\" = names(data),\n    \"num_missing\" = sapply(data, function(x) {sum(is.na(x))})\n  )\n  missing_value[, feature := factor(feature, levels = feature[order(-rank(num_missing))])]\n  missing_value[, pct_missing := num_missing / nrow(data)]\n  missing_value[pct_missing < 0.05, group := \"Good\"]\n  missing_value[pct_missing >= 0.05 & pct_missing < 0.4, group := \"OK\"]\n  missing_value[pct_missing >= 0.4 & pct_missing < 0.8, group := \"Bad\"]\n  missing_value[pct_missing >= 0.8, group := \"Remove\"][]\n  ## Set data class back to original\n  if (!is_data_table) class(missing_value) <- data_class\n  ## Create ggplot object\n  output <- ggplot(missing_value, aes_string(x = \"feature\", y = \"num_missing\", fill = \"group\")) +\n    geom_bar(stat = \"identity\") +\n    geom_text(aes(label = paste0(round(100 * pct_missing, 2), \"%\"))) +\n    scale_fill_manual(\"Group\", values = c(\"Good\" = \"#1a9641\", \"OK\" = \"#a6d96a\", \"Bad\" = \"#fdae61\", \"Remove\" = \"#d7191c\"), breaks = c(\"Good\", \"OK\", \"Bad\", \"Remove\")) +\n    scale_y_continuous(labels = comma) +\n    coord_flip() +\n    xlab(\"Features\") + ylab(\"Number of missing rows\") +\n    ggtitle(title) +\n    ggtheme + theme_linedraw()+\n    do.call(theme, theme_config)\n  ## Print plot\n  print(output)\n  ## Set return object\n  return(invisible(missing_value))\n}\nplot_missing(train)","metadata":{"_uuid":"ee7fc9fa880fc7b9f482dd334e13199262ab3c9f","execution":{"iopub.status.busy":"2022-07-09T18:51:46.532045Z","iopub.execute_input":"2022-07-09T18:51:46.534405Z","iopub.status.idle":"2022-07-09T18:51:48.306323Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> **Therefore there is NO MISSING VALUES in our data**","metadata":{"_uuid":"ee687e85d258c4214c9b63f835bde27cf004a67e","trusted":true}},{"cell_type":"markdown","source":"<a id=\"section4\"></a>\n# 4. Individual feature visualisations:","metadata":{"_uuid":"a59ef789b94e749b94bbb9d4c4de9ba39a218455"}},{"cell_type":"markdown","source":"<a id=\"section4.1\"></a>\n# 4.1. Histogram of Sale Price","metadata":{"_uuid":"793e2e14dcf6d5a0c0caedfc4aaf7dce7089f5c7","trusted":true}},{"cell_type":"code","source":"gbp1<-wes_palette(\"GrandBudapest2\")[1]\n\nggplot(train, aes(x=sales))+\n  geom_histogram(fill=\"#a6d96a\", alpha=.9)+\n  labs(x=NULL, y=NULL, title = \"Histogram of Sale Price\")+\n  # scale_x_continuous(breaks= seq(0,600000, by=100000))+\n  theme_minimal() + theme(plot.title=element_text(vjust=3, size=15) )","metadata":{"_uuid":"9aa89255ce2cc7563d81f0cf48ef398c8f724a92","execution":{"iopub.status.busy":"2022-07-09T18:52:20.703569Z","iopub.execute_input":"2022-07-09T18:52:20.705833Z","iopub.status.idle":"2022-07-09T18:52:22.485706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Therefore our Sales Price follows Positively Skewed Distribution.**","metadata":{"_uuid":"af4cd1c40f463191433d56a68aac7a4c27cf7372","trusted":true}},{"cell_type":"markdown","source":"> **Let's check the Growth & Growth Rate by Date****","metadata":{"_uuid":"0b0ec439705730f0e905947d10285eb58b63a51b","trusted":true}},{"cell_type":"markdown","source":"<a id=\"section4.2\"></a>\n# 4.2.Growth by date","metadata":{"_uuid":"9cbb2f615bd74ca16c6c8f599e7031e8bdb014c1","trusted":true}},{"cell_type":"code","source":"MSP <- aggregate(sales ~date, train, mean)\n# MSP <-na.omit(ddply(data, 'date', summarise, mean(Sale_Prices, na.rm=T)))\n\nsl1 <-ggplot(MSP, aes(x=as.factor(date), y=sales))+\n  geom_line(color=gbp1, aes(group=1), size=1.5)+\n  geom_point(colour=gbp1, size = 3.5, alpha=0.5)+\n  labs(title=\"The Growth of Sale Prices by date\", x=NULL, y=\"Sale Price\")+\n  theme( plot.title=element_text(vjust=3, size=15) ) + theme_minimal()\n\nMSP$rate = c(0, 100*diff(MSP$sales)/MSP[-nrow(MSP),]$sales)\n\nsl2 <-ggplot(MSP, aes(x=as.factor(date), y=rate))+\n  geom_line(color= \"gray50\", aes(group=1), size=1)+\n  #geom_point(colour=gbp1, size = 3.5, alpha=0.5)+\n  labs(title=\"Change rate of Sale Price\", x=\"date\", y=\"rate of change\")+\n  geom_hline(yintercept = 0, color = gbp1 )+\n  theme(plot.title=element_text(size=15))+ theme_minimal()\n\ngrid.arrange(sl1,sl2)\n","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","execution":{"iopub.status.busy":"2022-07-09T18:52:44.728819Z","iopub.execute_input":"2022-07-09T18:52:44.730420Z","iopub.status.idle":"2022-07-09T18:52:48.126472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"1. *The Growth of the Sales Price are Multiplicative with incresing TREND and SEASONALITY.*\n\n2. *The Change in Rate of Sales Price is looking constant by Date. But is the Growth Rate Change TRUE at Yearly / Monthly Level. Lets check it out*","metadata":{"_uuid":"45bc99f803b5fbd34d364e10e12512b4887abdaa"}},{"cell_type":"markdown","source":"**########### Let's check the Growth of Sales Price by Month of different Year #################**","metadata":{"_uuid":"8edb0319a94a4af59f99a8a38170b4954d6ebeb3"}},{"cell_type":"markdown","source":"<a id=\"section4.3\"></a>\n# 4.3. Growth by Month of Different Year","metadata":{"_uuid":"727c559801a8896a92c42d122ddf5f94b607fd55"}},{"cell_type":"code","source":"\nMSP <- aggregate(sales ~Month, train, mean)\n# MSP <-na.omit(ddply(data, 'date', summarise, mean(Sale_Prices, na.rm=T)))\n\nsl1 <-ggplot(MSP, aes(x=as.factor(Month), y=sales))+\n  geom_line(color=gbp1, aes(group=1), size=1.5)+\n  geom_point(colour=gbp1, size = 3.5, alpha=0.5)+\n  labs(title=\"The Growth of Sale Prices by Month of Year\", x=NULL, y=\"Sale Price\")+\n  theme( plot.title=element_text(vjust=3, size=15) ) + theme_minimal()\n\nMSP$rate = c(0, 100*diff(MSP$sales)/MSP[-nrow(MSP),]$sales)\n\nsl2 <-ggplot(MSP, aes(x=as.factor(Month), y=rate))+\n  geom_line(color= \"gray50\", aes(group=1), size=1)+\n  #geom_point(colour=gbp1, size = 3.5, alpha=0.5)+\n  labs(title=\"Change rate of Sale Price\", x=\"Month\", y=\"rate of change\")+\n  geom_hline(yintercept = 0, color = gbp1 )+\n  theme(plot.title=element_text(size=15))+ theme_minimal()\n\ngrid.arrange(sl1,sl2)\n","metadata":{"_uuid":"1a8d3e7a9fb0c33f6a3e87f685ab07b2bcb48591","execution":{"iopub.status.busy":"2022-07-09T18:53:40.391906Z","iopub.execute_input":"2022-07-09T18:53:40.394420Z","iopub.status.idle":"2022-07-09T18:53:49.604727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"** Therefore, our assumption is quite TRUE that the sales prices are multiplicative with increasing TREND & SEASONALITY.\n** \n\n**But there is a change in scenario for the monthly rate of change of sales price. The rate is almost constant at daily level. But it is highly fluctuating at monthly Level.**","metadata":{"_uuid":"cfef473160347b918162ccb0cb2402aa8ad206f4"}},{"cell_type":"markdown","source":"<a id=\"section4.4\"></a>\n# 4.4. Growth by Year","metadata":{"_uuid":"9c5959644095f05e420738ce5de7a32f0728d9a2"}},{"cell_type":"code","source":"MSP <- aggregate(sales ~Year, train, mean)\n# MSP <-na.omit(ddply(data, 'date', summarise, mean(Sale_Prices, na.rm=T)))\n\nsl1 <-ggplot(MSP, aes(x=as.factor(Year), y=sales))+\n  geom_line(color=gbp1, aes(group=1), size=1.5)+\n  geom_point(colour=gbp1, size = 3.5, alpha=0.5)+\n  labs(title=\"The Growth of Sale Prices by Year\", x=NULL, y=\"Sale Price\")+\n  theme( plot.title=element_text(vjust=3, size=15) ) + theme_minimal()\n\nMSP$rate = c(0, 100*diff(MSP$sales)/MSP[-nrow(MSP),]$sales)\n\nsl2 <-ggplot(MSP, aes(x=as.factor(Year), y=rate))+\n  geom_line(color= \"gray50\", aes(group=1), size=1)+\n  #geom_point(colour=gbp1, size = 3.5, alpha=0.5)+\n  labs(title=\"Change rate of Sale Price\", x=\"Year\", y=\"rate of change\")+\n  geom_hline(yintercept = 0, color = gbp1 )+\n  theme(plot.title=element_text(size=15))+ theme_minimal()\n\ngrid.arrange(sl1,sl2)","metadata":{"_uuid":"b96b88990693149d6a3584cc252cdda6a78260dc","execution":{"iopub.status.busy":"2022-07-09T18:54:14.871336Z","iopub.execute_input":"2022-07-09T18:54:14.874510Z","iopub.status.idle":"2022-07-09T18:54:15.970857Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" 1.The Growth of Sales Price is increasing at Yearly  level.\n \n 2.But there is a change in scenario for Rate of Change of Sales Price by Year.\n \n 3.From 2013 to 2014 the rate has monotonically increased, after thaat there is a drop in rate till 2015. \n \n 4.ALthough the rate increased upto 2016, but the highest incre occur in 2014.\n \n 5.Therefore we can expect that in 2018 the rate of change of sales will increse if it follows the simililar pattern.\n ","metadata":{"_uuid":"e076ea320568bab00529d21b5399861033612a20","trusted":true}},{"cell_type":"markdown","source":"<a id=\"section4.5\"></a>\n# 4.5 Growth By Store","metadata":{"_uuid":"0bf73a8c32b9d2f52cf646760f8ccaee7f621fe2"}},{"cell_type":"code","source":"unique(train$store)\nYear_state<-aggregate(sales ~store+Year, train,mean)\npal<-rep(brewer.pal(10, \"BrBG\"),5)\n\nggplot(Year_state, aes(group = store ))+\n  geom_line(aes(x=Year,y=sales,color=store), alpha=0.5, show.legend=F)+\n  labs(title=\"The Growth of Sales Price by Store from 2013 - 2017\", x=NULL\n  )+\n  theme(panel.background=element_rect(fill = \"Black\"),\n        plot.title=element_text(vjust=3, size=15),\n        panel.grid.major=element_line(color = pal))","metadata":{"_uuid":"88555c9c1f81b0d6b39bdf2010bebfbf5201efb9","execution":{"iopub.status.busy":"2022-07-09T18:54:44.083757Z","iopub.execute_input":"2022-07-09T18:54:44.089308Z","iopub.status.idle":"2022-07-09T18:54:44.969984Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Therefore there are 10 unique stores(1-10) and the Store-3 has the highest Yearly Growth of Sales whereas Store-7 has the lowest.**","metadata":{"_uuid":"8b6c28ddd961270add402a007d1d709316c9da01"}},{"cell_type":"markdown","source":"<a id=\"section4.6\"></a>\n# 4.6. Yearly Growth By Item","metadata":{"_uuid":"53b286f0ec66473bcdaa2721a1833c8774cc8443"}},{"cell_type":"code","source":"unique(train$item)\nYear_state<-aggregate(sales ~item+Year, train,mean)\npal<-rep(brewer.pal(10, \"BrBG\"),5)\n\nggplot(Year_state, aes(group = item ))+\n  geom_line(aes(x=Year,y=sales,color=item), alpha=0.5, show.legend=F)+\n  labs(title=\"The Growth of Sales Price by Store from 2013 - 2017\", x=NULL\n  )+\n  theme(panel.background=element_rect(fill = \"Black\"),\n        plot.title=element_text(vjust=3, size=15),\n        panel.grid.major=element_line(color = pal))","metadata":{"_uuid":"60225a92fca5c079d38d689cee4b011c0c6d78b6","execution":{"iopub.status.busy":"2022-07-09T18:55:06.665588Z","iopub.execute_input":"2022-07-09T18:55:06.672874Z","iopub.status.idle":"2022-07-09T18:55:08.238032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Therefore there are 50 unique stores(1-50) and the ITEM-17 has the highest Yearly Growth of Sales whereas Store-26 has the lowest.**","metadata":{"_uuid":"db603e1bdfe419489e0654fd6a14ed04587c086d"}},{"cell_type":"markdown","source":"<a id=\"section5\"></a>\n# 5. What is the Prophet model?","metadata":{"_uuid":"c5889bec4a7ff48853d3be512cc8c2fd23ae8164"}},{"cell_type":"markdown","source":"<img src=\"https://raw.githubusercontent.com/uselessskills/auto_ml/master/Prophet_analysis/pngs/final_formula_new.png\" width=\"100%\">","metadata":{"_uuid":"7f83ee57cf71f6946ceaff1f39e1949917c2ee4f"}},{"cell_type":"markdown","source":"Prophet model has 5 training parameters: \n- $k$ - base trend, \n- $m$ - offset parameter, \n- $\\delta = \\{\\delta_i\\}_{i = 0}^{\\#changepoints}$ - changing of trend in changepoints, \n- $\\beta = \\{ \\beta_i\\}_{i = 0}^{Z}$ - parameters for Fourier series (seasonal component) and extra regressors (exogenous features), where $Z = 2 * \\sum_{i: \\text{seasonal-features}}\\left(\\text{Seasonal-order}(i) \\right)+ \\# \\{\\text{exogenous regressors}\\}$, \n- $\\sigma$ - level of noise,\n- $mu_{\\text{train}}$, $std_{\\text{train}}$ - mean and std value, calculated for each the $i$-th regressor on train. \n\nFormula above doesn't take into account some detailes, for example that time is standartized.\n","metadata":{"_uuid":"9de7bd0e5eb7f19906bfc7ee548c9d3ac5d4c1b9"}},{"cell_type":"markdown","source":"  <a id=\"section5.1\"></a>\n  # 5.1. Background: ","metadata":{"_uuid":"aa896f726255afae8cd8e01e76f18f3ca2c97f94"}},{"cell_type":"markdown","source":"**1. When you want to forecast the time series data in R, you typically would use a package called ‘forecast’, with which you can use models like ARIMA.\n2. But then, beginning of 2017, a team at Facebook released ‘Prophet’, which utilizes a Bayesian based curve fitting method to forecast the time series data.\n**","metadata":{"_uuid":"bd7cc2b23655fae4162a33faa6cb2c7e614555cf"}},{"cell_type":"markdown","source":"<a id=\"section5.2\"></a>\n# 5.2.Advantages over other Time Series Models:","metadata":{"_uuid":"1064a2c176b777bcefc87ba72e303f3b79c2eca6"}},{"cell_type":"markdown","source":"**a.\tThe cool thing about Prophet is that it doesn’t require much prior knowledge or experience of forecasting time series data since it automatically finds seasonal trends beneath the data and offers a set of ‘easy to understand’ parameters.  Hence, it allows non-statisticians to start using it and get reasonably good results that are often equal or sometimes even better than the ones produced by the experts.**\n\n**b.\tProphet modelling can be able to detect the Change Points in time series data.**\n\n**c.\tWe can include the holidays (play-offs & super-bowls) in our data. Details has been added later.**\n\n**d.\tWe can regularise the parameters by means of Bayesian oprimisation with cross-validation.**\n\n**e.\tWe can incorporate the multiplicative-seasonality and determine the uncertainty intervals in the data.**\n\n**f.\tAdditional regressors can be added to the linear part of the model using the add_regressor method or function. A column with the regressor value will need to be present in both the fitting and prediction data-frames.**\n\n","metadata":{"_uuid":"582449cbc8d2f9bf240693cd8ccfdaa5a429399f"}},{"cell_type":"markdown","source":"<a id=\"section5.3\"></a>\n# 5.3. Modeling Framework & Theory:","metadata":{"_uuid":"91827b972fb34239347840cd9b1a47894cac040c"}},{"cell_type":"markdown","source":"<img src=\"https://raw.githubusercontent.com/uselessskills/auto_ml/master/Prophet_analysis/pngs/stan.png\" width=\"100%\">\n\nWhere by $s_m$ and $s_a$ we denote coefficients (fourier coefficients and standartized regressor values) for multiplicative and additive terms respectively.","metadata":{"_uuid":"0b246616a92c5486c8fb0ffcd88e213c4e317f1a"}},{"cell_type":"markdown","source":"<img src=\"https://raw.githubusercontent.com/uselessskills/auto_ml/master/Prophet_analysis/pngs/changepoints_process.png\" width=\"100%\">\nThis image was taken from original Facebook Prophet documentation, see https://facebook.github.io/prophet/docs/trend_changepoints.html","metadata":{"_uuid":"ce82cfd24808c9721489a728078b803bc4cacfb1"}},{"cell_type":"markdown","source":"<img src=\"https://raw.githubusercontent.com/uselessskills/auto_ml/master/Prophet_analysis/pngs/changepoints_prior_scale.png\" width=\"100%\">\nThis image was taken from original Facebook Prophet documentation, see https://facebook.github.io/prophet/docs/trend_changepoints.html","metadata":{"_uuid":"371377cba494df8e50405d19be542ab593cd10fb"}},{"cell_type":"markdown","source":"Parameter **changepoint_prior_scale** affects only in MCMC mode. \nChangepoint_prior_scale is used as initial parameter $\\tau$ for prior for $delta$  ~ $Laplace(0, \\tau)$. Stan optimization uses priors for optimization procedure, that obtains estimation on parameters. \n\nIn the case when we set changepoint_prior_scale is equal to 0.5, we have more flexibility for prior for $\\delta$ in opposite to changepoint_prior_scale=0.001 . And this will influence on estimation of this parameter. This effect we see on the plot.","metadata":{"_uuid":"864f80d320a74df18ad6a919e59463ff24ffbf2c"}},{"cell_type":"markdown","source":"<img src=\"https://raw.githubusercontent.com/uselessskills/auto_ml/master/Prophet_analysis/pngs/trend_middle.png\" width=\"100%\">","metadata":{"_uuid":"6cb1e8f6d2faf60fa661d3ce9f486add326980da"}},{"cell_type":"markdown","source":"<img src=\"https://raw.githubusercontent.com/uselessskills/auto_ml/master/Prophet_analysis/pngs/multiplicative_terms_new.png\" width=\"100%\">","metadata":{"_uuid":"15414689812462e5ce3c047e6060c17f3f0b1966"}},{"cell_type":"markdown","source":"<img src=\"https://raw.githubusercontent.com/uselessskills/auto_ml/master/Prophet_analysis/pngs/additive_terms_new.png\" width=\"100%\">","metadata":{"_uuid":"1f79cb8a41263fcbfe98a0f05c56a84f7ee21b27"}},{"cell_type":"markdown","source":"Let us denote by $r_1$, $r_2$ standartized value for regressor $f_1$ anf $f_2$ respectively.\n$\\phi_{sin_Y}$, $\\phi_{cos_Y}$ denotes the fourier coefficient for yearly seasonal component.\n$\\phi_{sin_M}$, $\\phi_{cos_M}$ denotes the fourier coefficient for monthly seasonal component.\n\n- **Then the multiplicative term in this case will be just $\\phi_{sin_Y} * \\beta_{0} + \\phi_{cos_Y} * \\beta_{1}$ + $\\phi_{sin_M} * \\beta_{2} + \\phi_{cos_M} * \\beta_{3} + r_1 * \\beta_4$**\n\n- **And additive will be $r_2 \\beta_5$**\n\n<img src=\"https://github.com/uselessskills/auto_ml/blob/master/Prophet_analysis/pngs/additive_product_new.png?raw=true\" width=\"100%\">","metadata":{"_uuid":"e23f6704a1d83316c2b5f1f5bd46eb42cce27022"}},{"cell_type":"markdown","source":"<a id=\"section5.4\"></a>\n# 5.4. Model building using PROPHET:  store= 1, Product_ID=1","metadata":{"_uuid":"9f48c634afa6db0125272f87dbee81f93424c005"}},{"cell_type":"markdown","source":"\n> I have created a baseline model without optimising the prophet parameters. To do this I have taken the time series based on logarithmic sales.","metadata":{"_uuid":"e8d4fff34294d5600f159f036a4a404bad94bf34"}},{"cell_type":"code","source":"train_final_store1_item1=subset(train,train$store==1 & train$item==1)\n\nstats=data.frame(y=log1p(train_final_store1_item1$sales)\n                 ,ds=train_final_store1_item1$date)\nstats=aggregate(stats$y,by=list(stats$ds),FUN=sum)\nhead(stats)\ncolnames(stats)<- c(\"ds\",\"y\")\n\nmodel_prophet = prophet(stats)\nsummary(model_prophet)\nfuture = make_future_dataframe(model_prophet, periods = 90)\nforecast = predict(model_prophet, future)","metadata":{"_uuid":"864b8308bf15ec1837633d4a79313e622ec6791f","execution":{"iopub.status.busy":"2022-07-09T18:56:48.789525Z","iopub.execute_input":"2022-07-09T18:56:48.791704Z","iopub.status.idle":"2022-07-09T18:56:59.524013Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's  write a function to visualise the **CHANGE POINTS** in our prophet model. ","metadata":{"_uuid":"6939d921307b11ba2dfcbe1eca75d09f7ea72555"}},{"cell_type":"code","source":"add_changepoints_to_plot <- function(m, threshold = 0.01, cp_color = \"red\",\n                               cp_linetype = \"dashed\", trend = TRUE, ...) {\n  layers <- list()\n  if (trend) {\n    trend_layer <- ggplot2::geom_line(\n      ggplot2::aes_string(\"ds\", \"trend\"), color = cp_color, ...)\n    layers <- append(layers, trend_layer)\n  }\n  signif_changepoints <- m$changepoints[abs(m$params$delta) >= threshold]\n  cp_layer <- ggplot2::geom_vline(\n    xintercept = as.integer(signif_changepoints), color = cp_color,\n    linetype = cp_linetype, ...)\n  layers <- append(layers, cp_layer)\n  return(layers)\n}\nplot(model_prophet, forecast)+ add_changepoints_to_plot(model_prophet)","metadata":{"_uuid":"855ed7e24fe826b9a60b1d6cdbb2869fe1292749","execution":{"iopub.status.busy":"2022-07-09T18:57:07.594329Z","iopub.execute_input":"2022-07-09T18:57:07.595321Z","iopub.status.idle":"2022-07-09T18:57:08.268977Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Inspecting Model Components:**","metadata":{"_uuid":"7cb780a82b89a9e778c92d2fcc53a6a1d8ec028f"}},{"cell_type":"markdown","source":"Such a bad baseline forecastong on train data.Isn't it!!\n\nModel is completely overfitting and there are so many changing points( Points marked by the RED LINES)  which we need to remove.Let’s inspect the model components:\n","metadata":{"_uuid":"a501832b787dce9e0ad21ef653277c356d26fa77"}},{"cell_type":"code","source":"prophet_plot_components(model_prophet, forecast)","metadata":{"_uuid":"80e25434aa54c4b00a44c1c9d8e3e3ce6f536c75","execution":{"iopub.status.busy":"2022-07-09T18:57:21.077151Z","iopub.execute_input":"2022-07-09T18:57:21.078616Z","iopub.status.idle":"2022-07-09T18:57:22.572767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"** Therefore we can see there is a Drop of Sales from Sunday to Monday. Thefore there MUST be an HOLIDAY effect on our sales data. \n There is peak in sales in July that means those may the festive times or seasonal sales with high discount prices.\n Let's optimise our Prophet Parameters, try to Exclude Change Points and inlcude Holiday Effects with Additional Regressors.**","metadata":{"_uuid":"20415479cd72ef57d4eedd1696dbeb07f2130a7b"}},{"cell_type":"markdown","source":"<a id=\"section5.5\"></a>\n# 5.5. Customizing holidays and events:","metadata":{"_uuid":"4b2faf41652cb954ce071b5cb5ff836ba13f72ea"}},{"cell_type":"markdown","source":" ######Holiday Effects########\nIf you have holidays or other recurring events that you’d like to model, you must create a dataframe for them. \nIt has two columns (holiday and ds) and a row for each occurrence of the holiday.\nIt must include all occurrences of the holiday, both in the past (back as far as the historical data go) and in the future (out as far as the forecast is being made).\nIf they won’t repeat in the future, Prophet will model them and then not include them in the forecast.\nYou can also include columns lower_window and upper_window which extend the holiday out to [lower_window, upper_window] days around the date. For instance, if you wanted to included Christmas Eve in addition to Christmas you’d include lower_window=-1, upper_window=0. If you wanted to use Black Friday in addition to Thanksgiving, you’d include lower_window=0, upper_window=1. \nYou can also include a column prior_scale to set the prior scale separately for each holiday, as described below.\n\nThere are two types of Holidays:\n\n 1. **Playoffs**: These are the Less-Important public holidays and weekends.\n 2. **Superbowls**: These are the Festive Holidays with high importance. There may occurs a high increase or decrease in sales due to these holidays e.g. New Year, Christmas etc.\n","metadata":{"_uuid":"3d3f01b805c3e5aad588d24454820c8bc698a49b"}},{"cell_type":"markdown","source":" I have inlcuded the Holiday Sales in the month of July for different Years in Playoffs.","metadata":{"_uuid":"340e58e239c12ea0c175a740a017e79490ef2ae6"}},{"cell_type":"code","source":"playoffs <- data_frame(\n  holiday = 'playoff',\n  ds = as.Date(c('2013-07-12', '2014-07-12', '2014-07-19',\n                 '2014-07-02', '2015-07-11', '2016-07-17',\n                 '2016-07-24', '2016-07-07','2016-07-24')),\n  lower_window = 0,\n  upper_window = 1\n)\n","metadata":{"_uuid":"506573cd1527e97545fe59c8b64b4372864c7d39","execution":{"iopub.status.busy":"2022-07-09T18:57:46.810842Z","iopub.execute_input":"2022-07-09T18:57:46.814998Z","iopub.status.idle":"2022-07-09T18:57:46.830201Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I have inlcuded the Holiday Sales for Festive seasons like New Year & Christmas for different Years in Superbowls.","metadata":{"_uuid":"6725da378d23a32e810e68234513a483d40536d4"}},{"cell_type":"code","source":"superbowls <- data_frame(\n  holiday = 'superbowl',\n  ds = as.Date(c('2013-01-01', '2013-12-25', '2014-01-01', '2014-12-25','2015-01-01', '2015-12-25','2016-01-01', '2016-12-25',\n                '2017-01-01', '2017-12-25')),\n  lower_window = 0,\n  upper_window = 1\n)\nholidays <- bind_rows(playoffs, superbowls)","metadata":{"_uuid":"50f8651ff75c7d031a9e4cafb7e4f2fe78ab8f6e","execution":{"iopub.status.busy":"2022-07-09T18:58:02.969594Z","iopub.execute_input":"2022-07-09T18:58:02.971682Z","iopub.status.idle":"2022-07-09T18:58:02.997442Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"section5.6\"></a>\n# 5.6.IncludeFourier Order of Seasonality &  Additional Regressors","metadata":{"_uuid":"2e73fa0c762b8c779901e9b71086fa5afa6e65f4"}},{"cell_type":"markdown","source":"Seasonality is estimated using partial Fourier sum which can approximate an arbitrary periodic signal. \nThe no of  terms in the partial sum (the order) is a parameter that determines how quickly the seasonality can change. \nDefault is 10, but using Bayesian optimisation  with cross validation, I have  taken fourier order=5.\nProphet will by default fit weekly and yearly seasonalities, if the time series is more than two cycles long. \nIt will also fit daily seasonality for a sub-daily time series. You can add other seasonalities (monthly, quarterly, hourly) using the add_seasonality in R.\n\nThe inputs to this function are a name, the period of the seasonality in days, and the Fourier order for the seasonality. \nFor reference, by default Prophet uses a Fourier order of 3 for weekly seasonality and 10 for yearly seasonality. \nAn optional input to add_seasonality is the prior scale for that seasonal component.\n\nAdditional regressors can be added to the linear part of the model using the add_regressor method or function. A column with the regressor value will need to be present in both the fitting and prediction dataframes.\nI have added an additional effect on Sundays during the NFL season (National Football League ).\nThe add_regressor function provides a more general interface for defining extra linear regressors, and does not require that the regressor be a binary indicator. Another time series ould be used as a regressor, although its future values would have to be known.\nThe extra regressor must be known for both the history and for future dates. It thus must either be something that has known future values (such as nfl_sunday), or something that has separately been forecasted elsewhere. Prophet will also raise an error if the regressor is constant throughout the history, since there is nothing to fit from it.\nExtra regressors are put in the linear component of the model, so the underlying model is that the time series depends on the extra regressor as either an additive or multiplicative factor (see the next section for multiplicativity).\n\n","metadata":{"_uuid":"e8f5961a2910431906a6afd73a72629815167997"}},{"cell_type":"markdown","source":"<a id=\"section5.7\"></a>\n# 5.7  Including Additional Regressors e.g. NFL Sundays","metadata":{"_uuid":"500eb60cb0158034b9833354b06a0224880d4e23","trusted":true}},{"cell_type":"code","source":"nfl_sunday <- function(ds) {\n  dates <- as.Date(ds)\n  month <- as.numeric(format(dates, '%m'))\n  as.numeric((weekdays(dates) == \"Sunday\") & (month > 8 | month < 2))\n}\nstats$nfl_sunday <- nfl_sunday(stats$ds)\n\n\nmodel_prophet <- prophet()\nmodel_prophet <- add_regressor(model_prophet, 'nfl_sunday')\nmodel_prophet <- add_seasonality(model_prophet, name='daily', period=60, fourier.order=5)\nmodel_prophet <- prophet(stats, holidays = holidays,holidays.prior.scale = 0.5, yearly.seasonality = 4,\n                         interval.width = 0.95,changepoint.prior.scale = 0.006,daily.seasonality = T)\nfuture = make_future_dataframe(model_prophet, periods = 90, freq = 'days')\nforecast = predict(model_prophet, future)","metadata":{"_uuid":"f6cb17360aa624eeda3f74d7b630d846da24723c","execution":{"iopub.status.busy":"2022-07-09T18:58:11.888069Z","iopub.execute_input":"2022-07-09T18:58:11.890972Z","iopub.status.idle":"2022-07-09T18:58:25.300132Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot(model_prophet, forecast) + add_changepoints_to_plot(model_prophet)\nprophet_plot_components(model_prophet, forecast)","metadata":{"_uuid":"168f5d8398519048baa57ab2a92b445f68a89858","execution":{"iopub.status.busy":"2022-07-09T18:58:39.382594Z","iopub.execute_input":"2022-07-09T18:58:39.386923Z","iopub.status.idle":"2022-07-09T18:58:42.930407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we can see the model has improved a lot after optimising the prophet parameters and including holidays, seasonality and additional regressors in the data. \nThe  trend is NOT much fluctuating like the baseline models and there is NO CHANGE POINTS of sales as well after fitting a better model. \nThe model is NOT much overfitting as well. \n  We can colnculde that Holidays has an effect on Sales Price and we have taken care of it in our optimised models.","metadata":{"_uuid":"133b456a4e14b8a11f2ba63e365f0f24e450cee2"}},{"cell_type":"code","source":"predict_store1_item1=data.frame(date=forecast$ds,forecast=expm1(forecast$yhat))\npredict_store1_item1$yearmonth=as.yearmon(predict_store1_item1$date)\n\ncolnames(predict_store1_item1)<-c(\"ds\",\"forecast\",\"yearmonth\")","metadata":{"_uuid":"5f7423067bcb4b9e325a7338ced9c9fb423b26d2","execution":{"iopub.status.busy":"2022-07-09T18:58:57.128849Z","iopub.execute_input":"2022-07-09T18:58:57.131086Z","iopub.status.idle":"2022-07-09T18:58:57.152220Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"section5.8\"></a>\n# 5.8.SMAPE CALCULATION\n\nUsed to estimate sMAPE","metadata":{"_uuid":"541ca84c05cd625c2f5d1fe94962c7b3afb094aa"}},{"cell_type":"markdown","source":"**Symmetric Mean Absolute Percent Error (SMAPE)** is an alternative to Mean Absolute Percent Error (MAPE) when there are zero or near-zero demand for items. SMAPE self-limits to an error rate of 200%, reducing the influence of these low volume items. Low volume items are problematic because they could otherwise have infinitely high error rates that skew the overall error rate.\n\nSMAPE is the forecast minus actuals divided by the sum of forecasts and actuals as expressed in this formula:","metadata":{"_uuid":"67b1a6a809184eca89f81c5fe5610b111de62e5c"}},{"cell_type":"markdown","source":"<img src=\"http://www.vanguardsw.com/vg-content/uploads/2017/01/Screen-Shot-2017-01-10-at-3.56.47-PM-copy.png\" width=20%>","metadata":{"_uuid":"c07e7288860f6ff9cc8db9ea6d51c5f5c6689dea"}},{"cell_type":"code","source":"smape_cal <- function(outsample, forecasts){\n  outsample <- as.numeric(outsample)\n  forecasts<-as.numeric(forecasts)\n  smape <- (abs(outsample-forecasts))/((abs(outsample)+abs(forecasts))/2)\n  return(smape)\n}\n\nstr(stats)\nstats$ds=as.Date(stats$ds)\npredict_store1_item1$ds=as.Date(predict_store1_item1$ds)\n\ntrain_predict=merge(stats,predict_store1_item1,by=\"ds\",all.x=T)\nSMAPE_ERR <- smape_cal(outsample=train_predict$y, forecasts=train_predict$forecast)\nSMAPE<-mean(SMAPE_ERR,na.rm = T)\nsprintf(\"The value of SMAPE for Store-1 & Item-1 is %f \", SMAPE )","metadata":{"_uuid":"3e922de4f6a470aa27d31f5f13b51e2e7125d313","execution":{"iopub.status.busy":"2022-07-09T18:59:08.588638Z","iopub.execute_input":"2022-07-09T18:59:08.590895Z","iopub.status.idle":"2022-07-09T18:59:08.642128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The value of SMAPE for Store-1 & Item-1 is 1.448148 '","metadata":{"_uuid":"b71f940ff58bdbd623fc6f9e66fb2623a6aad411"}},{"cell_type":"markdown","source":"<a id=\"section5.9\"></a>\n# 5.9. Automation for Prophet: Splitting data by Store and Item","metadata":{"_uuid":"2371ac34222e4228d13f2fa1f5d82a782db80292"}},{"cell_type":"code","source":"train$Year=NULL\ntrain$Month=NULL\nhead(train)\ntrain$sales=log1p(train$sales)\n\ncolnames(train)<- c(\"ds\",\"store\",\"item\",\"y\")\ntrain_splitting= split(train, by=c('store', 'item'), keep.by=FALSE)\nclass(train_splitting)\n\nprediction<-function(df)\n{\n  playoffs <- data_frame(\n    holiday = 'playoff',\n    ds = as.Date(c('2013-07-12', '2014-07-12', '2014-07-19',\n                   '2014-07-02', '2015-07-11', '2016-07-17',\n                   '2016-07-24', '2016-07-07','2016-07-24')),\n    lower_window = 0,\n    upper_window = 1\n  )\n  \n  #######  I have inlcuded the Holiday Sales fofr Festive seasons like New Year & Christmas for different Years in Superbowls.\n  superbowls <- data_frame(\n    holiday = 'superbowl',\n    ds = as.Date(c('2013-01-01', '2013-12-25', '2014-01-01', '2014-12-25','2015-01-01', '2015-12-25','2016-01-01', '2016-12-25',\n                   '2017-01-01', '2017-12-25')),\n    lower_window = 0,\n    upper_window = 1\n  )\n  holidays <- bind_rows(playoffs, superbowls)\n  \n  \n  \n  \n  model_prophet <- prophet()\n  model_prophet <- add_seasonality(model_prophet, name='daily', period=60, fourier.order=5)\n  model_prophet <- prophet(df, holidays = holidays,holidays.prior.scale = 0.5, yearly.seasonality = 4,\n                           interval.width = 0.95,changepoint.prior.scale = 0.006,daily.seasonality = T)\n  \n  \n  future = make_future_dataframe(model_prophet, periods = 90)\n  forecast = predict(model_prophet, future)\n  forecast_final<-  xts::last(forecast[, c(\"ds\",\"yhat\")],90)\nreturn(forecast_final)\n\n}\n\n\n","metadata":{"_uuid":"dcab08f467f58fbcfd391c662c6b18647d41f525","execution":{"iopub.status.busy":"2022-07-09T18:59:14.745897Z","iopub.execute_input":"2022-07-09T18:59:14.749228Z","iopub.status.idle":"2022-07-09T18:59:14.950639Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prediction_final=as.data.frame(sapply(train_splitting[c(1,2)],prediction))","metadata":{"_uuid":"0def529fd9317cd8084585bb38eb0b6ee0d5fb88","execution":{"iopub.status.busy":"2022-07-09T18:59:28.562634Z","iopub.execute_input":"2022-07-09T18:59:28.565332Z","iopub.status.idle":"2022-07-09T18:59:54.024699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"library(reshape)\ndim(prediction_final)\nmd <- melt(prediction_final)\ndim(md)\ncolnames(md)<-c(\"store\",\"date\",\"sales\")\nmd$sales=expm1(md$sales)\nhead(md)","metadata":{"_uuid":"075b02595dad62830f74f124d9d0e253e4568204","execution":{"iopub.status.busy":"2022-07-09T19:00:22.556202Z","iopub.execute_input":"2022-07-09T19:00:22.558493Z","iopub.status.idle":"2022-07-09T19:00:22.603645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"section5.10\"></a>\n# 5.10 Final prediction using Prophet","metadata":{"_uuid":"9e68acf3151afe737ec7afa2313fa0ab2498b794"}},{"cell_type":"markdown","source":"<img src=\"https://github.com/uselessskills/auto_ml/blob/master/Prophet_analysis/pngs/final_product.png?raw=true\" width=\"100%\">","metadata":{"_uuid":"7ee6c8969be6cc0c8b8ad95091e05d06968db751"}},{"cell_type":"markdown","source":".**NOTE: HOW TO GET & SUBMIT THE OUTPUT: **\n\n\n\nThere are 10 STORES and 50 ITEMS. Therefore there are 500 elements in \"train_splitting\" list. I have ran it for 2 elements.  Ran the whole list to get the full prediction and then ran the below code. You will get the final output for submission. ","metadata":{"_uuid":"9050fd013ba136c36ce798e04df1b8ebd0326235"}},{"cell_type":"code","source":"#sub=fread(\"../input/sample_submission.csv\")\n#sub$sales=md$sales\n#head(sub)\n#fwrite(sub,\"sub_prophet_v1.csv\")","metadata":{"_uuid":"6a39155a26050d62ee1eb51f1ea2329f37ea39c6","execution":{"iopub.status.busy":"2022-07-09T19:00:38.882549Z","iopub.execute_input":"2022-07-09T19:00:38.885304Z","iopub.status.idle":"2022-07-09T19:00:38.905203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"section5.11\"></a>\n# 5.11.HOW TO IMPROVE PROPHET MODELS FURTHER #","metadata":{"_uuid":"d20623d9be53e8b407aded6c713bbd9319f7b02c"}},{"cell_type":"markdown","source":"1. More Optimising the parameters using Bayesian Optimisation.\n2. Include Multiplicative Seasonality in the model.\n3. Remove the outliers and adjust Trend Flexibility.\n4. By defualt Prophet uses LINEAR GROWTH for training the model. But we can use the LOGISTIC GROWTH in case of Multiplicative Seasonality.\n\n","metadata":{"_uuid":"6bbcf7e9d7fedee0a0fe75e1f96f6dc9d8f4f20e"}},{"cell_type":"markdown","source":"<a id=\"section6\"></a>\n# 6. PART 2: Forecasting Time Series with Neural Networks: ","metadata":{"_uuid":"e9e1b702832a3729bb7a2578f8e82893a7ba5d26"}},{"cell_type":"markdown","source":"<a id=\"section6.1\"></a>\n# 6.1. Background","metadata":{"_uuid":"46d6cd8ed80ac05418e909e45a869c0fae871e2a"}},{"cell_type":"markdown","source":"Artificial neural network (ANN) is a widely used pattern-recognition methodology for machine learning. ANN is an  emulation of biological neural network, which is composed of many interconnected neurons. However, it only utilized a very limited set of concepts from its biological counterpart. An ANN could have one or more layer of neurons. They could be fully or partially connected. Each connection between two nodes has a weight, which encapsulate the “knowledge” of the system. By processing existing cases\nwith inputs and expected outputs, these weights would be adjusted based on differences between actual and expected outputs. Because of the nonlinear fashion of ANN, they could be used in a lot of business applications. ","metadata":{"_uuid":"11e0f2b53fb9eccc2873ea7b51884c696b414621"}},{"cell_type":"markdown","source":"<img src=\"https://i.imgur.com/0tpEWj3.png\" width=100%>","metadata":{"_uuid":"481bf8a22a436fe818c9c04a5c324f16d1e9dcc9"}},{"cell_type":"markdown","source":"<img src=\"https://i.imgur.com/VKSe7EG.png\" width=100%>","metadata":{"_uuid":"82c8523d26f50ff526f1556e5a2bb0485650358b"}},{"cell_type":"markdown","source":"<img src=\"https://i.imgur.com/9xti9Ao.png\" width=100%>","metadata":{"_uuid":"3f7808e51643abf143e94b937c48f7fd5fe83f34"}},{"cell_type":"markdown","source":"<img src=\"https://i.imgur.com/53Xf7tb.png\" width=100%>","metadata":{"_uuid":"367e8744933c7e7236db47ca8a5a99ce4b8bbd63"}},{"cell_type":"markdown","source":"<img src=\"https://i.imgur.com/VevpkhG.png\" width=100%>","metadata":{"_uuid":"207fbcbc895649fae7bf11734acd172bdfbbec5e"}},{"cell_type":"markdown","source":"<img src=\"https://i.imgur.com/9pHnxl1.png\" width=100%>","metadata":{"_uuid":"f0ba7e6aab5fc8be1b4d292b7c372a2d418ce9b3"}},{"cell_type":"markdown","source":"<img src=\"https://i.imgur.com/GnpPKou.png\" width=100%>","metadata":{"_uuid":"6f778024b524b271cf064ad8af2fd76b8ae2c1bc"}},{"cell_type":"markdown","source":"# We are going to use nn for package developed by Nikolaos Kourentzes. \n\n**NOTE: Added the theoretical background in the Comments section.****","metadata":{"_uuid":"9d2b624229ace14d87eade19624275d611cdcfc2"}},{"cell_type":"markdown","source":" Currently there are two types of neural network available, both **feed-forward:**\n\na. **multilayer perceptrons** (use function **MLP**); \n\n b. **extreme learning machines** (use function elm)","metadata":{"_uuid":"9424fb04dbc559942f329fea9e58bdd4224377ef"}},{"cell_type":"code","source":"y<-ts(stats$y,frequency=365,start = 2013,end=2017)\nhead(y)\n\nplot(y)","metadata":{"_uuid":"bd91ddd6891e59b9a0e027a84d53be09a0933976","execution":{"iopub.status.busy":"2022-07-09T19:04:10.248495Z","iopub.execute_input":"2022-07-09T19:04:10.251033Z","iopub.status.idle":"2022-07-09T19:04:10.447239Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"section6.2\"></a>\n# 6.2.Baseline Predcition with 1 Hidden Layer \n\nAlthough I have used **MLP(Multi Layer Perceptron) function, for basline prediction I have used a single hidden layer, No differencing, No lags**, to make a simpler model.\n\n We will add more layers of complicated Network with **Diffrencing order and lags in ELM (Extreme Learning Machines) **","metadata":{"_uuid":"a92c22bc9155a491f53c66ace079e0c214684bea"}},{"cell_type":"code","source":"h <- 90   # We will predict for 90 days sales starting from 01-JAN-2018.\ntt <- cbind(c(1:(length(y)+h),rep(0,2*h)))\n# Observe that the deterministic trend ends with zeros\n#print(tt)\n\n# Fit a network with no differencing, no univariate lags, and fixed deterministic trend\nfit1 <- mlp(y,difforder=0,lags=0,xreg=tt,xreg.lags=list(0),xreg.keep=TRUE)\n","metadata":{"_uuid":"40ef72377130a430f41cb443b58cba12c0ea916b","execution":{"iopub.status.busy":"2022-07-09T19:04:30.510260Z","iopub.execute_input":"2022-07-09T19:04:30.512805Z","iopub.status.idle":"2022-07-09T19:08:14.673932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(fit1)\nplot(fit1)","metadata":{"_uuid":"44e41bf19c87b482ff3cc0b91dbde7dea337d45b","execution":{"iopub.status.busy":"2022-07-09T19:09:34.152345Z","iopub.execute_input":"2022-07-09T19:09:34.154586Z","iopub.status.idle":"2022-07-09T19:09:34.311662Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Therefore we can observe our MLP with single hidden layer. Let's see how the fitted graph looks on our training data of Store -1 & Item-1\n\nThis is the basic command to fit an MLP network to a time series. This will attempt to automatically specify autoregressive inputs and any necessary pre-processing of the time series. With the pre-specified arguments it trains 20 networks which are used to produce an ensemble forecast and a single hidden layer with 5 nodes. You can override any of these settings. The output of print is a summary of the fitted network:","metadata":{"_uuid":"6cc60ac1f2310af195f3da6a3a60ab8d21170f7d"}},{"cell_type":"markdown","source":"# The light red inputs represent the binary dummies used to code seasonality, while the grey ones are autoregressive lags. To produce forecasts you can type:","metadata":{"_uuid":"f787e8fff8fc67270eb59810e24d123d12987b58"}},{"cell_type":"markdown","source":"<a id=\"section6.3\"></a>\n#  6.3.Plotting the Forecast & Calculate the MSE (Mean Squared Error)","metadata":{"_uuid":"1f016ef7a961504f65d1f72da3f58bee612a4bc4","trusted":true}},{"cell_type":"code","source":"plot(forecast(fit1,h=h,xreg=tt))\n\nprint(\"The MSE for Store-1 & Item -1 is\")\nprint(round(fit1$MSE,4))","metadata":{"_uuid":"55fa696443d5c4b384f94d7dd28065c6874eac61","execution":{"iopub.status.busy":"2022-07-09T19:10:32.727803Z","iopub.execute_input":"2022-07-09T19:10:32.730188Z","iopub.status.idle":"2022-07-09T19:10:33.292936Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**MSE  for Store-1 & Item- 1: is 0.0674 (Not bad as of the baseline right)**","metadata":{"_uuid":"7bcdf42879d9d92d7b0e52f1d4813de44f837ee9"}},{"cell_type":"markdown","source":"<a id=\"section6.4\"></a>\n# 6.4. Shifting Input","metadata":{"_uuid":"db5e2e78642e9ee3d242acb457a05bcc036ff170"}},{"cell_type":"markdown","source":"# Now let us shift the input so that the zeros are in the forecast period","metadata":{"_uuid":"4295dbc7af434d2cf08efca1cf6a1485904e98e8"}},{"cell_type":"code","source":"tt2 <- tt[-(1:h),,drop=FALSE]\nplot(forecast(fit1,h=h,xreg=tt2))\n","metadata":{"_uuid":"961463410582cb225b821bd7be9ac67cd49c30ae","execution":{"iopub.status.busy":"2022-07-09T19:10:55.939976Z","iopub.execute_input":"2022-07-09T19:10:55.942218Z","iopub.status.idle":"2022-07-09T19:10:56.604231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" The seasonality is there, **but there is zero trend**, as the inputs suggest. \n\nAlso note that the **mlp modelled multiplicative seasonality on its own**. NNs are cool. ","metadata":{"_uuid":"e79c41b29f7ee5d14b429de0dcd52d07aa3d56da"}},{"cell_type":"markdown","source":"","metadata":{"_uuid":"2f7ae94c6b6ab0bdcedd9e8accb9253529c20996"}},{"cell_type":"markdown","source":" <a id=\"section6.5\"></a>\n # 6.5. Outplot & Model Fitting:\n  Now let us fit a network on the shifted inputs  I will ask for outplot=1 to see the model fit\n","metadata":{"_uuid":"cd680a37b198a4ab381caedae41707e11548d3cf"}},{"cell_type":"code","source":"fit2 <- mlp(y,difforder=0,lags=0,xreg=tt2,xreg.lags=list(0),xreg.keep=TRUE,outplot=1)","metadata":{"_uuid":"01c460ea00000a44e6c5758f098e4b8136b86620","execution":{"iopub.status.busy":"2022-07-09T19:11:12.715478Z","iopub.execute_input":"2022-07-09T19:11:12.716893Z","iopub.status.idle":"2022-07-09T19:16:11.794914Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"** This fitted model has taken the Multiplicative seasonality into account.**","metadata":{"_uuid":"09cf564161e502e7238d21c4ed2e851b34b9f50a"}},{"cell_type":"code","source":"print(fit2)\nplot(forecast(fit2,h=h,xreg=tt2))","metadata":{"_uuid":"b8fd5560642ddfc1cbc6e856e6207f8a5245ef36","execution":{"iopub.status.busy":"2022-07-09T19:16:26.587499Z","iopub.execute_input":"2022-07-09T19:16:26.589649Z","iopub.status.idle":"2022-07-09T19:16:27.106667Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Let's check the MSE.\nprint(round(fit1$MSE,4))  #0.0671 ","metadata":{"_uuid":"d4c7b50577a7b6a0eab1278903ed4ea4a10ceb7f","execution":{"iopub.status.busy":"2022-07-09T19:16:34.055595Z","iopub.execute_input":"2022-07-09T19:16:34.057763Z","iopub.status.idle":"2022-07-09T19:16:34.072516Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"You can also let it **choose the number of hidden nodes**. There are various options for that, but all are **computationally expensive**.","metadata":{"_uuid":"606dee90d82f5132b8501f0a79c4db7245cfbe3a"}},{"cell_type":"markdown","source":"<a id=\"section7\"></a>\n# 7. Extensions:\n\n**1. Choosing Oprimal No of Hidden Layer based on min(MSE).**\n\n**2. Fitting ELM Feed Forward Network for forecasting.**\n\n**3.  Temporal Hierarchies forecasting**\n\n**4. Automation based on MLP/ ELM NN Forecasting.**","metadata":{"_uuid":"33da475463f3601bfe7d9de5a0636b4ee197c563","trusted":true}},{"cell_type":"markdown","source":"<a id=\"section8\"></a>\n# 8. References:\n\n-----------------\n**PART 1 : PROPHET**\n\na.[Facebook Prophet Github](https://github.com/facebook/prophet)\n\nb. [Favebook Prophet Blog](http://research.fb.com/prophet-forecasting-at-scale/)\n\nc. [Prediction using Prophet in Python by Conor](https://www.kaggle.com/conorm97/prediction-using-prophet-with-us-holidays)\n\nd. [Prophet Analytics Vidhya](https://www.analyticsvidhya.com/blog/2018/05/generate-accurate-forecasts-facebook-prophet-python-r/)\n\n---------------------\n**PART 2: Forecasting uisng ANN:**\n\na.[Neural network Forecasting: Theory](http://www.neural-forecasting.com/Downloads/EVIC05_tutorial/EVIC%2705%20Slides%20-%20Forecasting%20with%20Neural%20Networks%20Tutorial%20SFCrone.pdf)\n\nb. [NNFOR Github](https://github.com/trnnick/nnfor)\n\nc. [introduction to Neural Network Forecasting](http://www.neural-forecasting.com/)","metadata":{"_uuid":"8149d4a112e9af88fbd46ba2f37b5ef251efdb45"}},{"cell_type":"markdown","source":"Thank you for keeping your patience and come till the end.\n\nHappy Kaggling :) ","metadata":{"_uuid":"b096e68dad67b72512b5244609570e75817c44b5","trusted":true}}]}