1. INTRODUCTION
Debates on agriculture in Vietnam have been familiar with buzzwords such as rice or coffee, however, in recent years, the cacao production is slowly becoming prominent. In 2016, Vietnam joined Indonesia to become the only two fine flavour cocoa nations in Asia, granted by The International Cocoa Organization (ICCO). This can be considered as one big advantage since the consumption on premium cacao products, which is mainly made from flavour cacao, is increasing globally. According to the VCC, in Vietnam, cacao crops cover the area of 11,559 ha. There are projects to improve imoverished soil resulted from poor coffee and rubber plantations and turn them into cacao crops accross the Southeast and the Central Highlands.
In 2017, it is reported by the ICCO that cacao production had suffered a surplus due to the downturn in price. The surplus within the period of 20162017 reached 371,000 tones, in spite of the fact that the previous season was a success to this industry with the deficit of 174,000 tons of output. Nevertheless, leading financial institutes expect to observe a recovery of this market regarding the 20172018 season. The drop in price should play an incentive role to boost the demand and restore the marker equilibrium.
In the other words, the manufacturing and crops in the fields of agricultural products is one of the most strength of Vietnam, however, the challenges came from the globalization and technological revolution (Hao et al., 2018). Moreover, the other disadvantages of famers in rural area in Vietnam are that the manufacturers and factories are mostly small and mediumsized enterprises (SMEs), and “it is difficult for them to take quick actions and relocate resources flexibly based on turbulent requirements. Therefore, it is critical to move towards sustainable development in terms of business models and business practices.” (Hao et al., 2018).
Demand forecasting is claimed as the main trigger function of the production planning and control process (Kazemi et al., 2018). It attempts to calculate future scenarios, conditions, and contexts, offering the possible estimation of accessible market information (Rabbani et al., 2017). Besides, the global supply network is exposed to a variety of risks, including supply disruption, supply delays, demand fluctuations, price, fluctuations and exchange rate fluctuations (Soni and Kodali, 2013). Most of the time, the effect of mentioned risks on supply chain network is carried out after the realization of facilities and resources (Prakash et al., 2017). On the other hand, in present context of globalization, customer requirements are changing continuously and product life cycle is shortening (Singh et al., 2017). Moreover, recent studies demonstrated that an accurate demand forecast has a direct impact on customer satisfaction, safety stock levels, and total cost (Gharaei et al., 2018a, 2018b, 2019a, 2019b, 2019c).
According to a report on food wastage footprint published by the Food and Agriculture Organization of the United Nations, roughly 35% of all food produced for consumptionis lost or wasted (reported by FAO, 2013). Annually, this wastage economically costs us about USD 1 trillion. Nevertheless, it also contributes to the degradation of our environment, affects the society’s wellbeing and livelihoods. Particularly, in South and Southeast Asia, the per capital food loss was 125kg/year, of which the food loss due to the supply chain accounted for 90% (Figure 1). It also comes to light that in lowincome countries, a significant amount of food is wasted during the supply chain process (FAO, 2013). On the other hand, the Department for Environment Food and Rural Affairs stated clearly that one of the reasons that create food waste is poor demand forecasting (Rabbani et al., 2018).
In 2017, Vietnam experienced the surplus of pork, fruits and many other kinds food which caused decreases in prices vary from 20%50% and several damages to farmers and manufactures. The Government had to provide significant financial support to farmers in order to deal with this surplus. Either at a micro or macro point of view, overproduction creates losses and violates the equilibrium of the market, then eventually has negative impact on other aspects. In order to maintain good conditions of the systems, we have to decide which components should be given maintenance activities in the limited time allotted between missions (Dao et al., 2014). In particularly, the cacao crops can produce fruits all year round. This could be risky as if the famers and the manufacturer do not work on prediction for the upcoming season. Failure in forecasting may cause the famers to overproduct and consequently create considerable damage to the economy as a whole. Hence, an inappropriate tool for forecasting demand would jeopardize the supply chain and production outcomes (Gharaei and Pasandideh, 2017).
This research aims to examine different forecasting methods for a cacao manufacturer then picking out the most suitable one for longterm usage. It can also be considered as suggestions for a point of view on production planning for agriculture product in general.
Puratos Grand – Place Indochina is a joint venture between Puratos Group and GrandPlace Holding that was established in 2012. Within the period of 20172019, PGPI plans to extend its business to the North region of Vietnam. Despite the fact that it aims to manage the large market of Indochina area, there is no accurate forecasting tool implemented to improve the efficiency of supply chain and production. Regularly, they purchase an amount of cacao beans from its own cacao plantation located in Ben Tre twice a month. This amount can be varied depends on the productivity of the crops. However, the productivity level is complicated to measure and not stable since there are several external factors that could affect the harvest, for instance, the weather, and the quality of fertilizer. The problem was addressed seriously when the power cut happened more than usual. For cacao manufacturer, storing the finish goods is probably the most important process within the SCM since running the cold room for storage is costly. Besides, it is not environment – friendly to use the cooling machine because of the CFC.As the market is volatile and competitive, the firm calls for cacao powder demand prediction in order to reduce any losses occur from inaccurate production planning.
2. LITURATURE REVIEW
2.1 Demand Forecasting
According to Sobhanallahi et al. (2016a, b), the changes and challenges that maybe introduced in small increments could be handled just by alone experience. Nevertheless, business environment is characterized as a highly competitive, dynamic and volatile market (Kumar Sharma and Bhat, 2014); comparing with the past, recent supply chain (SC) becomes more complex and expanding dramatically, that is the reason why “the industries and businesses need the scientific planning, modeling and optimizing in order to remain competitive which can be achievable by supply chain management (SCM)”.
A forecast of demand is the prediction of seasonal data and trend based on historical data. Accordingly, many manufacturing firms have tried various approaches to remain competitive (Mahmood et al., 2011). Based on the above mentioned issues , if companies and organizations intend to adapt their business models and environmental strategies to the mutable and competitive market, they must be extremely flexible and employ appropriate SCM which can be attained by optimizing accurate models in the business environment. It sets the foundation for Supply Chain Management, helping the manufacturers to identify the level of the product’s consumption by the customers within a specific period, then come out with an optimal production plan. Due to the increasing complexity exists within the marketplaces; manufacturers are required to develop their awareness as to future scenarios. Thus, over the last few decades, several theoretical and heuristic methods had been tested and emerged to improve the ability of having an accurate forecast (Hoseini Shekarabi et al., 2018).
Otherwise, flexibility in SC means to maintain the customer service levels adapting disturbances in supply and sudden changes in demand (Kumar et al., 2008). Research studies on demand forecast have often presumed that the level of accuracy is crucial in order for one to select the ideal model. A failure in forecasting can lead to several possible damages including high cost of storage, purchasing material, disposal goods (Nguyen and Tran, 2019).
2.2 Demand Forecasting in Food Industry
Food manufacturers have to overcome the challenge to provide the right quantity of goods at an excellent quality and reasonable prices.Based on the paper of Pinçe et al. (2008), it analysed an inventory system facing stochastic external demands and an autonomous supply under a continuous review for the purpose of replenishment– disposal policy. In addition, they learned that customer satisfaction would be increased under circumstances where the inventory is adjusted rationally. On the other hand,optimizing the supply chain is the key in cutting cost but to decide a level for production is complicated. In addition, Sharma and Henriques (2005) have argued that we would reduce spending with the significant savings in cost. Then, manufacturers must have an accurate prediction on the demand in order not to put themselves in a position where they are either running out of stock or overstocked (Pasandideh et al., 2015).
However, the industrial logistics for food and beverage is different from other products due to certain specific characteristics. The expiration day of the good, especially perishable good, is one of those most important characteristics that must be tackled (Prusa and Chocholac, 2015). Studies have been conducted to examine different demand prediction methods, yet by no means shall the findings recognize the existence of a consensus. In order to determine an optimal forecasting model, one must implement performance indicators to evaluate the pros and cons of several methodologies giving initial conditions. In this dissertation, the researcher applies forecasting methods to suggest the suitable directions for demand forecasting in the context of a dominant cacao manufacturer (Gharaei et al., 2017;Awasthi and Omrani, 2019;Duan et al., 2018). The production line operates under the maketo stock system since there is low level of customization.
Some traditional time series techniques such as moving averages, multiple regression analysis (MRA) and autoregressive integrated moving average (ARIMA) are referable by researchers and companies because of the fact that they are particularly friendly in terms of interpretation and implementation. For instance, Shukla and Jharkharia (2013) used ARIMA model to predict the daily sales for a whole sale agriculture market and the mean absolute percentage error (MAPE) for the case was approximate 30%.
On the other hand, not using familiar methods, a recent investigative study in the field of food industry indicated that the Neutral Networks models effectively presents a superior forecasting result for demand (Gharaei and Pasandideh, 2017b). Meanwhile, Nguyen and Tran (2018) presented the use of Grey Model – GM (1, 1) in analyzing the trend of quality index for food in the future. In addition, a research on prediction of grain yield conducted by a group of researchers from China demonstrate that among other methods such as DGM(1, 1), the combined model, a multiple linear regression model, GM(1, N) model is claimed to have a highest application value for this case. Therefore, within the scope of agriculture product, the GFMs has shown a significant contribution in term of forecasting.
2.3 Methods of Forecasting Demand
2.3.1 Arima
The Autoregressive Integrated Moving Average (ARIMA) models for forecasting time series are essentially agnostic. The ARIMA methods are not assumable with background of any underlying economic model or structural relationships. It is assumed that past values of the series plus previous error terms contain information for the purposes of forecasting (Nguyen and Tran, 2018).
The number of times in the series must be differenced to induce the stationary of the integrated component with the (p, d, q)(P, D, Q) as the notations following.

