--- title: "Analysis of aggregated item series" author: "Andriy Levitskyy" date: '`r Sys.Date()`' output: html_document: number_sections: true fig_caption: true toc: true fig_width: 7 fig_height: 4.5 theme: cosmo highlight: tango code_folding: hide --- ```{r setup, include=FALSE, echo=FALSE} knitr::opts_chunk$set(echo=TRUE, error=FALSE) ``` # Introduction Being my first take in EDA, here I try to visualise the dataset by plotting moving avaerages to deal with weekly, monthly (I guess payday effect) and yearly (holiday effects) seasonality. I use 2 years worth of data, because this is the minimal time period one can actually observe yearly seasonal patterns (for example, ice cream is consumed less in winter). Based on this moving averages, I caluclate sample statistics and apply hierarchical clustering analysis on them, to see if it can spot any similarities in behaviour of the time series. Finally, I check the cannibalisation hypothesis (items in the same class compete with each other) by using correlation analysis. Another good test for cannibalisation hypothesis would be running regression on own onpromotion columns and another item onpromotion series, but it is not implemented here. # Preparations {.tabset .tabset-fade .tabset-pills} ## Loading libraries ```{r, message=FALSE, results='hide'} library(tibble) library(data.table) library(lubridate) library(dplyr) library(R.utils) library(ggplot2) library(gplots) library(foreach) library(doMC) library(forecast) library(Kendall) library(kableExtra) ``` ## Loading dataset ```{r, message=FALSE, results='hide'} train <- as.tibble(fread('../input/train.csv',nrows = 3)) true_names<-colnames(train) l2keep <- 70e6 nL <- countLines("../input/train.csv") train <- fread('../input/train.csv',header=FALSE,nrows = l2keep,skip=nL-l2keep) colnames(train)<-true_names rm(l2keep) rm(nL) rm(true_names) #oil = fread('input/oil.csv') test <- fread('../input/test.csv') #transactions <- fread('input/transactions.csv') stores <- fread('../input/stores.csv') items <- fread('../input/items.csv') ``` ## Reformatting data ```{r, message=FALSE, results='hide'} test <- test %>% mutate(date = ymd(date), store_nbr = as.factor(store_nbr), item_nbr = as.factor(item_nbr)) train <- train %>% mutate(date = ymd(date), store_nbr = as.factor(store_nbr), item_nbr = as.factor(item_nbr)) test$dow<-weekdays(test$date) train$dow<-weekdays(train$date) train %>% group_by(item_nbr,date) %>% summarize(unit_sum=sum(unit_sales),sumpromo=sum(onpromotion))%>% arrange(item_nbr,date)->item_series ``` ## Calculating moving averages ```{r, message=FALSE, results='hide'} timseries<-dcast(item_series, item_nbr ~ date, value.var = "unit_sum", fill = 0) true_names<-timseries[,1] timseries<-timseries[,c(-1)] ts_length<-length(timseries[1,]) registerDoMC(7) start_time <- Sys.time() #don`t forget to change here length(rownames(train3)) ma_weekly <- foreach(i=1:nrow(timseries),.combine=rbind, .packages=c("pracma")) %dopar% { ma_weekly <- movavg(as.numeric(timseries[i,]),n=7,type="s")[8:ts_length] } print("ma weekly done") ma_monthly <- foreach(i=1:nrow(timseries),.combine=rbind, .packages=c("pracma")) %dopar% { ma_monthly <- movavg(as.numeric(timseries[i,]),n=30,type="s")[31:ts_length] } print("ma monthly done") ma_yearly <- foreach(i=1:nrow(timseries),.combine=rbind, .packages=c("pracma")) %dopar% { ma_yearly <- movavg(as.numeric(timseries[i,]),n=365,type="s")[366:ts_length] } print("ma yearly done") print("the run time") print(Sys.time()-start_time) ``` ## Plotting function ```{r, message=FALSE, results='hide'} ma_plot<-function(item_nubmer){ par(mfrow=c(2,2)) plot(as.numeric(ma_yearly[item_nubmer,]),type="l",ylab="yearly MA") plot(as.numeric(ma_monthly[item_nubmer,]),type="l",ylab="monthly MA") plot(as.numeric(ma_weekly[item_nubmer,]),type="l", ylab="weekly MA", main=paste0(item_nubmer,"/", true_names[item_nubmer], "/", subset(items,item_nbr==true_names[item_nubmer])$family, "/", subset(items,item_nbr==true_names[item_nubmer])$class)) plot(as.numeric(timseries[item_nubmer,]),type="l", ylab="daily plot") } ``` ## Generating statistics ```{r, message=FALSE, results='hide'} stats <- foreach(i=1:nrow(ma_weekly),.combine=rbind, .packages=c("moments")) %dopar% { stats <- c(mean(as.numeric(timseries[i,])), sd(as.numeric(timseries[i,])), skewness(as.numeric(timseries[i,])), kurtosis(as.numeric(timseries[i,])), MannKendall(timseries[i,])$tau[1], mean(as.numeric(ma_weekly[i,])), sd(as.numeric(ma_weekly[i,])), skewness(as.numeric(ma_weekly[i,])), kurtosis(as.numeric(ma_weekly[i,])), MannKendall(ma_weekly[i,])$tau[1], mean(as.numeric(ma_monthly[i,])), sd(as.numeric(ma_monthly[i,])), skewness(as.numeric(ma_monthly[i,])), kurtosis(as.numeric(ma_monthly[i,])), MannKendall(ma_monthly[i,])$tau[1], mean(as.numeric(ma_yearly[i,])), sd(as.numeric(ma_yearly[i,])), skewness(as.numeric(ma_yearly[i,])), kurtosis(as.numeric(ma_yearly[i,])), MannKendall(ma_yearly[i,])$tau[1], sum(timseries[i,]==0)) } stats<-as.data.frame(stats) colnames(stats)<-c("mean_daily","sd_daily","skewness_daily","kurtosis_daily","trend_daily", "mean_weekly","sd_weekly","skewness_weekly","kurtosis_weekly","trend_weekly", "mean_monthly","sd_monthly","skewness_monthly","kurtosis_monthly","trend_monthly", "mean_yearly","sd_yearly","skewness_yearly","kurtosis_yearly","trend_yearly", "num_zeros") rownames(stats)<-1:4036 stats$sd_mean_ratio_daily<-stats$mean_daily/stats$sd_daily stats$sd_mean_ratio_weekly<-stats$mean_weekly/stats$sd_weekly stats$sd_mean_ratio_monthly<-stats$mean_monthly/stats$sd_monthly stats$sd_mean_ratio_yearly<-stats$mean_yearly/stats$sd_yearly stats<-subset(stats,select=-c(mean_daily,mean_weekly,mean_monthly,mean_yearly)) ``` # Plots of the first 10 items {.tabset .tabset-fade .tabset-pills} In this section, one can have a look at the moving averages plots. Some of the trends may not be apparent by looking at the daily data due to weekly seasonality which makes data look noisy, however monthly and yearly moving averages provide smoother representation of the trend. The title of the plots are of the following format "position in my data frame/ item number/ family/ class" ## Item 1 This is an example of an item which was just recently introduced and only has 200 days worth of transcations history. I do not think that one can rely on predictive time series modelling for this particular item. ```{r, message=FALSE} ma_plot(1) ``` ## Item 2 This is an example of a seasonal grocery good, which is only popular during the summer season and which is most likely to be on decling during the test period ```{r, message=FALSE} ma_plot(2) ``` ## Item 3 This is a cleaning product which has quite stable demand over the 2 year period. One can clearly see the monthly seasonal patters by looking atweekly moving average plot. ```{r, message=FALSE} ma_plot(3) ``` ## Item 4 Another good with stable demand. An interesting feature of time series behaviour of this good is a strong holiday effect observed near observation 100 and 500. (THe distance between the spikes should be around 365 days) ```{r, message=FALSE} ma_plot(4) ``` ## Item 5 Looking at the yearly moving average, it is possible to see a clear downward trend which is not that apparent when looking at the daily plot. ```{r, message=FALSE} ma_plot(5) ``` ## Item 6 The only interesting thing about this grocery item is a massive spike around obsevation 200. ```{r, message=FALSE} ma_plot(6) ``` ## Item 7 Strong monthly seasonality with some time of the year spikes. One of the spikes happened just now and may not persist up until the end of the test period. ```{r, message=FALSE} ma_plot(7) ``` ## Item 8 There is a holiday effects around 100 and 500 observations obeseved before, monthly sesonality and time of the year demand surges. ```{r, message=FALSE} ma_plot(8) ``` ## Item 9 The behaviour is quite similar to item 6 and it is not surprising because this items are of the same class. ```{r, message=FALSE} ma_plot(9) ``` ## Item 10 ```{r, message=FALSE} ma_plot(10) ``` #Some plots based on statistics {.tabset .tabset-fade .tabset-pills} In 2.6 there is ascript which calculates first 4 central moments of each of the 3 moving averages and time series itself. I also calculate Mann-Kendal statistic which measure if the series are trending for each of the series. Finally, I calculate number of zero observations in the item series. Given this infomration I was curios to see the behaviour of the items with the largest/smallest calculated statistics, and this section presents my findings. In the analysis I exclude the recently discontinues or launched items as defined in the next sub section. ## Discontinued and recently introduced items The code below finds the discontinued items (the criterion for this is that last 15 observations are 0. My guess is that for these observations it is safe to assign 0s in the test set), and items which were launched in the last 2 years (i.e. if the first 100 observations are 0). The final table outputs items with the largest number of zero observations in the increasing order. For the items in the bottom of the table, predictive analysis will not do well at all. ```{r, message=FALSE} library(knitr) row_sum_last_15<-apply(timseries[,(ts_length-15):ts_length],1,sum) row_sum_first_100<-apply(timseries[,1:100],1,sum) items_discontinued<-rownames(stats)[row_sum_last_15==0] items_discontinued<-as.numeric(items_discontinued) print("Discontinued items") print(true_names[items_discontinued]) items_recently_started<-rownames(stats)[row_sum_first_100==0] items_recently_started<-as.numeric(items_recently_started) number_of_zeros<-c() for(i in items_recently_started){ number_of_zeros<-c(number_of_zeros,sum(timseries[i,]==0)) } number_of_zeros<-as.data.frame(cbind(true_names[items_recently_started],number_of_zeros)) colnames(number_of_zeros)<-c("item_nbr","num_zeros") number_of_zeros %>% arrange(num_zeros) -> number_of_zeros items_with_no_history<-number_of_zeros$item_nbr[number_of_zeros$num_zeros>630] stats_no_discontinued<-stats[-c(items_discontinued,items_with_no_history),] stats_for_clustering<-stats[-c(items_discontinued,items_recently_started),] kable(number_of_zeros, "html") %>% kable_styling() %>% scroll_box(width = "500px", height = "200px") ``` ## Strongest downward trend Based on minimal values of the Mann-Kendal statistic. This beverage lost a lot of its popularity over the time ```{r, message=FALSE} ma_plot(as.numeric(rownames(stats_for_clustering)[stats_for_clustering$trend_yearly== min(stats_for_clustering$trend_yearly)])) ``` ## Strongest upward trend Based on maximal values of the Mann-Kendal statistic. The popularity of this produce grown a lot over the last year, however was recently on decline. I would be very careful predicting this particular item. ```{r, message=FALSE} for(item_nbr in as.numeric(rownames(stats_for_clustering)[stats_for_clustering$trend_yearly== max(stats_for_clustering$trend_yearly)])){ ma_plot(item_nbr) } ``` ## The least stable demand Based on the largest standard deviation/mean ratio. This beverage shows some interesting cyclical behaviour (which is probably induced by promotions, but can not say for sure) ```{r, message=FALSE} ma_plot(as.numeric(rownames(stats_for_clustering)[stats_for_clustering$sd_mean_ratio_yearly== max(stats_for_clustering$sd_mean_ratio_yearly)])) ``` ## The most stable demand Based on the smalles standard deviation/mean ratio. I guess this definition of stability is not very good. This item was not classified as recently started due to very small purchases in the start of the sampling period but most of the observations are 0. I guess the statistic shows that this item has very stable 0 demand. ```{r, message=FALSE} ma_plot(as.numeric(rownames(stats_for_clustering)[stats_for_clustering$sd_mean_ratio_yearly== min(stats_for_clustering$sd_mean_ratio_yearly)])) ``` ## Highest monthly MA kurtosis One need to be very careful with this item, because of the recent jump. The test period behaviour is highly unpredictable. ```{r, message=FALSE} ma_plot(as.numeric(rownames(stats_for_clustering)[stats_for_clustering$kurtosis_monthly== max(stats_for_clustering$kurtosis_monthly)])) ``` ## Lowest monthly MA kurtosis This item was trialed around observation 100 and then only started being introduced around observation 300 ```{r, message=FALSE} ma_plot(as.numeric(rownames(stats_for_clustering)[stats_for_clustering$kurtosis_monthly== min(stats_for_clustering$kurtosis_monthly)])) ``` ## Highest yearly MA kurtosis There were 2 big spikes, but otherwise no demand. It could be one of the items which are only really purchased during holidays. ```{r, message=FALSE} ma_plot(as.numeric(rownames(stats_for_clustering)[stats_for_clustering$kurtosis_yearly== max(stats_for_clustering$kurtosis_yearly)])) ``` ## Lowest yearly MA kurtosis 4 big spikes in demand observed. ```{r, message=FALSE} ma_plot(as.numeric(rownames(stats_for_clustering)[stats_for_clustering$kurtosis_yearly== min(stats_for_clustering$kurtosis_yearly)])) ``` ## Highest monthly MA skewness While the demand is historically stable, there was a very recent jump which may have long-term effect but may not. ```{r, message=FALSE} ma_plot(as.numeric(rownames(stats_for_clustering)[stats_for_clustering$skewness_monthly== max(stats_for_clustering$skewness_monthly)])) ``` ## Lowest monthly MA skewness Very stable demand following the launch ```{r, message=FALSE} ma_plot(as.numeric(rownames(stats_for_clustering)[stats_for_clustering$skewness_monthly== min(stats_for_clustering$skewness_monthly)])) ``` ## Highest yearly MA skewness Moving averages based forecasts may suffer on this one. There was a very recent surge in demand which is not yet over. While I have no idea what is going to happen, I guess I would bet the current value to hold for next 3-4 days but then return to the historical mean. ```{r, message=FALSE} ma_plot(as.numeric(rownames(stats_for_clustering)[stats_for_clustering$skewness_yearly== max(stats_for_clustering$skewness_yearly)])) ``` ## Lowest yearly MA skewness ```{r, message=FALSE} ma_plot(as.numeric(rownames(stats_for_clustering)[stats_for_clustering$skewness_yearly== min(stats_for_clustering$skewness_yearly)])) ``` # Clustering on Statistics {.tabset .tabset-fade .tabset-pills} In this section I am applying hierarchical clustering on the sample statistics described and calcualted above. I am cutting the tree to get 20 biggest clusters. Then I am plotting 3 random items in 9 of the clusters trying to guess what criterion the algorithm uses to assign the items to each of the clusters. I found this task very hard to accomplish, it could be that the assignment to clusters is almost random. For example, clusters 2,4 and 9 are very hard to explain. Cluster 1 (based on the history of 3 items in it) reperesents items which has discontinuities in its demand (0 purchases during extended periods of time). Cluster 3 and 5 show items which experienced very large spike/spikes in demand with otherwise low demand relative to them. I think cluster 6 has some inherent form of seasonality but I am not sure. In cluster 8 all 3 plots have downward trends but can not say much more about it. ```{r, message=FALSE} clusters <- hclust(dist(stats_for_clustering)) clusterCut <- cutree(clusters, 20) print("Number of items in each cluster") print(table(clusterCut)) ``` ```{r, message=FALSE, results='asis'} stats_for_clustering$clust<-clusterCut for(i in 1:9){ cat(" \n##", "cluster",i," \n") random_sample<-sample(rownames(stats_for_clustering)[stats_for_clustering$clust==i],3) for(z in random_sample){ ma_plot(as.numeric(z)) cat(" \n") } } ``` # Canibalization {.tabset .tabset-fade .tabset-pills} Canibalization occurs when one item in a class takes market share of another item in the same class. One way canibalisation effects can be realized on data is through negative correlation in movement of two items in the same class. In this section I use correlation analysis to find such items. Canibalization however may occur even if the correlation is positive. For example, if both items have similar day of the week pattern the correlation may still be positive even if one item takes market share of another item. In this case promotions act as natural experiments. Regressing unit_sales on own onpromotion variable and the other item onpromotion item, and then checking for significance may provide a formal test of this hypothesis. In the next subsection, I check the classes which has more than 2 items in them. Then for each of these classes I take 3 random items and calculate all pairwise correlations and record the smallest value. The table of these values are reported in the next subsection. Other subsections plot item series for 3 of the classes with negative inter-class correlations. In the first two, items were discarded in favour of other items. In the last case with correlation -0.3 one item was slowly losing it popularity to another item. ## Table of interclass correlations ```{r, message=FALSE} class_numbers<-list() for(i in true_names){ class_numbers[length(class_numbers)+1]<-items$class[which(items$item_nbr==i)] } class_numbers<-unlist(class_numbers) class_numbers<-cbind(true_names,class_numbers) class_numbers<-as.data.frame(class_numbers) colnames(class_numbers)<-c("item_nbr","class_nbr") big_classes<-names(table(class_numbers$class_nbr))[table(class_numbers$class_nbr)>2] class_numbers2<-data.frame(item1=numeric(),item2=numeric(),item3=numeric()) for(i in big_classes){ class_numbers2[nrow(class_numbers2)+1,]<-sample(subset(class_numbers,class_nbr==i)$item_nbr,3) } minimum_pairwise_correlations<-function(class_number_row){ correlation1<-cor(x=t(timseries[as.integer(class_number_row[1]),]), y=t(timseries[as.integer(class_number_row[2]),])) correlation2<-cor(x=t(timseries[as.integer(class_number_row[1]),]), y=t(timseries[as.integer(class_number_row[3]),])) correlation3<-cor(x=t(timseries[as.integer(class_number_row[2]),]), y=t(timseries[as.integer(class_number_row[3]),])) return(min(correlation1,correlation2,correlation3)) } vector_of_minimum_correlations<-apply(class_numbers2,1,minimum_pairwise_correlations) class_numbers2<-cbind(class_numbers2,vector_of_minimum_correlations,big_classes) colnames(class_numbers2)<-c("item1","item2","item3","correlat","class_nbr") class_numbers2 %>% arrange(correlat) -> class_numbers2 kable(class_numbers2,"html") %>% kable_styling() %>% scroll_box(width = "500px", height = "200px") ``` ## Smallest correlation plot ```{r, message=FALSE} ma_plot(430) ma_plot(1198) ma_plot(3959) ``` ## Second small correlation plot ```{r, message=FALSE} ma_plot(606) ma_plot(238) ma_plot(3958) ``` ## Correlation -0.3 ```{r, message=FALSE} ma_plot(2895) ma_plot(2812) ma_plot(3047) ```