The topic I am studying is the revenue or sales for Indigo Books and Music, which is Canada’s largest book, gift, and specialty toy retailer. It is important to study this topic because there are larger uncertainty and a more rapid change in today’s business environment, and therefore there is a heavier role to play within predicting future sales (sales forecasting). Prediction becomes more important in order to not lose market shares and has become quite a common practice in the retail industry. Moreover, as competition grows stronger, companies continuously need to come up with new advantages to compete and survive. By forecasting sales, Indigo will be able to better strategize how it wants to operate in the future, in order to meet demands, perhaps realize that demand is decreasing and change needs to be employed. In order to conduct this study, the Bloomberg terminal was used to find quarterly revenue data for the years 2002 quarter one to 2017 quarter three.

To begin with, it is key to conduct statistical analysis and data visualization to get a good understanding of the properties of the data. To get a good insight it is important to determine whether the data has these three components of interest; a trend component, a seasonal component and an irregular component.

Load Library and Indigo sales data, indicate that it is a time series, and create a name variable
library(fpp)
## Loading required package: forecast
## Loading required package: fma
## Loading required package: expsmooth
## Loading required package: lmtest
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: tseries
library(knitr)
library(forecast)
library(vars)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following objects are masked from 'package:fma':
## 
##     cement, housing, petrol
## Loading required package: strucchange
## Loading required package: sandwich
## Loading required package: urca
library(nnet)
indigo_data <- read.csv("~/Dropbox (Personal)/ECON4210/Assignment Two/indigo_data.csv")
indigo_data <- ts(indigo_data, start = 2002, frequency = 4)
data = indigo_data[,"sales"]

Considering trend and seasonal influences are very important for the results of time series analysis. When observing the series against time, it is evident that Indigo sales fluctuate quite dramatically. There have been several periods where the sales peaked, reaching nearly 400 million, however, in a quarter or two they plummet down to just over 150 million. This could be because one or two quarters are more profitable than the other due to various factors, but overall indigo sales stay over the years stays nearly constant. This points towards seasonality. Moreover, there was an upward trend from 2002 to around 2010, but then it shifted down before slightly picking up again in around 2014. This proves that the data is quite random, and there is no evident trend.

Plot Autocorrelation
Acf(data, main="Indigo Sales Autocorrelation")

The autocorrelation of Indigo Sales graph shows strong evidence of autocorrelation in sales. The autocorrelation graph shows it peaking upward at when we analyze the residuals and the ACFs of the residuals, it is apparent that there are large spikes on lag 4, 8, 12 and 16 going upwards, and large spikes going downward at 2, 6, 10, and 14. This shows that there is a high degree of both pattern (which appears to be seasonal) and autocorrelation of the residuals. More important, indigo sales shows repeating patterns of positive and negative autocorrelations, typical for seasonality.

Therefore, it is safe to assume that indigo sales for this time period show a strong seasonal pattern but no trend. The season pattern is not confused with a cyclical behavior because the period is unchanging and associated with an aspect of the calendar (quarters).

Next we will conduct a time series decomposition which determine if there are irregular components in the data set, basically determine is there is anything else in the time series.

TSDfit <- stl(data,s.window="periodic", robust=TRUE)
plot(TSDfit)

Due to the fact that there is high seasonality, as evident in both the data and seasonal graph, it is safe to assume that the seasonal pattern is constant through time, and to reflect that the parameter has been set to be infinite, forcing the seasonal component to be periodic (s.window=”periodic”). The entire data set will be used to perform an analysis. However, we are more interested in the remainder component, which is what is left over when the seasonal and trend-cycle components have been subtracted from the data. It is apparent that near around 2012, sales increased in certain quarters, reaching a high at the end of 2017 (corresponding to some large positive values in the remainder component). This should be taken into consideration as it may indicate that a drift component may be necessary to accurately forecast the data.

Now that we have looked at the data set, it is time to make forecasts using different methods. We will begin first by setting up our training period and looking as the mean, naive and seasonal naive forecasts for indigo sales.

Set the training and testing windows
train <- window(data,start=c(2002, 1),end=c(2014, 2))
test <- window(data, start=c(2014,2),end=c(2017, 3))
both <- window(data,start=c(2002, 1))
h=length(test)
Forecast using Average, Naïve, and Seasonal Naïve Method
Indigofit1 <- meanf(train, h=h)
Indigofit2 <- naive(train, h=h)
Indigofit3 <- snaive(train, h=h)
Plot forecasts for Average, Naïve, and Seasonal Naïve Method
plot(Indigofit1, PI=FALSE,
     main="Forecasts for quarterly Indigo sales")
