Introduction

Vanilla (Vanilla planifolia) is an aromatic orchid native to Mexico, whose species lives in tropical regions of America, although its main production area is located in the area known as Totonacapan, which covers 20 municipalities in the state of Veracruz and 19 in the state de Puebla (^{Damirón, 2004}; ^{OMPI, 2009}; Díaz et al. 2017). From the fruit benefited from this species, vanillin is obtained, the most popular and demanded flavoring and flavoring in the world, widely used in the soft drink, wine and culinary industry, placing this crop as the most profitable, only after saffron, with an estimated demand by specialists of the flavoring industry of around 14 000 t for the year 2012 (^{Lamas, 2012}).

However, despite its importance, only around 100 tons of natural vanilla extracts are produced worldwide, the rest corresponds to synthetic substitutes based on eugenol, coumarin and other aromatic fermentation by-products (^{Retes et al., 2015}). National vanilla production increased from 257 t in 2003 to 494.69 t in 2018 and registered an average annual growth rate (TCMA) of 4.78%, higher than that of the population (1.9%) (^{SIAP, 2017}). Even though Mexico is the center of origin and diversity of vanilla and has the ideal agroecosystems for its growth and development, together with the extensive technical and traditional knowledge, widely developed by the Totonaca and Nahuatl culture (Díaz et al., 2017).

The production of fruits in green obtained in 2017 was 514.94 tons, placing the country as the fourth world producer, after Madagascar, Indonesia and China, which in the same year produced 3 227, 2 402 and 662 t respectively (^{FAO, 2017}). Therefore, Madagascar is the country that dominates the world market in volume of production, making other producing countries have to adapt to the circumstances of that country. The main producing entities of this orchid in Mexico are Veracruz (76.08%), Puebla (11.55%), Oaxaca (10.5%), San Luis Potosí (1.87%) (SIAP, 2018).

In the three main states, there are 11 municipalities with the highest number of vanilla producers, among which are: Papantla, Cazones de Herrera, Gutierrez Zamora, Tecolutla, Tuxpan and Tihuatlan, all of them belonging to the Totonacapan Zone in the State of Veracruz , and by the State of Puebla, the municipality of Venustiano Carranza, Pantepec and Ayotoxco de Guerreo, also highlights the production of the municipalities of San Juan Bautista Valle Nacional and San Miguel Soyaltepec in Oaxaca (SIAP, 2018).

The destination of the vanilla benefited is the international market, approximately 90%, the remaining 10% is distributed in the national territory in the handicraft market and in the extract industries (^{CESPV, 2007}). There are important factors that influence the decision-making of the members of the agricultural sector, such as the predictions and production estimates established by government agencies, so they must be efficient and not biased, ^{Sanders and Manfredo (2002)}. For this reason it is convenient, both for producers and for designers of agricultural public policy, to have means that help to know the behavior of the market and that allow access to updated information that makes it possible to predict the direction of the market (^{Myers et al., 2010}).

Given the variability of the vanilla offer and its level of growth in a production is influenced by biological and climatic factors, ^{Allen (1994)}, it is relevant to know in depth the behavior of this offer. For this purpose, the selection and estimation of a statistical model with adequate predictive capacity that allows representing the information generating process (PGI) is essential. Forecasting methods are used to make statements about the future value of a variable under study ‘since having indications of what will happen allows reducing uncertainty’ (^{Montemayor, 2013}).

According to ^{Evans (2003)}, since the sixties many economic analyzes have been carried out with the help of multivariate linear regression models, which seek to obtain forecasts although in many cases with particularly disappointing results, because it is based on a causal relationship established by theory or logic and it is assumed that the future values of the variable under study can be determined by projecting the values of influence variables that may or may not be controlled by the researcher or in some cases it is impossible to have the information of explanatory variables of importance (^{Montemayor, 2013}), so to achieve better predictions, most economists and practical forecasters have resorted to the use of statistical methods that do not involve economic theory such as univariate time series models (^{Evans, 2003}), but with the disadvantage that the construction of these is based on the use of behavior of the same variable to make the prediction (^{Asteriou and Hall, 2007}).