p: the number of autoregressive terms,

p: the number of moving average terms,

d: the number of times a series.

P: the number of seasonal autoregressive components

Q: the number of seasonal moving average terms

D: the number of seasonal differences”
This may be written as:
where,
${X}_{t}={\nabla}^{d}{\nabla}_{s}^{D}{Y}_{t}$ is a stationary series,
${\nabla}^{d}={(1B)}^{d}$ represents the number of regular differences, and ${\nabla}_{s}^{D}={(1{B}^{s})}^{D}$ represents the number of seasonal differences required to induce stationarity in Y_{t}, s is the seasonal span (hence for quarterly data s = 4 and for monthly data s = 12), B is the backshift operator (such that ${B}^{0}{X}_{t}={X}_{t},\hspace{0.17em}{B}^{1}{X}_{t}={X}_{t1},\hspace{0.17em}{B}^{2}{X}_{t}={X}_{t2},\hspace{0.17em}\dots .$), $\theta \left(B\right)=1+{\theta}_{1}B+{\theta}_{2}{B}^{2}+\dots +{\theta}_{q}{B}^{q}$ is a qorder polynomial in the backshift operator,
Followings are some come models with their parameters which can help the forecasting process accurately.
ARIMA(0,1,0) shows random elements:
ARIMA(1,1,0) (AR(1)) shows differenced firstorder autoregressive model:
or
ARIMA(0,1,1) (MA(1)) without constant means simple exponential smoothing:
where e(t1) denotes the error at period t1. Note that this resembles the prediction equation for the ARIMA(1,1,0) model.
ARIMA(0,2,1) or (0,2,2) without constant stands for linear exponential smoothing:
or
2.3.2 Exponential Smoothing
The SES model is given by the model equation
Or can be explained as
where ε_{t} is the forecast error (actual forecast) for period t.
In other words, the new forecast is the previous one adjusted for the error that occurred in the last forecast.
2.3.3 Basic Grey Model – GM (1,1)
In 1982, Deng proposed a system theory that is claimed to overcome the problem of limited data. He named a limited valuable sequence data as grey information. He also came up with a technique called grey forecasting model that can use grey information as input to calculate future value. The model handles original data through an accumulating generation operator (AGO) to find the inherent structure in a system. Moreover, it is not necessary to make any other assumptions on the statistical distribution when employing this method.
When the Grey Forecasting method was first mention, there were two models established by Professor Deng. One was the GM (1,1) which is first–order and single–variable grey dynamic model, and the other is first – model and multivariable grey dynamic model or GM (1,N). After that, more and more novel Grey Forecasting Models have been developed. The most applications of GFMs are agriculture, environment and marketing. Yet there are some drawnbacks, for instant, discrete, larger than 4 elements and positive data only
The Figure 2 below was created by Wang et al. (2012) can be considered as a good explanation on how the general process of GMFs solves the data input. It is summarized in seven steps: variable selection, modeling sequences develop, model choosing, parameters calculation, error analysis, property analysis and forecasting.
The reason behind choosing the GFMs for the given case of demand forecasting at Puratos GrandPlace Indochina is the ability of dealing with poor data set. In fact, the Vietnamese market contains a considered amount of uncertainty, means that there are a lot of unpredicted factors can affect the trend and causing outliers. In addition, the original data was in the form of cacao powder mixture, which means it not the pure cacao powder that here we are trying to predict the demand. Therefore, data has gone through a process where it is converted into the sale of pure cacao powder based on the exact proportion of ingredients. This processmight generate some violation to the numbers.
The basic grey model, GM(1,1), is known as the most conventional and basic approach in grey system theory. It is developed for time series forecasting with limited scale of data.
The procedure of GM(1,1) is described below:

Collect a time series data sequence

Convert the raw series X^{(0)} into a series X^{(1)} using the AGO
$${X}^{(1)}=\{{x}^{(1)}(1),\hspace{0.17em}{x}^{(1)}(2),\hspace{0.17em}\mathrm{...},\hspace{0.17em}{x}^{(1)}(1)\},\hspace{0.17em}\hspace{0.17em}{x}^{(1)}(1)={x}^{(0)}(1)$$where ${x}^{(1)}(k)={k}_{i=1}{x}^{(0)}(i),\hspace{0.17em}k=2,\hspace{0.17em}3,\hspace{0.17em}\mathrm{...},\hspace{0.17em}n$

Determine the background values z^{(1)} :

Estimate the developing coefficient a by equation (3), then calculate the gray input b
In above, [a, b]^{T} is a sequence of parameters that can be evaluated as follows:
where
Theorem 2: Use the result for estimated coefficients a and b together with the condition
to solve equation (3)
then, we obtain the predicted value at time (k+1) as equation:
2.4.4 Verhulst
As we noticed, Grey model usually cause errors due to exponential sequence. In order to increase the quality of simulated outcomes, a modified Grey Verhulst Model was applied, using the same prediction model mechanism as Basic Grey Model. For the differences, Verhulst model finds an inherent error, then it recommends a new method to construct an optimized derivative, then news parameters are used or old parameters could be reused by equations which meet the property of un_biasedness.