lines(Indigofit2$mean,col=2)
lines(Indigofit3$mean,col=3)
legend("topleft",lty=1,col=c(4,2,3),
       legend=c("Mean method","Naive method","Seasonal naive method"),bty="n")

Mean method: This method of forecasting takes an average of all the data and says that is to forecast the future performance of the firm. However, we determined that the data has a seasonal component, which this forecast does not account for.

Naïve method: This method takes the last data point observed, and sets that as all the future predictions. All predicts are one number, and because of this the seasonal component is not accounted for.

Seasonal Naïve method: The seasonal naïve method is very useful for seasonal data and sets each forecast to be equal to the last observed value from each quarter. This is especially useful is the dataset has a strong seasonal component, which indigo sales data demonstrates it has. This method is very versatile and can be used for data with nonlinear relationships, however, majority of the emphasis is put on seasonality which could prove to be an issue if there are other outside factors impacting the dataset.

Forecast using Holt-Winters smoothing methods
Indigofit4 <- hw(train, seasonal="multiplicative", h=h) 
Indigofit5 <- hw(train, seasonal="additive", h=h)

Smoothing Methods: It is evident that the best smoothing method to use as a prediction model is the Holt (1957) and Winters (1960) extended Holt’s method, which captures seasonality. There are two variations to this method, the multiplicative and additive method. For this time series the additive model is preferred because we observed earlier that we see roughly the same size peaks and troughs throughout the time series. However, for this research paper, I have forecasted using both methods to further verify that the additive method is more accurate.

Forecasting using Seasonal Decomposition of Time Series by Loess
y.stl <- stl(train, t.window=15, s.window="periodic", robust=TRUE)
Indigofit6 <- forecast(y.stl, method="rwdrift",h=h)

Seasonal Decomposition of Time Series by Loess Method: To forecast a decomposed time series, we have to determine if we want to include the seasonal component (which is unchanging) or adjust for the seasonal component. Earlier, we found that there were a fair bit of irregularities apparent in the data from 2012 to 2017, so the forecast has to adjust for seasonality. For time series decomposition, a random walk with drift model is the most appropriate (Seasonal and Trend decomposition using Loess).

Forecasting using ARIMA modelling and forecasting
adf.test(train, alternative = "stationary")
kpss.test(train)
ndiffs(train)

y.arima <- auto.arima(train)
Indigofit7 <- forecast(y.arima, h=h)

#Forecast by taking logs to reduce variability 
y.arima.lambda <- auto.arima(train, lambda=0)
Indigofit8 <- forecast(y.arima.lambda, h=h)
tsdiag(y.arima.lambda)

Auto-Regressive Integrated Moving Average Method This method is a class of model that captures a suite of different standard temporal structures in time series data. The ARIMA forecasting method shows residuals - both the ACF and the Ljung-Box statistic. Both of which show evidence of correlated residuals.When looking at the graphs, we notice only one lag as the line goes above the dotted blue line, but majority of the data shows autocorrelation.

Forecast using ANN
net <- nnetar(train)
Indigofit9 <- forecast(net,lambd=0,h=h)

Artificial Neutral Networks Method The ANN model is a mathematical model inspired by the function of the human brain and it is useful because it has powerful pattern classification and pattern recognition abilities. There is a seasonal pattern within our data, and this model may help to explain it further.

Forecast using BATS
tbats = tbats(train)
Indigofit10 <- forecast(tbats(train), h=h)

Box-cox form, ARIMA , Trend, and Seasonal components Method This is quite a good model to use for this data because it takes into consideration seasonality, which this data has a great deal of.

Forecast using VAR
# load Amazon and Indigo Date 
amazon_indigo_data <- read.csv("~/Dropbox (Personal)/ECON4210/Term Paper Stage Two/indigo_amazon_data.csv")

# indicate as time series
amazon_indigo_data = ts(amazon_indigo_data, start=2002, frequency=4)

y = amazon_indigo_data[,"indigo"]
cor (amazon_indigo_data[,2:3])

# Vector Autoregression
vardata = log (amazon_indigo_data[,c(3,2)])

VARselect(vardata, lag.max = 9, type = "both", season=4)
VARselect(vardata, lag.max = 9, type = "const", season=4)

Indigofit11 = VAR(vardata, p=3, type = "const", season =4)

serial.test(Indigofit11, lags.pt = 16, type = "PT.adjusted")
acf(residuals(Indigofit11), type="partial", lag.max=10)

