Introduction

The global cosmetic products market was valued at USD 532.43 billion in 2017 and is
expected to reach a market value of USD 805.61 billion by 2023, registering a Compound
annual growth rate (CAGR) of 7.14% during 2018-2023.

Hong Kong is the gateway to the mainland Chinese market. The beauty product sales
on the local market is not only limited to local consumers, but is driven by tourists, especially
those from the Chinese mainland (PR Newswire Association LLC,2018). According to
statistics, there are over 60 million tourists coming to Hong Kong every year, cosmetics and
skincare products are on top of their shopping lists. Brand availability, good services, product
authenticity and competitive prices are the main reasons why tourists buy such items in Hong
Kong. As no duty is applied on cosmetics and skincare products, prices are quite attractive to
visitors, especially for shoppers from the mainland. Thus, the skincare and cosmetics market
in Hong Kong is forecast to register high growth in value terms during 2015-2020 supported
by high disposable income of the consumer (Normans Media Ltd., 2016)

Majority of the beauty business companies in Hong Kong are traders. Many of them
act as agents for overseas manufacturers and they sell to wholesalers, retailers or re-export
the products to neighbouring countries like the Chinese mainland, Macau, Japan and
Southeast Asia. Traders in Hong Kong are beneficial partners for foreign brands as they have
extensive market knowledge, strong business skills, experience in trades with international
countries and strong networks. Therefore, Hong Kong continues to be a thriving cosmetics
and personal care market, with its traders benefiting from its strategic location as an
important gateway to mainland China (HPCi Media Ltd., 2017).

With a total national level sale estimated at $2.7 billion in 2015, Hong Kong’s
personal care and cosmetic products market is dominated by imports, as it has a very small
manufacturing sector of its own (HPCi Media Ltd., 2017). Thus, imported brands dominate
Hong Kong’s market for cosmetics, toiletries and skincare products. Hong Kong’s market for
such products is about 2.650 billion US dollars. The import statistics of Hong Kong identified
a steady increase over the last five years, with an increase of 15% between 2016 and 2017
alone. This unique phenomenon is also a matter of concern and active intervention by the
government, which has formed our basis to study the cosmetics sales index in Hong Kong.

Given the fluctuations of imported cosmetics, it creates a volatile environment in the

cosmetics sales index. Such a volatile environment in the sales of cosmetics creates several
problems. This could cause large fluctuations in the sale of cosmetics and being unaware of
them suggest that investors are exposed to extensive risk and losses. This is because,
investors would not be able to gauge the demand of cosmetics in Hong Kong and implement
strategical plans to profit from it. Furthermore, a volatility environment often creates
speculative bubbles that can cause investors and traders to panic and lead to destabilizing the
financial market and the economy.

Therefore, the motivation to forecast cosmetics sales index stems from the problems
mentioned above. By obtaining forecast of cosmetic sales, traders, investors and the
government can make informed decisions and discover policies to prevent huge financial
losses and instability in the economy. Thus, the forecast obtained in this report can solve
these problems and help with the decision-making process to improve the economy.
Phase 1: Model Identification

Firstly, package fpp2 and forecast was loaded by the library function. The data was imported
into R by read.csv function. One of the first tasks is required is plot the FULL data set with
the autoplot function as shown below

Figure 1: The monthly HK cosmetic retail sales value index

Based on time series plot, the x -axis is label as year and y -axis is label as retail sale index
from Oct 2014 to Sep 2015 = 100. The sales index expresses the current year's sales revenue
as a percentage of a base year's sales revenue. A sales index of over 100 indicates a year in
which sales exceeded the base period 's totals, while a number of less than 100 shows that the
current year's sales under-performed in comparison to the base period which is Consumer
Price Index of Oct2014 to Sep2015=100 for this data set. The graph reveals some interesting
characteristics such as there is a clear and increasing upward non-linear trend. There is also a
strong seasonal pattern that increases the variation increase overtime. The sudden increase of
sales at the start of each year was caused by dry weather in Hong Kong due to the winter
season. January is the coldest month where the average highest temperature is 19 °C and the
lowest is 14 °C. Therefore, average humidity in January is 77%, which is the lowest
throughout the year and to be considered uncomfortable (Boland,2018). Besides, air quality is
worse in winter because of temperature inversion, the warm air acts as a lid, covering air
pollutants beneath the cloud. Hence, more dust and air particles increase the probability of
getting pimples and blemishes, this increases the demand of cosmetics. Besides, after 2015
there is a slight decrease and after 2017 it increases again because Hong Kong’s economy
grew 2.4 percent in 2015, about half the pace of 2011, as a slowdown in mainland China and
a weaker yuan curbed Chinese spending, while a volatile stock market also hit domestic
consumption(Chatterjee & Pan, 2016)