Denoted x^{(0)} as nonnegative series, x^{(1)} as the 1 AGO series of x^{(0)}, z^{(1)} as the mean series of x^{(1)}, we call ${x}^{\left(0\right)}\left(k\right)+a{z}^{(\text{1})}\left(k\right)=b{\left({z}^{(\text{1})}\left(k\right)\right)}^{\text{2}}$ be the grey differential equation of Verhulst model.

Let x^{(0)} be nonnegative series, x^{(1)} be the 1 AGO series of x^{(0)}, we call $\frac{d{x}^{(1)}}{dt}+a{x}^{(1)}=b{({x}^{(1)})}^{2}$ be the whitened differential equation of grey Verhulst model.
Theorem 1: Let x^{(0)} be nonnegative series, x^{(1)} be the 1 AGO series of x^{(0)}, z^{(1)} be the mean series of x^{(1)} and α = (a, b)^{T} be the parameter of the equation, let
then the least squares estimate of the parametric sequence $\widehat{\alpha}={\left[a,b\right]}^{T}$ of the equation ${x}^{\left(0\right)}\left(k\right)+a{z}^{(\text{1})}\left(k\right)=b{\left({z}^{(\text{1})}\left(k\right)\right)}^{\text{2}}$ is $\widehat{\alpha}={({B}^{\text{T}}B)}^{1}{B}^{\text{T}}Y$
The solution of x^{(0)}(k) + a z^{(1)}(k) = b(z^{(1)}(k))^{2} is
The time response sequence of the grey Verhulst model is:
2.3.5 Discrete Grey Model
Dealing with exponential sequence usually cause errors with GM(1,1). The application and characteristic of Discrete Grey Model  DGM(1,1) are similar to GM(1,1). But being modified with optimization technique applied, DGM(1,1) could be more appropriate for simulate an exponential pattern set of data.
Discrete grey model abbreviated as DGM Theorem 1:
If $\widehat{\beta}={({\beta}_{1},\hspace{0.17em}{\beta}_{2})}^{\text{T}}$ is a sequence of parameters, and
then the least squares estimate sequence of the grey equation
satisfies
Proof
Substitute data values into the grey differential equation
gives that
That is, $Y=B\widehat{\beta},$ for the evaluated values of β_{1} and β_{2} , substitute ${x}^{(1)}(k+\text{1}),\hspace{0.17em}k=\text{1},\hspace{0.17em}\text{2},\hspace{0.17em}\dots ,\hspace{0.17em}n1$ with ${\beta}_{\text{1}}{x}^{(1)}(k)+{\beta}_{\text{2}}$, which gives the error sequence $\epsilon =YB\widehat{\beta}$. Let
The values of β_{1}, β_{2} make S minimum and should satisfy
Solve the equations, then we can get
and
from $Y=B\widehat{\beta}$, it follows that ${B}^{\text{T}}B\widehat{\beta}={B}^{\text{T}}Y,\hspace{0.17em}\widehat{\beta}={({B}^{\text{T}}B)}^{1}{B}^{\text{T}}Y$.
But
and
Therefore,
“Theorem 2
Assume that $B,\hspace{0.17em}Y,\hspace{0.17em}\widehat{\beta}$ are the same as Theorem 1. That is $\widehat{\beta}=[{\beta}_{1},\hspace{0.17em}{\beta}_{2}]={({B}^{\text{T}}B)}^{1}{B}^{\text{T}}Y$.
Then the following holds true