Under this criterion ^{Luis et al. (2019)} he used a univariate time series of monthly white egg prices paid to the producer from 2000 to 2016 and through the Box Jenkins methodology he demonstrated that the series behaves as a seasonal model Seasonal Autoregressive Integrate Mobile Averange (SARIMA) of coefficients (0,1,1)X(1,0,1)s= 12 and with the results of the model they were able to make estimates regarding the price paid to the producer, until December 2019.

However, unlike the economic variables necessary for the construction of econometric models, in the study of the behavior of agricultural productions over time, the theoretical foundation of the model is not indispensable, which increases the applicability of a time series model (^{Barreras et al., 2014}). In the present investigation the methodology of integrated autoregressive processes of moving averages (ARIMA) was used. These models are more appropriate for long-term predictions than for short-term ones that present seasonal patterns; however, these are analyzed on a stationary series and at least 50 data are needed (^{Box and Jenkins, 1976}), hence their prediction is more reliable.

The ARIMA model is useful in situations where time series data show periodic seasonal fluctuations that are repeated with almost the same intensity each year (^{Box and Jenkins, 1976}). This characteristic makes this model suitable for studies of agricultural products and therefore for the national production of vanilla (^{García et al., 2003}) stated that the biological process of production and the climatic factors, cause that most of the production of most agricultural products goes out in certain months of the year, causing a marked seasonality of production and, consequently, of prices.

In the case of vanilla in Mexico, the harvest is marked by the Law of Promotion and Protection of Vanilla issued on September 13, 1941 and in force to date, since it establishes that the harvest is made from November 15 each year and this does not last more than a month (^{Hernández, 2014}). The objective of this study was to develop a time series model to forecast the production of green vanilla in Mexico (PVAINI) based on the production of green vanilla available between 1961 and 2016. With the data collected from 2015 and 2016, the model was validated. Given the uncertainty of the market, the PVAINI forecast provides useful information to facilitate planning in the decision-making of the different actors that are part of this sector. The hypothesis was that the production of green vanilla in Mexico in the future is explained by the production that occurred 4 years earlier.

Materials and methods

To know the behavior of PVAINI and make forecasts, a historical series of green vanilla production was used, provided by the Agrifood and Fisheries Information System (SIAP, 2019) and by the database of the United Nations Organization for the Food and Agriculture (FAO, 2019). The PVAINI time series was divided into two parts: data from 1961 to 2014 expressed in metric tons (t), which were used to develop the time series model; also, to validate the model, the 2015 and 2016 production data were used.

If it is assumed that Y'=(Y1,Y2, …,Yn) is a time series, then an ARIMA model is denoted as (p, d, q) (^{Guerrero, 1991}) and is determined by:

 ϕpB ∇dYt=θqB at

Where: B = is the delay operator (ΒrYt=Yt-r) ; Yt = is the time series of the variable Y at time t;  ϕpB = is the autoregressive polynomial (AR), of order ‘p’, it develops like this,  ϕpB= 1-ϕ1B-ϕ2Β2-…-ϕpΒpθqB = is the polynomial of moving averages (MA), of order ‘q’, where: θqB= 1-θ1B-θ2Β2-…-θqΒq ; ∇d=Yt-Yt-rD=1 it is the number of differences needed to park the series and is symbolized ‘d’; at = refers to the terms of random error or also called white noise that are assumed to be random variables independently distributed identically sampled from a distribution preferably with a mean equal to zero and variance at∼ N0 ,δ2 ( ^{Box and Jenkins, 1976}).

In the analysis of time series, according to ^{Box and Jenkins (1976)}, the variable at is commonly known as white noise and is interpreted as an exogenous effect that the model does not explain. Taking into account the PVAINI time series, this white noise can be an effect of climatic variables, diseases, the peso-dollar exchange rate, the country risk situation, natural disasters, etc. The PROC ARIMA procedure of the statistical analysis system (SAS) software version 9.4 was used. ^{SAS (2014)} , to estimate the ARIMA model to the PVAINI series of the 1961-2014 period using the Box and Jenkins methodology for the construction and adjustment of the forecast model (identification, parameter estimation, assumption verification and prediction) (^{Box and Jenkins, 1976}).