To build the ARIMA model, stationarity in data is crucial for model identification step, that
helps us to build ARIMA(p,d,q) model. A weak stationary time series means;

• 𝐸(𝑦𝑡 ) = 𝜇 for all t → mean is time invariant

• 𝑉𝑎𝑟(𝑦𝑡 ) = 𝐸[(𝑦𝑡 − 𝜇)2 ] = 𝛾0 (a finite constant) for all t and j

• 𝐶𝑜𝑣(𝑦𝑡 , 𝑦𝑡−𝑗 ) = 𝛾𝑗 for all

• No patterns predictable in the long run

Figure 2 The training HK cosmetic sales value index data graph

Firstly, we partition the data into a training and test set. The training set contain 100 data
from Jan 2008 to April 2016 and test set contain 39 data from May 2016 to July 2019. This is
done by allocating 72% of the total number of observations into the training set and allocate
28% of the total number of observations into the test set. Based on observation of training set
data plot, there is an upward trend and seasonality is present in the data. Besides that, the
variation of the seasonality is increasing with the level which means the variance is not
constant and mean is not constant. This means the data is non stationary and transformation is
necessary.
To further prove on the data is non-stationarity, we used the ggtsdisplay function. This shows
an ACF graph that is decreasing slowly and the value of 𝑟1 is large and positive shown in
appendix 1. To change the data into stationary series, we identified the lambda of the training
data using the BoxCox.lambda function. The lambda we have obtained is 0.4054653 as
shown in appendix 2. Then we stabilize the variance using the Box-Cox transformation with
the lambda we have obtained.

Next, both seasonal and first differences are applied, it is suggested that if the seasonality is
strong seasonal differencing is done first because sometimes after seasonal differencing
stationary series is obtained and first differencing is no needed. Based on function nsdiffs
and ndiffs which serve the purpose of using a sequence of KPSS tests to determine the
appropriate number of first differences is carried out by the function, both differing is
required once to stabilize the mean as shown by the raw output by R in appendix 3.

Eyeballing on the plot in appendix 4, the data looks stationary now because it has a constant
variance and a constant mean.

To further prove it, we conducted KPSS test to test if our data is stationary. The test statistic
(0.0483) obtained is lower than the 1% critical value shown in appendix 5, indicating that we
do not reject the null hypothesis and conclude that the data is stationary. Since the data is
stationary, we can move forward to the model identification step.

From the figure 4 above, we analyzed the ACF and PACF of the model and obtain the
following:
1. We identified an AR (2) process because there is significant spike at lag 2 in
PACF
2. We observed a MA (1) process because there is a significant spike at lag 1 in the
ACF
3. We identified a Seasonal AR (1) process because there is a significant spike at lag
12 in the PACF. The spike in 12 also means monthly data is used.
4. We observed a Seasonal MA (1) process because there is a significant spike at lag
12 in the ACF. The spike in 12 also means monthly data is used.

We identified the 𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 process by using the MA(1) and seasonal MA(1)
process. Other than that, we identified the 𝐴𝑅𝐼𝑀𝐴(2,1,0)(1,1,0)12 by using the AR(2) and
seasonal AR(1) process. Furthermore, we identified the 𝐴𝑅𝐼𝑀𝐴(2,1,1)(1,1,1)12 by using the
AR(2), MA(1), seasonal AR(1) and seasonal MA(1) process. All specified ARIMA models
should not have a constant because if c ≠0 and d + D = 2, the long-term forecasts will follow
a quadratic trend, this is not accurate because we would expect the long-term forecasts will
follow a straight line with intercept and slope determined by the last few observations where
c = 0 and d + D = 2.

Phase 2: Estimation and testing

Based on appendix 6, we estimated the identified ARIMA models from Phase 1 using the
training set data by specifying the orders of each ARIMA model, constant should not be
included and the lambda that R estimated on the training set. Besides that, we used the
function auto.arima by setting the stepwise and approximation argument as FALSE. With this
function and arguments, R will return the best ARIMA model according to either AIC, AICc
or BIC value. The function conducts a search over possible model within the order
constraints provided. Furthermore, the arguments allow R to do selection more accurate with
excessive computation times. The parameters estimate of the ARIMA models and relevant
statistic is reported on a Table 1 below