Set x^{(1)}(1) = x^{(0)}(1), then recursive function is given by ${\widehat{x}}^{(1)}(k+1)={\beta}_{1}^{k}{x}^{(0)}(1)+\frac{1{\beta}_{1}^{k}}{1{\beta}_{1}}\cdot {\beta}_{2},\hspace{0.17em}k=1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}n1$ or ${\widehat{x}}^{(1)}(k+1)={\beta}_{1}^{k}\left({x}^{(0)}(1)\frac{{\beta}_{2}}{1{\beta}_{1}}\right)+\frac{{\beta}_{2}}{1{\beta}_{1}},k=1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}n1$.

The restore value of ${\widehat{x}}^{(0)}(k)$ can be given by
Proof:
Substitute ${\beta}_{\text{1}},\hspace{0.17em}{\beta}_{\text{2}}$ into ${x}^{(1)}(k+\text{1})={\beta}_{\text{1}}{x}^{(\text{1})}(k)+{\beta}_{\text{2}}$, then:
Let ${x}^{(1)}(1)={x}^{(0)}(1),$ then
2.4 Forecasting Error
After being specified, each model’s performance should be verified or validated by comparison of its forecast with historical data for the process it was designed to forecast. This is no compromise among researchers as to which measure is best for determining the most appropriate forecasting method. Accuracy is the criterion that determines the best forecasting method; thus, accuracy is the most important concern in evaluating the quality of a forecast. Some of the common indicators used to evaluate accuracy are MAD (Mean Absolute Deviation), MSE (Mean squared error), RMSE (Root mean squared error) and MAPE (Mean absolute percentage error). The smaller accuracy indicators are, the better models perform.
3. CASE ANALYSIS
3.1 Data source
The cacao powder is not either the endproduct or the input of the production at PGPI. Cacao bean, which is the input, can be produced all time of the year. It is obvious to consider that forecast for the growing state is not suitable. In the other hand, one semifinish product like the cacao powder seems to be more promising. Cacao powder takes approximate 70% of the resources for production and R&D. Moreover, it is also an agricultural product, so it can take advantage of other available research on forecasting demand in agriculture.
Historical data comprising the period from January 2016 to December 2017 was collected to support the analysis of the demand forecast. Database was exported directly from the official daily sale report provided by the Financial and Logistics Departments.
The Table 1 and the Figure 3 below illustrate the descriptive data variable. Figure 3 was exported from XLSTAT addin of Microsoft Excel, shows the distribution of data and gives us a summary on the trend.
3.2 Research Process
Figure 4 describes how the study is conducted.
3.2.1 ARIMA:
Figure 5 illustrates the process of fitting the ARIMA model to our data.
The first step in fitting an ARIMA model is to recognize the value p, d, q. As we mentioned above, these three values decide whether our series require AR or MA terms to correct any autocorrelation. The ACF plot is a bar chart graphs the coefficients of correlations between a time series and lags of itself. On the other hand, the PACF stays true to its name that is a plot of the0partial correlation coefficients between the series and lags of itself.
The ACF and PACF plots below are produced by using XLSTAT addin of Microsoft Excel.
It is stated that by observing the ACF and PACF plot of a sequence, we can recognize the value for p and q, or how many of AR and MA0 terms are necessary to correct the autocorrelation in our series. Here in the PACF, we see the correlation at lag 1 is quite significant and positive, in addition, it has a sharper cutoff line than the ACF. Hence, we conclude that an AR term should be added to the model. Moreover, the PACF has one significant spike, therefore our series displays an AR(1) signal, or p=1.
Next, we fit the ARIMA (1,1,0) model with our data by XLSTAT addin, then obtain a ACF and PACF plots for the residuals as below:
The autocorrelation at the crucial lags 1,2 has been abolished, and the plots show no significant pattern in higherorder lags. Therefore, we accept the forecasting results of fitting the ARIMA (1,1,0) (Figures 6,7, 8 and 9).
3.2.2 Exponential Smoothing
The general equation for Exponential smoothing as below:
where ε_{t} is the difference between the actual and forecast value for period t (Table 2).