# Granger Causality
causality(Indigofit11, cause= c("amazon"))
## $Granger
## 
##  Granger causality H0: amazon do not Granger-cause indigo
## 
## data:  VAR object Indigofit11
## F-Test = 2.0611, df1 = 3, df2 = 100, p-value = 0.1102
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: amazon and indigo
## 
## data:  VAR object Indigofit11
## Chi-squared = 3.1977, df = 1, p-value = 0.07374
# Impulse Responses
var1a.irf <- irf(Indigofit11,  n.ahead = 16, boot = TRUE,  runs=500, seed=99, cumulative=FALSE)

par(mfrow=c(2,2))
plot(var1a.irf, plot.type = "single")

par(mfrow=c(1,1))

par(mfrow=c(2,1))
plot( irf(Indigofit11, response = "indigo", n.ahead = 24, boot = TRUE,  runs=500) , plot.type = "single")

par(mfrow=c(1,1))

# Forcastings 
var.fc = predict(Indigofit11, n.ahead= 4)

VARS method

For the purpose of conducting a Vector auto regression forecasting model, I chose to compare indigo sales with amazon sales. That is because amazon has been known to booksellers like Borders Group Inc. out of business, as there has been a rise of e-reading. Also, it is important to note that Amazon controls about 7% of the Canadian online retail market, mainly due to its core line of book products. Indigo has reported to be struggling under the pressures facing all physical bookstores in the age of digital retail, as everything is going online. With this in mind, I want to see if amazon has been able to steal or possibly deter sales from indigo, which is why I chose to compare the two.

The null hypothesis of no Granger-causality shows a large p-value (0.1102> 0.05) which indicates weak evidence against the null hypothesis, so you fail to reject the null hypothesis. The null hypothesis of non-instantenous causality between amazon and indigo is rejected because it shows a large but still quite small p-value of (0.007373> 0.05) which indicates a strong evidence for the null hypothesis. Therefore, although you reject both null hypothesis, it is evident that “amazon go not granger-cause indigo” is much stronger, so it is more likely that amazon sales do not affect indigo sales. We will however, consider this further by analyzing the impulse responses.

Looking at the impulse response graphs, at no point does a change in the amazon sales seem to significantly shock indigo sales. Therefore, we can confidently conclude that although amazon and indigo are competitors, one’s performance does the significantly affect the others.

Indigo has had a hard time, however, they have been able to stand strong because indigo’s main attraction is not just books, they have a variety of neat gifts and generic merchandise. Moreover, they are able to offer product that amazon cannot, as many books are not available at a certain time on Amazon, especially new releases that have limited initial availability. Essentially, Indigo has been able to adapt well to change and stay competitive.

Plots of Forecasts
par(mfrow=c(2,3))
#####################
# Plot of Forecasts using Holt-Winters smoothing methods
plot(Indigofit5,ylab="Indigo sales",
     plot.conf=FALSE, type="o", fcol="white", xlab="Year", main= "Holt-Winters Seasonal Smoothing Forecasts")
lines(fitted(Indigofit4), col="red", lty=2)
lines(fitted(Indigofit5), col="green", lty=2)
lines(Indigofit4$mean, type="o", col="red")
lines(Indigofit5$mean, type="o", col="green")
legend("topleft",lty=1, pch=1, col=1:3, 
       c("data","Holt Winters' Additive","Holt Winters' Multiplicative"))

#####################
# Plot of Forcasts using Seasonal Decomposition of Time Series by Loess
plot(Indigofit6, main= "Forcast using Seasonal Decomposition of Time Series by Loess Method")

#####################
# Plot of Forecasts using ARIMA modelling
plot(Indigofit7, main="Forecasts from ARIMA Method")
plot(Indigofit8, main="Forecasts from ARIMA Method(log)")

######################
# Plot of Forecasts using ANN 
plot(Indigofit9, main="Forecasts using ANN Method")

######################
# Plot of Forecasts using BATS 
plot(Indigofit10, main="Forecasts using BATS Method")

par(mfrow=c(2,1))
######################
# Plot of Forecasts using VAR 
plot(var.fc, main = "Forecast using Vars Method")

Accuracy Measures
a1Indigofit1 = accuracy(Indigofit1, test)
a2Indigofit2 = accuracy(Indigofit2, test)
a3Indigofit3 = accuracy(Indigofit3, test)
a4Indigofit4 = accuracy(Indigofit4, test)
a5Indigofit5 = accuracy(Indigofit5, test)
a6Indigofit6 = accuracy(Indigofit6, test)
a7Indigofit7 = accuracy(Indigofit7, test)
a8Indigofit8 = accuracy(Indigofit8, test)
a9Indigofit9 = accuracy(Indigofit9, test)
a10Indigofit10 = accuracy(Indigofit10, test)
a.table <- rbind(a1Indigofit1, a2Indigofit2, a3Indigofit3, a4Indigofit4, a5Indigofit5, a6Indigofit6, a7Indigofit7, a8Indigofit8, a9Indigofit9, a10Indigofit10)