Table 1: Parameters estimated and information criteria

Rank the ARIMA models

RANK AICc
𝐴𝑅𝐼𝑀𝐴 (2,1,0) (0,1,1)_12 1 27.58
𝐴𝑅𝐼𝑀𝐴 (0,1,1) (0,1,1)_12 2 30.14
𝐴𝑅𝐼𝑀𝐴 (2,1,1) (1,1,1)_12 3 32.04
𝐴𝑅𝐼𝑀𝐴 (2,1,0) (1,1,0)_12 4 33.98

Table 2: Ranking of ARIMA models based on AICc

We tabulated and ranked the models base one information criteria AICc. Based on table 2
above, the 𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 model and the 𝐴𝑅𝐼𝑀𝐴(2,1,0)(0,1,1)12 are our top 2
models estimated based on the AICc because it has a lower AICc compared to the other
ARIMA models. Although AICc sometimes being much more difficult to compute than AIC
it is much preferred because AICc has the advantage of tending to be more accurate than AIC
especially for small samples.

2𝑘 2 + 2𝑘
𝐴𝐼𝐶𝑐 = 𝐴𝐼𝐶 +
𝑛−𝑘−1

k = the number of parameters

n = the sample size

If the candidate models have the same k and the same formula for AICc, then AICc and AIC
will give identical (relative) valuations; hence, there will be no disadvantage in using AIC,
instead of AICc. However, based on the ARIMA models we are comparing models with
different numbers of parameter hence AICc is preferred. Furthermore, we select AICc rather
than AIC or BIC because refer to the Help function for the auto.arima() function, the option
'aicc' appears as the first option for the argument ic. Therefore, unless an alternative
information criterion is specified for the model selection, the default setting will be 'aicc'.

Moreover, in order to check the top 2 models are white noise we performed diagnostic checks
on the residuals with the function of checkresidual which produces a time plot of the
residuals, the corresponding ACF, and a histogram and automatically perform Ljung-Box
test.

Based on the Ljung-Box test in appendix 7.1, the 𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 model has a p-
value higher than alpha. Therefore, we do not reject the null and conclude that the residuals
are white noise. This makes sense because at least 95% of the ACF lags are not significant in
the ACF lags of the residuals.

Based on the Ljung-Box test in appendix 7.2, the 𝐴𝑅𝐼𝑀𝐴(2,1,0)(0,1,1)12 has a p-value
higher than alpha. Therefore, we do not reject the null and conclude that the residuals are
white noise. This makes sense because there are no significant spikes in the ACF lags of the
residuals.

Thus, both 𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 model and 𝐴𝑅𝐼𝑀𝐴(2,1,0)(0,1,1)12 models are adequate.

Also, both 𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 model and 𝐴𝑅𝐼𝑀𝐴(2,1,0)(0,1,1)12 model does not need
to be reidentified because the residuals of both models are white noise.
Since reidentification is not needed, we use the accuracy function to obtain range of summary
measures of the forecast accuracy as shown in a table below.

ARIMA MODELS SETS MAPE MASE RMSE

Training set 2.889028 0.270137 3.38646
𝐴𝑅𝐼𝑀𝐴 (0,1,1) (0,1,1)_12
Test set 4.037557 0.509737 5.270988
Training set 2.78516 0.257853 3.189602
𝐴𝑅𝐼𝑀𝐴 (2,1,0) (0,1,1)_12
Test set 4.207067 0.520298 5.167254

Table 3: Accuracy measures of ARIMA models

Based on Table 3, we can identify that 𝐴𝑅𝐼𝑀𝐴 (2,1,0)(0,1,1)12 model has a better goodness
of fit compared to 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 model because the MAPE and MASE of the
training set for 𝐴𝑅𝐼𝑀𝐴 (2,1,0)(0,1,1)12 model is lower than 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 model.
However, the 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 model has a better forecast accuracy compared to
𝐴𝑅𝐼𝑀𝐴 (2,1,0)(0,1,1)12 because the 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 model has a lower MAPE and
MASE for the test set compared to 𝐴𝑅𝐼𝑀𝐴 (2,1,0)(0,1,1)12 model. On the other hand,
RMSE concluded that 𝐴𝑅𝐼𝑀𝐴 (2,1,0)(0,1,1)12 is a better in goodness of fit and has a better
forecast accuracy with lower RMSE for training and test set. Nevertheless, the mean absolute
percentage error (MAPE) is preferred because it is unit-free. There are no values close to
zero. Therefore, we do not need to worry about extreme values and values being infinite.
Thus, 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 model is a better forecasting model.