Choosing the smoothing weight α = 0.02

Assume the forecast error of period 1 is 0. We then start to calculate the forecasting value for period 2
3.2.3 GM(1,1)
For GM(1,1), we can only use the positive number for calculation. In order to not let the randomness affect the result, the first step is using AGO to calculate X^{(1)} from X^{(0)}
Apply to our set of data X^{(0)} = {1994.897; 8065.967; 7599.073; …; 15301.5409}
(n = 24)
Calculate X^{(1)} by performing the accumulated generating operation:
X^{(1)} = {1994.87; 20060.864; 27659.937; …; 312739} Determine the background values z^{(1)} :”
z^{(1)} (2) = ½ (1994.87 +2006.864) = 16027.88
z^{(1)} (3) = ½ (20060.864 + 27659.937) = 23860.4
z^{(1)} (4) = ½ (27659.937 +35476.443) = 31568.19
…
z^{(1)} (24) = ½ (297437.94 +312739.48) = 156369.7407
Then, a and b are found by subtituting the X^{(0)} sequence of value into the Grey differential equation:
Convert the linear equation into the form of a matrix Let
By using the least square method, we obtain the two parameters a and b
Then, conduct the whitening equation of the differential equation by using the value of a and b
The prediction model is found based on the equation below:
Replace the value of k in to the equation, we obtain the correspondent value of x^{(1)} to x^{(0)}
Derivethe predicted value of the original series according to the AGO and obtain the forecasting value for X^{(0)}(t). The values are presented in Table 3.
3.2.4 DGM (1,1)
The process for applying and calculation of the DGM(1,1) is describe as below

Collect a time series data sequence

Convert the raw series X^{(0)} into a series X^{(1)} using the AGO
where $\begin{array}{l}{X}^{(1)}=\{{x}^{(1)}(1),\hspace{0.17em}{x}^{(1)}(2),\hspace{0.17em}\mathrm{...},\hspace{0.17em}{x}^{(1)}(1)\},\hspace{0.17em}\hspace{0.17em}{x}^{(1)}(1)={x}^{(0)}(1)\\ {x}^{(1)}(k)={k}_{I=1}{x}^{(0)}(i),\hspace{0.17em}\hspace{0.17em}k=2,\hspace{0.17em}3,\hspace{0.17em}\mathrm{...},\hspace{0.17em}n\end{array}$

The firstorder and single variable DGM(1,1) equation is describe as below:
β_{1} and β_{2} are parameters