The suitability of each model was verified by graphs of the residuals against time. ^{Ljung and Box (1978)}, if the model is adequate, the residuals would be expected to be distributed randomly around zero without showing any pattern or trend; in case there is some kind of irregular behavior, this would be a reason to assume that the residuals do not have a normal distribution and are not purely random. In this work, we used the parsimoniously model choice suggested by ^{Burnham and Anderson (2002)} using the akaike information criterion (AIC) and the Schwartz bayesian criterion (SBC) which compare the goodness of fit of the different models.

Both criteria are based on the use of the sum of the squares of the errors, seeking to minimize it from various combinations of p and q. Lower AIC and SBC values indicate a better fit in the model.

Results and discussion

The average annual growth rate (TCMA) of the PVAINI series for the period 1961-2017 was 1.43%, which indicates that it has a slight upward trend, however; it is also observed that for the period 1961-1979, vanilla production declined to a TCMA of -16.9%, said production behavior was due to the fact that less hectares of land were allocated to this crop for those times, when passing from land agricultural producers of this orchid to land for livestock production, for oil production, extreme weather and discouragement in the international price of vanilla, which caused losses in vanilla production.

According to ^{Damiron (2004)}, in 1980 it was estimated that there were at least 800 organized and 700 unregistered producers that produced 30 tons of vanilla, in this context, it would be expected that the production of this orchid will stagnate; however, for the 2012-2017 period, a positive growth rate (5.72%) is presented, so production would be expected to increase in the short term.

The highest peaks of the time series correspond to the years of 1993, 1998 and 2007 and are explained by good weather conditions and the speculation of the price that occurred three years earlier, while the lowest production of vanilla in Mexico corresponds to the period 1978-1980, explained in the previous paragraph, regarding the fall in production in 2002 are explained by extreme droughts that occurred in 1996 and 1997 and the tropical depression number 11 occurred in October 1999, while the fall in production in 2011, was heavily influenced by excess rainfall (hurricanes Stan in 2005 and Dean in 2007) that affected the Totonacapan area and damaged a large part of the crops (Figure 2a).

Figure 2 Graphic diagnostics of the original series (PVAINI): (a) Original behavior of the PVAINI series in (t), 1961-2017; (b) Behavior of the PVAINI series transformed to natural logarithms (LPVAINI); and (c) LPVAINI series differentiated, transformed with natural logarithms and without apparent tendency. Elaboration with data from ^{SIAP (2017)}; ^{FAO (2017)}. 

In the Box-Jenkins methodology, it is necessary that the series studied be stationary in mean and variance, the first one is obtained with at least one difference, while to transform the variance a transformation of the Cox-Box family is applied, these in addition to stabilize the variance, they improve the approximation to the normal distribution of the process Zt , where λ is the transformation parameter. The Cox-Box test showed a λ= 0.5, so the PVAINI series had to be transformed into natural logarithms to have constant variance.

The magnitude of this coefficient is consistent with the theory of time series, proposed by (^{Cox and Box, 1964}) for the application of any of the ARMA family models. Then Y'= (Y₁, Y₂, ..., Y_n) is the vector of natural logarithms, now renamed the series as LPVAINI. It is observed that the series Y₁, Y₂, ..., Y_n against time there is still a certain tendency, but through the first difference (𝛻) a stationary series is obtained. Therefore, d= 1 is considered (Figure 2b).

The PVAINI series differentiated and transformed to natural logarithms has been renamed LPVAINI (1), and it is a priori intuited that it is already stationary on average (there is no trend) and in variance, a situation that agrees with what is established with ^{Box and Jenkins (1976)} since the series to study must have a constant variability over time and without tendency because in the methodology worked it is a necessary condition to obtain better forecasts (Figure 2c). The stationarity or non-stationarity of the time series was statistically verified, proposed by ^{Dickey and Fuller (1981)} by means of the augmented dickey fuller unit root test (ADF) which consists of including lags of the first difference of Yt in the test regression to admit the possible existence of serial autocorrelation.