row.names(a.table)<-c('Mean training','Mean test', 'Naive training', 'Naive test', 'S. Naive training', 'S. Naive test' , 'Holt Winters multiplicative training', 'Holt Winters multiplicative test', 'Holt Winters additive training', 'Holt Winters additive test', 'STL training','STL test', 'ARIMA training', 'ARIMA test', 'ARIMA training (log)', 'ARIMA test (log)','ANN training', 'ANN test', 'BATS training', 'BATS test')

#Order the table according to MASE
a.table<-as.data.frame(a.table)
kable(a.table, caption="Forecast Accuracy",digits=3, bootstrap_options = c("striped", "hover", "condensed"),full_width = F, position = "left")
Forecast Accuracy
ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
Mean training 0.000000e+00 5.958800e+01 4.787200e+01 -6.432000e+00 2.141200e+01 4.819000e+00 -0.101 NA
Mean test 3.090800e+01 8.570600e+01 6.026200e+01 4.626000e+00 2.034200e+01 6.066000e+00 -0.153 8.400000e-01
Naive training 8.980000e-01 8.840700e+01 6.832700e+01 -5.869000e+00 2.934400e+01 6.878000e+00 -0.381 NA
Naive test 6.830800e+01 1.051480e+02 6.830800e+01 2.111000e+01 2.111000e+01 6.876000e+00 -0.153 1.019000e+00
S. Naive training 2.891000e+00 1.211400e+01 9.935000e+00 1.217000e+00 4.837000e+00 1.000000e+00 0.445 NA
S. Naive test 2.023100e+01 2.844800e+01 2.100000e+01 7.896000e+00 8.165000e+00 2.114000e+00 0.547 2.900000e-01
Holt Winters multiplicative training -2.045000e+00 9.742000e+00 8.100000e+00 -1.148000e+00 3.852000e+00 8.150000e-01 0.063 NA
Holt Winters multiplicative test 4.466900e+01 5.771500e+01 4.466900e+01 1.665100e+01 1.665100e+01 4.496000e+00 0.326 5.900000e-01
Holt Winters additive training -1.608000e+00 9.445000e+00 7.352000e+00 -1.090000e+00 3.549000e+00 7.400000e-01 0.047 NA
Holt Winters additive test 4.050600e+01 4.949000e+01 4.050600e+01 1.629000e+01 1.629000e+01 4.077000e+00 0.593 5.020000e-01
STL training 0.000000e+00 1.454100e+01 9.727000e+00 -4.310000e-01 4.499000e+00 9.790000e-01 -0.552 NA
STL test 2.880800e+01 3.972600e+01 2.933000e+01 9.966000e+00 1.024900e+01 2.952000e+00 0.059 4.010000e-01
ARIMA training -2.001000e+00 9.417000e+00 7.092000e+00 -1.269000e+00 3.378000e+00 7.140000e-01 -0.020 NA
ARIMA test 3.847400e+01 4.732300e+01 3.847400e+01 1.548400e+01 1.548400e+01 3.873000e+00 0.597 4.810000e-01
ARIMA training (log) -2.540000e+00 1.033700e+01 7.800000e+00 -1.353000e+00 3.598000e+00 7.850000e-01 -0.087 NA
ARIMA test (log) 4.005100e+01 5.116200e+01 4.005100e+01 1.508100e+01 1.508100e+01 4.031000e+00 0.357 5.220000e-01
ANN training 9.000000e-03 8.750000e+00 7.048000e+00 -2.090000e-01 3.433000e+00 7.090000e-01 0.281 NA
ANN test -1.251209e+76 4.232192e+76 1.251209e+76 -3.304156e+75 3.304156e+75 1.259423e+75 0.007 4.272033e+74
BATS training -1.276000e+00 1.095500e+01 8.624000e+00 -7.400000e-01 3.995000e+00 8.680000e-01 0.030 NA
BATS test 5.711600e+01 7.214500e+01 5.711600e+01 2.134300e+01 2.134300e+01 5.749000e+00 0.253 7.410000e-01

Comparing the accuracy for those model, it is evident that the seasonal naive model has the lowest amount of error and therefore is the best predictor for this data set.

Although it can be surprising, some forecasting methods are very simple and surprisingly effective. This method is the best one because it is versatile and can be used for data with nonlinear relationships. The data set have such a strong seasonal component which is the strongest factor, enough that our results have found that simply taking the previous years sales is the most accurate forecast.

From this research paper, I have learned that the sales data for Indigo shows seasonality, no apparent trend, a remainder component, and is best forecasted using the seasonal naive method. This forecast shows that Indigo sales stay relatively constant throughout time.