Using the least squares method to calculate parameters which should satisfy:
$${B}^{\text{T}}B\widehat{\beta}={B}^{\text{T}}Y,\hspace{0.17em}\hspace{0.17em}\widehat{\beta}={({B}^{\text{T}}B)}^{1}{B}^{\text{T}}Y$$where $Y=\left[\begin{array}{c}{x}^{(1)}(2)\\ {x}^{(1)}(3)\\ \vdots \\ {x}^{(1)}(n)\end{array}\right],B=\left[\begin{array}{cc}{x}^{(1)}(1)& 1\\ {x}^{(1)}(2)& 1\\ \vdots & \vdots \\ {x}^{(1)}(n1)& 1\end{array}\right]$
With B^{T} is the transpose matrix of B

Derive the simulative value of x^{(1)}(k+1) as
Thanks to the R studio, researchers can construct the code following the R studio tutorial then input data and fitting the models. Here the the parameters used for Cacao Powder Demand calculation are β_{1} =1.0344, β_{2} = 8213.8631 , therefore, ${x}^{(0)}(1)({\beta}_{1}1)+{\beta}_{2}=8626.6234$
3.2.5 Verhulst
As the step by step tutorial of identify the Verhulst model for a sequence of data is given in Section – 2.4.2. Verhulst model, however, it is still complicated to fitting the model. Again, by using the R studio application, research can compute the parameters easily with high accuracy level. Given the parameters of Verhulst in calculation, we have a = 0.0415 and b = 0.00, ax^{(1)} = 498.0000, bx^{(1)} = 0.0249 and a – bx^{(1)} (0) = 0.0166.
4. FINDINGS AND DISCUSSION
4.1 Result
Following the calculation tutorial that has been demonstrated in the section 3, the Table 3 below describes the detail result.
4.2 Evaluation
Table 4 summarizes the errors occur using MAPE, MSE, RMSE, and MAD. This is an evaluation on how well each model performed.
MAPE
With y1 is the model prediction value, yt’ is the actual value of observation

 if MAPE < 10%, model has an “excellent” forecasting ability

 if MAPE < 20%, model has “good” forecasting ability

 if MAPE < 50%, model has a “reasonable” forecasting ability

 if MAPE > 50%, model has a “poor” forecasting ability