For this, the following hypothesis is proposed, Ho: the series is NOT stationary (ρ= 1) and has a unit root vs Ha: the series is stationary (ρ≠1) and has no unit root. Decision rule: Ho is rejected if │tau calculated│≥│tau of tables│ (^{Pankratz 1983}). Since the absolute value of the calculated tau (7.81) is greater than the tau of tables, the null hypothesis is rejected and it is concluded that the PVAINI series does not have a unit root, therefore, it is stationary; that is, it has a mean and constant variance over time (Table 1).

The estimated autocorrelation function (ACF) and the partial autocorrelation function (PACF) of the ∇ LPVAINI series show a high ACF in the lag = 13. With respect to the FACP there are high peaks in the correlations for the lags = 3 and 4, so it is intuited that an ARIMA process describes and explains the behavior of PVAINI (Figure 3).

Figure 3 Autocorrelation function ACF (a) and partial autocorrelation function PACF (b), estimated for the series Yt=1-BLPVAINI . Elaboration with the output results of PROC ARIMA, ^{SAS (2014)}. 

The stage of identification of the model is a visual inspection of the ACF and the PACF plotted although they can hide information of interest for the choice of the appropriate model, according to what was proposed by ^{Gujarati and Porter (2009)}, for which 10 models were proposed. Of these, the four best ones that meet the significance of parameters and white noise were identified. AR, MA coefficients were chosen since the ACF and PACF simultaneously present non-zero coefficients, also, they were estimated considering different choices of p and q, as well as the values of SBC, AIC and the variance δε2 for the four best models ARIMA adjusted to the LPVAINI series.

It is considered that the first difference in the non-seasonal part is usually sufficient (d= 1) d was set to 1 in all models. The model with the lowest SBC and AIC value for this data set was the ARIMA (1, 1, 1) (Table 2).

Using the PROC ARIMA process (^{SAS 2014}), the parameters AR1,1 ( ϕ4 ) and the moving average component MA1,1 ( θ13 ) were calculated, using maximum likelihood, because this method starts with its estimators being asymptotically optimal and if the size of the series is large, it can be considered that they are centered or unbiased, efficient and their distribution is normal (^{Montemayor, 2013}). For the interpretation and significance of the estimators (^{Box et al., 1994}), they state that for a moderate model to be considered, those parameters whose absolute t-statistic is greater than 2 or the p-values must be less than 0.05 must be included. Therefore, the ARIMA model (1, 1, 1) is considered moderate, since it fits sufficiently with the old data without using any unnecessary parameters (Table 3).

After estimating the model parameters, it was validated through the analysis of residues (^{Yafee and McGee, 1999}). Standardized residues were analyzed through a histogram, the respective ACF graph and the p-values for the white noise tests. It is clarified that the estimated standardized residuals, from this model should behave as an independent and identically distributed sequence, with a mean of zero and a constant variance (Figure 4a). The distribution of the waste approached a normal distribution, indicating a high affinity of the data, regardless of the magnitude of the data (Figure 4b).

Figure 4 Graphic diagnostics to evaluate the adjustment of the ARIMA model (1, 1, 1): a) standardized waste; b) histogram of standardized residues; c) ACF of the residuals; and (d) P values for white noise tests. Elaboration with the results of the best estimate of the ARIMA model (1, 1, 1). 

The ACF of the residues suggests that autocorrelations are within the confidence band; that is, they are close to zero. This result indicates that the residues did not deviate significantly from a zero white noise process, and are purely random, so there is no longer any dependence information of some data with others over time (Figure 4c). Given the high p-value, associated with Chi-square statistics, there is no reason to reject the null hypothesis that the residuals are white noise (p value greater than or equal to 0.05). Therefore, the ARIMA model (1, 1, 1) conforms to the behavior of the LPVAINI data (Figure 4d).