Backshift notation of 𝑨𝑹𝑰𝑴𝑨 (𝟎, 𝟏, 𝟏)(𝟎, 𝟏, 𝟏)𝟏𝟐

12 )
(𝑦𝑡𝜆 − 1)
(1 − 𝐵)(1 − 𝐵 ( ) = (1 + 𝜽1 B)(1 + 𝞡𝟏 𝐵12 )𝜀𝑡
𝜆
Phase 3: Application

In phase 3, we re-estimated the full data set base on 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 model that we
concluded in phase 2. Arima function was used again with the same codes as estimating the
train set with 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 except the data set was FULL data set as shown in
appendix 8. Backshift notation is used to report the estimated model as shown below.

BACKSHIFT NOTATION FOR ARIMA (0,1,1)(0,1,1) of entire DATA

(𝑦𝑡0.405 − 1)
(1 − 𝐵)(1 − 𝐵12 ) ( ) = (1 − 0.5464B)(1 − 0.8834𝐵12 )𝜀𝑡
0.405

Furthermore, the estimated model is used to produce a h-step ahead forecast, along with their
corresponding 95% forecast intervals, h = 24 because monthly data is used, and it is plotted
as shown below

Figure 3 above is the plot of the 24 steps ahead forecast of 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 model
using the entire dataset.

Based on the graph above the pattern of the forecast follows the seasonal pattern of the entire
data set. Additionally, the forecasted result shows that retail sales are increasing steadily with
a stable variance which is reasonable because the real data set shows deterministic upward
trend.

Also, the prediction interval is relatively tight and does not variate far. Thus, the forecast
interval is reasonable. Although our prediction interval is relatively tight, there are some
concerns that arise from our forecast. This is because, due to the Anti-Extradition Law
Amendment Bill movement in Hong Kong, this will affect the sales index of cosmetic
products in Hong Kong. Therefore, the forecast and prediction intervals we have concluded
might not be accurate.

Comparison with an ETS model

ETS function is used in this section to let R select the ETS model automatically the model.
Next, the estimated ETS model is used to forecast 12 months ahead forecast as shown in
appendix. Both 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 and ETS model forecast is plotted in a graph below.

Figure 4 above is the plot of the 24 steps ahead forecast of 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 model
and ETS model using the entire dataset.

In comparison between the 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 model that we selected with the ETS
model that R selected. Based on figure 8, we can conclude that Arima model provide a more
reasonable forecast because the it follows the upward trend from the original data set.
Besides, the seasonality is captured by 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 with a reasonable variation.
On the other hand, ETS estimated a rapid decrease in sales for the following two years which
is inaccurate based on the FULL data set which shows a strong upward deterministic trend.
Conclusion

In conclusion, we have selected Hong Kong cosmetic sales index as the data to
forecast and converted it to a time series. Next, we plot the data to identify the features of the
plot. After that, we partitioned the data into training and test set and plot the training set.
From the plot of the training set, we established that the data was not stationary and needed to
undergo transformations to stabilize the mean and variance. We used a Box-Cox
transformation to stabilize the variance and we used a seasonal difference and a first
difference to stabilize the mean. After that, we used the ggtsdisplay function in R to plot the
transformed data along with the ACF and PACF of the data. From the ACF and PACF of the
transformed training data, we identified 3 possible ARIMA models which are
𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 , 𝐴𝑅𝐼𝑀𝐴(2,1,0)(1,1,0)12 and 𝐴𝑅𝐼𝑀𝐴(2,1,1)(1,1,1)12. Next, we
estimated the 3 ARIMA models and used the auto.arima function in R to estimate another
ARIMA model. R had generated an 𝐴𝑅𝐼𝑀𝐴(2,1,0)(0,1,1)12model from the auto.arima
function. Subsequently, we ranked the estimated ARIMA models based on the model’s AICc
and selected the top 2 models to compare. The top 2 models that was selected was the
𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 and 𝐴𝑅𝐼𝑀𝐴(2,1,0)(0,1,1)12 . We then proceeded to perform
diagnostic checks on the residuals of the selected models. Both models were adequate
because their residuals were white noise based on the Ljung-Box test. After that, we
forecasted both models and selected the model with the best forecast accuracy against the test
set. We compared both models using the MAPE and MASE of the test set to determine which
model had the best forecast accuracy. We selected the 𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 model
because it had a lower MAPE and MASE of the test set compared to the
𝐴𝑅𝐼𝑀𝐴(2,1,0)(0,1,1)12 model. After that, we re-estimated the 𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12
model using the entire dataset and produce a 24-step ahead forecast along with the
corresponding forecast interval. Subsequently, we plot the forecasted data and identified that
the prediction interval is relatively tight and does not variate far. Thus, the forecast interval is
reasonable. We then proceeded to forecast from the entire data set using ETS. We then plot
the forecast of the 𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 model and the ETS model on data to compare
them. From the plot, we have identified that the 𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 had a better fit
because it was following the upward trend of the data while the ETS model was decreasing
rapidly. Hence, we concluded that 𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 is a better model.
 Lastly, the limitations of using an ARIMA model is that it is high in cost. Because of
the large data requirements, the lack of convenient updating procedures, and the fact that they
must be estimated using nonlinear estimation procedures, the Box Jenkins models tend to be
high cost. Other than that, another limitation of using an ARIMA model is that it is hard to
interpret coefficients or explain “how the model works” and there is a danger of overfitting or
mis-identification if not used with care.
Reference

• Analysis of Hong Kong's Imports of Goods | Census and Statistics Department.

(2019). Retrieved 26 October 2019, from
https://www.censtatd.gov.hk/hkstat/sub/sp230.jsp?productCode=FA100016
• Beauty & Personal Care - Hong Kong | Statista Market Forecast. (2019). Retrieved 26
October 2019, from https://www.statista.com/outlook/70000000/118/beauty-personal-
care/hong-kong
• Beauty and Personal Care in Hong Kong, China. (2018).
• Boland, R. (2018, October 28). Average Humidity in Hong Kong Month by Month.
Retrieved from https://www.tripsavvy.com/average-humidity-in-hong-kong-month-
by-month-1535877.
• Consumer and Market Insights: Skincare in Hong Kong Citation metadata. (2016).
• Kastner, J. (2017). Can Hong Kong retain Chinese visitor interest? Hong Kong's
cosmetics and personal care sectors are thriving. But the territory risks losing
mainland Chinese consumers to the lure of Korea and Japan Citation metadata. Soap
Perfumery & Cosmetics, 90(11).
• Lam, J., & Leinders, K. (2018). COSMETIC AND SKINCARE PRODUCTS IN
HONG KONG. Retrieved 26 October 2019, from
https://www.flandersinvestmentandtrade.com/export/sites/trade/files/market_studies/F
ITcosmeticsReportHongKong_1.pdf.
• Pan, T., & Chatterjee , S. (2016, May 12). https://www.reuters.com/article/us-
hongkong-economy-gdp-idUSKCN0Y316X. Retrieved from
https://www.reuters.com/article/us-hongkong-economy-gdp-idUSKCN0Y316X.
APPENDIX

Appendix 1: ggtsdisplay() function for training set data

Appendix 2: Lambda estimated by R using the train data set

Appendix 3: Output of function nsdiffs and ndiffs

Appendix 4: ggtsdisplay function on stationary training data

Appendix 5: (KPSS) test results

Appendix 6.1 above shows the estimation of an 𝐴𝑅𝐼𝑀𝐴(2,1,1)(1,1,1)12 model using the
training set labelled fit1.
Appendix 6.2 above shows the estimation of an 𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 model using the
training set labelled fit2.

Appendix 6.3 above shows the estimation of an 𝐴𝑅𝐼𝑀𝐴(2,1,0)(1,1,0)12 model using the
training set labelled fit3.

Appendix 6.4 above shows the estimation of an 𝐴𝑅𝐼𝑀𝐴(2,1,0)(0,1,1)12 model obtained

using the auto.arima function on the training set labelled fit4
Appendix 7.1 above shows the diagnostic checks on the residuals for the
𝐴𝑅𝐼𝑀𝐴(0,1,1)(0,1,1)12 model.

Appendix 7.2 above shows the diagnostic checks on the residuals for
the 𝐴𝑅𝐼𝑀𝐴(2,1,0)(0,1,1)12.
Appendix 8 above shows the estimated 𝐴𝑅𝐼𝑀𝐴 (0,1,1)(0,1,1)12 model using the entire data
labelled fit 5.