MSE
MSE is the short form of Mean squared error. For calculating the MSE, we estimate the sum of squared of the errors, then multiply by 1/n with n is the number of observation.
RMSE:
Root Mean Square Error (RMSE) is the standard deviation of the residuals (prediction errors). It estimates how the distance between the regression line and data points is; RMSE is a measure of how spreads out these residuals are. In other words, it tells you how concentrated the data is around the line of best fit.
MAD
Mean absolute deviation (MAD) of a data set is the average distance between each data value and the mean. Mean absolute deviation is a technique to identify variation in a data set. Mean absolute deviation gives us an analysis on how spread out the values in a data set. The smaller/lower the value of MAD, the more accurate the forecast is.
4.3 Discussion
Significantly, DGM(1,1) receive an excellent MAPE at 9.19%. GM(1,1) and Verhulst are the two methods that also have a good MAPE, 10.83% and 12.45% respectively. Moreover, the three GFMs also have the low MSE, RMSE and MAD. ARIMA and Exponential smoothing however do not appear to work well since they obtain very high points regarding of the four indicators. In short, with given data and case, DGM(1,1) proves the ability to produce a best forecasting, comparing to other four approaches.
The real observation data is considered as a midlong data, spread along 2 years with corresponding 2 sequences. Therefore, any lack of relative information such as seasonal factors, market impact, consumer trend, weather impact, etc. may cause huge deviation. Despite of the fact that real data was fair and confirmed stationary, most of the model gave quite large error. In this situation, those complex models like ARIMA which was considerably accurate for full information cases, seems unusable with huge error. ARIMA figured out the trend of data but missed several of outliner, due to unexpected impacts.
As our main purpose is choose identify the best methods for the establish of demand forecast for PGPI, let use the three outstanding models including DGM(1,1), GM(1,1) and Verhulst to create an forecast for 2018.
Table 5 describes the future demand in 2018 with the use of our three chosen models
Then, for comparison, we have Figure 10 graphing the actual demand in 2017 and the forecasting demand in 2018.
We can see from the illustrator that DGM(1,1) and GM(1,1) produce quite similar forecasting result for the 2018. From the last section, we also concluded that DGM(1,1) and GM(1,1) are the two models that have the highest accuracy level, generating smaller errors comparing to other models. Nevertheless, three models all present a possible increase in the demand for cacao powder in 2018.
In the ending of 2017, the company observe a significantly decrease in sale, as this is claimed to be the general situation of cacao industry. According to the ICCO, within the period of 20162017, the cacao manufacturers suffers a huge loss as the price for cacao drop to the lowest value that has been seen within five year in a row. However, experts also believed that the surplus in production from 2017 would make the manufacturers become hesitate as they are aware of the instability of the market in recent years. Therefore, when it comes to production planning, it is expected that farmers and manufacturers would decrease the scale of production to avoid generating any surplus that can harm their business. Now, let take a look at the market within the first 5 months of 2018
The cacao price had been slowly increase, which demonstrate the recovery of the market. Started from the April of 2018, the price significantly decreased, however, it still not as low as the record from 2017. Compare to our forecasting result in Figure 10, we observe the similarity in trend of the forecast and actual values. As the price for cacao in the international market decreased in April 2018, the forecasting demand of cacao powder at PGPI also suffer a shrinkage. Base on this discussion and our evaluation from the last section, the models, especially DGM(1,1) and GM(1,1) has performed a good ability to predict future trend. In addition, both Table 1 and Figure 3 indicate a an upturn in the market. Therefore, manufacturers can expect to see a recover of cacao market in 2018. This would help the managers at PGPI as well as farmers and other firms within the cacao industry of Vietnam to schedule for the purchasing of resources and the production planning.
5. CONCLUSION AND RECOMMENDATIONS
Globalization has caused the manufacturing conditions that many firms operate under to change dramatically. It is thus very important for manufacturing firms to efficiently control their manufacturing systems. The early stages of a manufacturing system are especially critical for managers. Setting the foundation for manufacturing system are especially crucial to managers. If problems can be addressed at a premature stage and appropriate actions are taken to overcome them, the business will become more effective and efficient.
In order to be prepared for the future scenarios right in the first place, a company must adopt a good forecasting technique. Once a forecast is established, a manufacturer can schedule their production efficiently with resources are wellallocated. Regardless of how much a good forecast can do for SCM, its error also causes huge damage to a business. One fail prediction can lead to a sequence of wrong decisions. Therefore, selecting a method for prediction asks for a lot of effort and resource (Sayyadi and Awasthi, 2018b;Shukla and Jharkharia, 2013).
This research analyzed data taken from the main factory and office of PGPI located in Dong Nai Province, identified the suitable forecasting techniques for demand of cacao powder. The result after being evaluated and compared has shown that DGM(1,1) and GM(1,1) work efficiently on the data set and manage to provide a good prediction for production planning. These two methods also calculated the suggested demand for the upcoming period t as it is shown. Since we only include five forecasting methods for comparison in this study, there are still rooms for further research. DGM(1,1) and GM(1,1) might or might not be the final options for developing demand forecast. However, the managers of PGPI still can take into account if any of the forecasting result are seems to be not consistent with the current performance and resources distribution. In short, the result of this research can be considered as an indicator showing whether or not there is any bottleneck within the whole system. Also, it is encouraged for the company to invest in further research with the variety of models for a better result.
The lack of predication in concern of quantity of output and the demands of customers not only results in the low efficiency in the process of storage but also bring the negative effects in the figure of environment. To eliminate those minuses of forecasting, people should pay attention in the inventory management literature (Kazemi et al., 2018). Therefore, “it is required to jointly incorporate those two relevant aspects in a single research to support decisions, compare the results and obtain new insights for complexities in practice.” If researchers can deeply reach the exact order quantity as a sustainable point of view, the impact of emission costs and the problem of waste could be tackled immediately. Furthermore, those of them will give a hand to handle for the other fields such as solving the shortage and limited warehouse space (Hoseini Shekarabi et al., 2018).
Due to the limitation of the historical data and knowledge available, this research only touches an aspect within SCM for agriculture product. In particularly, yet this is not a tutorial for guiding the PGPI to choose one specific methods for demand prediction. However, I would consider the study as a reference to establish formal understanding on this matter because of its undeniable important role in the sustainable growth of the company (Singh, 2015;Tsao, 2015).