The predicted out-of-sample values for the years 2017 to 2024 taking into account the ARIMA model (1, 1, 1) were compared with those of the PVAINI series, the predictions are good, with an average absolute percentage error (MAPE) of 8.32 which indicates that our forecast is wrong by ±8.32%. However, the forecast should not be based solely on this indicator, the Chi-square test should also be examined to see if the series presents residuals that behave like white noise.

This model predicts PVAINI values that are very close to those observed, because they are all within the estimated confidence band of ± 95%, under this criterion that, if urgent measures are not taken by the actors that influence the public policies, vanilla production in Mexico will decline rapidly in the next seis years to a TCMA of -5.05%; that is to say, the estimated production for 2019 will be from 464 to 340 t in 2024. Likewise, the 95% confidence interval of the values outside the series was estimated and it was found that by 2019, the range of occurrence will be of 188.3 t minimum and 1 143 t maximum, while by 2024, with the same confidence interval, the range is at 40.31 t minimum and 2 868.2 t maximum (Figure 5).

Figure 5 Observed and predicted vanilla production in Mexico, 2017 (metric tons), obtained through the ARIMA model (1, 1, 1). Prepared by the authors with the results of the ARIMA model (1, 1, 1) of the ^{SIAP (2017)} and ^{FAO (2017)}. 

This model is used to predict outside the sample the PVAINI and was able to do it with relative precision for a later year. However, ^{Chu (1978)} argues that these estimates may not be credible for predicting prices in the medium and long term, so it should be noted, however, that the main drawback of this type of model is its own foundation, since being based on past events and learning from the history itself to make the predictions, these will be correct to the extent that the factors that determine evolution continue to act in the same way and are not noticeably altered.

Since this product is important for its high commercial value and appreciation in gourmet haute cuisine, in addition to the fact that its cultivation requires edaphoclimatic conditions typical of certain tropical regions of Mexico, it is important that the actors that influence the design of agri-food policies, urgently propose the rescue, promotion and increase of the area sown with vanilla, as well as the impulse and enforcement of laws that punish theft of the fruit in plantation, intermediary and speculation since if not, not only production will decline but also the spirit to continue sowing vanilla.

Conclusions

It is possible to model the PVAINI in Mexico using the methodology proposed by Box and Jenkins, since the behavior of these series is not completely random, and they are described as time series with a high probability of success in modeling through this methodology. The best model that represented the behavior of the PVAINI was the ARIMA (1, 1, 1) by demonstrating that the production of vanilla in green in a given year is estimated by the production of vanilla that occurred 4 years earlier. Parameter estimates and forecasts were statistically appropriate and significant, even though the series studied shows high variability.

The predictions of the model in the short term differ by ±8.32% of the observed data, minimizing the random error, finding that in the short term the vanilla production in Mexico will decline at a rate of 5.05% per year. The estimated PVAINIs provide useful information to plan, encourage and rescue in an urgent way the cultivation of vanilla in Mexico, through the design and implementation of public policies of an agri-food nature, since it is a product originating in our country and profitable to improve the income of the producers, because otherwise, vanilla production in Mexico in the medium term tends to decrease rapidly.

Recommendations

It is suggested that for future studies where the Box-Jenkins methodology is used, new variables that explain the random disturbance and that are directly associated with the increase in PVAINI are introduced into the model, although their relationship is not instantaneous or static, consider the temporal dimension of the observations, allowing to measure and evaluate the dynamic responses and the way in which the effects are transmitted between the variables (Box and Jenkins, 1976). The variables to be considered would be the annual average prices of artificial vanilla flavors and the annual production of the countries of Madagascar and Indonesia, which would require more precise estimates, since by including unidirectional causality variables (independent variables) since the entry towards the response variable (dependent) the MAPE of the forecast should decrease (^{Chigona et al., 2007}). For this, the future use of multivariate time series models also called transfer function models (TFM) is proposed, unlike the univariate ARIMA model (1,1,1) estimated in this work for the PVAINI series. According to ^{Keller (1987)}, these models allow us to simulate and study how certain ‘scenarios’ defined by possible evolutions of the explanatory variable affect the response variable, hence the prediction is more reliable than the methods that use the classical regression model.