Universidade Federal de Santa Maria

Ci. e Nat., Santa Maria, v. 43, e83, 2021

DOI: 10.5902/2179460X41221

ISSN 2179460X

Submitted: 19/11/2020 • Approved: 11/09/2021 • Published: 15/10/2021

Sessão

Econometric modelling of time series applied in the generation of subsidies in the milk production chain in Rio Grande do Sul

Modelagem econométrica da série tempo aplicada na geração de subsídios na cadeia de produção de leite do Rio Grande do Sul

Francisca Mendonça SouzaI

Brunna Aver de PaulaII

Daniel Knebel BaggioIII

Argemiro Luis BrumIV

Juliana da Fonseca Capssa Lima SausenV

I Doutora em Estatística e Análise de Dados, Professora, Investigadora Associada UNIDE-(IUL), Regional University of Northwestern Rio Grande do Sul, Ijuí, RS, Brazil

http://orcid.org/0000-0003-3011-0086 - kikamatcom@gmail.com

II Pós graduada em Finanças e Mercado de Capitais, Regional University of Northwestern Rio Grande do Sul, Ijuí, RS, Brazil

http://orcid.org/0000-0002-5720-1896 - brunnaaver@hotmail.com

III Mestre em Contabilidade e Finanças, Doutor em Contabilidade e Finanças, Regional University of Northwestern Rio Grande do Sul, Ijuí, RS, Brazil

http://orcid.org/0000-0002-6167-2682 - danibaggio@gmail.com

IV Doutor em Economia Internacional, Professor, Regional University of Northwestern Rio Grande do Sul, Ijuí, RS, Brazil

http://orcid.org/0000-0002-8763-9514 - argelbrum@unijui.edu.br

V Doutoranda em Desenvolvimento Regional, Mestre em Desenvolvimento Regional, Regional University of Northwestern Rio Grande do Sul, Ijuí, RS, Brazil

http://orcid.org/0000-0003-4146-8294 - jucapssa@gmail.com

ABSTRACT

The purpose of this article is to make short-term forecasting using the methodology Box & Jenkins, the Johansen method and the Granger causality, and the impulse-response function between variables price and milk production in Rio Grande do Sul’s market. The monthly price series of milk and its production in Rio Grande do Sul were analysed, in the period from January 1995 to December 2017. The model that suits best, for forecasts, the data of the series of the milk price paid to the producer was an ARIMA (1, 1, 1) and production was SARIMA (1, 1, 1) (1, 1, 1)6, which provided reasonable estimates of forecasts for the months from February to July of 2017. The use of Johansen methodologies identifies the existence of the one cointegration vector and a long-term equilibrium relation between variables price and production of the milk. When we analyse Granger’s causality, the results point to a two-way relationship, that is, prices influence milk production and vice versa. The analysis of the impulse-response function showed that the shocks present significant impacts between production and cost, both in terms of duration and intensity.

Keywords: Econometric Modelling; Milk Production; Cointegration vector

1 INTRODUction

The present study aims to analyze the production and the price paid to the bovine milk producer in Rio Grande do Sul, using temporal methodologies, looking for forecast models for milk production and the price paid to the cattleman. It also proceeds with the estimation of linear models of class Autoregressive Integrated Moving Average (ARIMA) with particular emphasis on the short term forecast. In this way, to diagnose and analyse the factors that determine supply and demand in this sector, highlighting the importance of milk in the agribusiness universe. In addition to the prediction that plays a crucial role in the planning of milk producers, this research analyses the degree of causality and cointegration between the price of milk paid to the producer and production.

The production of milk and its derivatives plays a relevant role in the supply of food and in generating employment and income for the population. Milk is among the six most important products of Brazilian agriculture. The agricultural sector generated 46,049 new jobs in May 2017, an increase of 2.95% compared to April. Since the beginning of 2017, this sector has accumulated 77 thousand jobs, the most favourable performance among the main areas of the Brazilian economy (EMBRAPA, 2017).

According to the Brazilian Institute of Geography and Statistics (IBGE), in the first quarter of 2017, agribusiness accounted for 23% of the Brazilian GDP. It grew 15.2% in comparison with the first quarter of 2016 and 13.4% comparing to the fourth quarter of the same year. This increase was mainly driven by the result of the agricultural sector, which grew 13.4% between January and March 2017, representing the most significant expansion in 20 years (IBGE, 2017).

Milk is responsible for 40% of jobs in the Brazilian rural environment (EMBRAPA, 2017). To have a more objective idea of the impact of this sector on the economy, the increase in final demand for dairy products by one million reais generates 195 permanent jobs. This impact surpasses those of traditionally important sectors such as automobile, construction, steel and textiles (EMBRAPA, 2017).

Brazil is the sixth largest milk producer in the world and grows at an annual rate of 4%, higher than in all the countries that occupy the first places. It accounts for 66% of the total volume of milk produced in Mercosur countries. The Gross Value of Production (VBP) of some segments of cattle production estimated for 2017 is R $ 536 billion, 4.3% higher than the previous year, which was R $ 514 billion, as reported by the Policy Secretariat (SPA), of the Ministry of Agriculture, Livestock and Supply (MAPA, 2017). From this estimate, livestock farming is responsible for R $ 168.4 billion. According to the UN Food and Agriculture Organization of the United Nations (FAO, 2016), over the past three decades, world milk production has increased by more than 50% to reach 769 million tonnes by 2013. By 2017, according to the Organization for Economic Co-operation and Development (OECD, 2018), world milk production reached 823 million tonnes, from a worldwide herd of 719 million dairy cows. By 2016, around 150 million households worldwide were involved in dairy production, with most developing countries characteristic of small-scale farmers as it provides rapid returns to small producers.

According to a survey by IBGE (2017) for the year 2014, milk production in Brazil reached 35.17 billion litres, a volume that continues in 2017, the South region being the largest producer in the country, representing 34.7 % of the national total. Minas Gerais remained the most significant Brazilian producer, followed by Paraná and Rio Grande do Sul. Besides, the country has the second largest cattle herd in the world, behind only India.

Given this context, it is essential to building a trend for the behaviour of the milk market, a fact that may significantly help the family farms, which are the vast majority in the production of liquid in Brazil in general and in Rio Grande do Sul in particular.

In this sense, the present article relies on the prediction methodology known as Box & Jenkins. These mathematical models capture the behaviour of serial correlation and autocorrelation between time series values, providing accurate predictions. According to Souza (2016), this refers to the systematic method of identification, adjustment and checking, that is, the integrated autoregressive models of ARIMA moving averages. To analyse the behaviour and relations of these variables, the econometric time-series analysis methods is used, based on autoregressive vector models (VAR), Johansen cointegration test, causality test and the impulse response function. Recent studies like Camelo et al. (2017), Machado et al. (2018) and Leal Júnior et al. (2018) were used to estimate the climate, temperature, energy and wind of a given region. We highlight the following previous studies who used this technique to forecast the price of other commodities, such as Dalzotto Artuzo et al. (2017), Kovačević et al. (2017), Pinheiro and Senna (2017) and Retes-Mantilla, Torres-Mancera and Hernández (2017).

The article is structured in five sections: introduction; literature review; the methodology used; presentation and discussion of results; and finally the main conclusions.

2 Literature Revision

In this section we will present: the evolution of the milk market in the world, Brazil and Rio Grande do Sul; aspects of the dairy chain in Rio Grande do Sul; the Box and Jenkins methodology; autoregressive vector models; the Johansen cointegration test; the causality test; and the impulse-response function, for a better understanding of this study.

2.1 Milk in the world, Brazil and Rio Grande do Sul

World milk production reached a total of 823 million tons in 2017, projected to reach 993 million tons in 2027 (OECD, 2018). The five largest world dairy producers in 2017, in order, were: India (161 million tonnes); European Union (155 million); United States (98 million); China (41 million); and Russia (30 million tonnes). Brazil came in sixth place with 29 million tons (ANUALPEC, 2017). The projection for 2027, according to the OECD (2018) is for Brazil to reach 36 million tons of milk.

Therefore, world milk production is expected to increase by 20.6% between 2017 and 2027, while Brazilian production would grow at a faster rate, reaching 24.1% in the period.

According to OECD / FAO studies (2018), developed and developing countries will produce respectively 10% and 34% more milk by 2027. The contribution of developed countries to world milk production will fall from 49% to 44% by the end of the decade studied. A drop of 0.2% is expected in the herd of developed countries, and productivity growth of 1.1% per year in the medium term. In developing countries, the increase in production will be due to the growth of 1.1% and 1.6%, respectively, of the dairy herd and productivity. Despite the expected improvement, yield will remain very low in many developing countries, and productivity growth will remain flat in absolute terms. In the medium term, however, milk production will increase further thanks to improved productivity.

During the reference period, the ranking of the five largest milk producers will change, with the European Union coming in first place, accounting for 20% of world production, followed by India (19%), the United States (12%), and China and Pakistan, both with 5%. Asia is the continent on which evolution will be most marked until 2027, with Pakistan and India accounting for most of the increase in production. At the end of the decade under review, it is in India that the progression will be more pronounced, passing gradually to the first place in the world, reaching 23% of the total world production in 2027. In progress, it will be followed by Pakistan, which will register a high average of 3.4% per year and will advance 6% on total world production by 2027. In both cases, production will be mostly consumed in the market in the form of fresh milk. In the European Union, as in the United States, the annual growth rate will be more modest than in the previous ten years, with a growth of 18% and 11%, respectively, less than the 20% and 12% recorded in the prior decade. However, they will continue to be leading players in the export of industrialised dairy products.

In the European Union, milk production is expected to grow at a rate of 0.8% a year over the ten years considered, i.e. a more modest percentage than the 1.2% observed in the previous decade, even with the end of quotas milk production in 2015. This is explained by the drop in the domestic demand for fresh dairy products, more precisely the low consumption of milk, which will compensate for the high local demand for cheese and exports. Milk production will increase as a result of increased productivity of 1.2% per year over the next decade, while the dairy herd will suffer a slight contraction (0.4% per year). The European Union will account for an essential part of the production of skimmed milk in the period studied, from 34% to 37%. On the other hand, cheese production will decline from 44% to 42%, while butter production will fall from 21% to 20%, and there will be no change in the production of whole milk powder, which will remain 14%. Dairy production will grow at a rate of 1.3% per year over the next decade. On the other hand, the growth rate of skimmed milk powder and cheese will fall to 3.8% and 0.9% per year, respectively while, the growth in whole powder milk and butter production, 2.4% and 1.5% per year, respectively.

US milk production is expected to increase by 1.1% a year over the next decade. This increase will come mainly from animal productivity (1% per year). Compared to the previous ten years, production growth will be less intense, at 2.1% per year for skimmed milk powder, 0.4% for whole milk powder, 1.9% for cheese and 2% for whole milk powder, 4% for butter. The percentage of United States production in world production of butter and cheese will increase slightly.

In turn, China’s milk output growth will be 1.8% per year, but its share of world production will remain unchanged (5%) by 2027. Much of the production will be consumed with fresh produce, which will an increase of 2.1% per year. China will continue to be a major importer of dairy products, with its imports increased over the next ten years, albeit at a more modest pace.

On the other hand, in Latin America and the Caribbean, milk production, compared to the reference period, will increase by 18%. However, the percentage of world production will continue at 10%. In Argentina, in the medium term, production is expected to increase 1.3% per year, with the resumption of the sector after the climatic difficulties of 2016/17 and 2017/18. In the same way, in Brazil, production is expected to increase by 2.2% per year, while the industrialisation of dairy products is expected to grow by 2.1% per year.

In Oceania, although New Zealand is not the largest producer of dairy products, the country is undoubtedly the leading exporter. Its production growth rate is expected to fall to 2.3% per year, down from 4% over the previous decade. This country is one of the leading producers and exporters of whole milk powder. It will account for 24% of world production and 53% of exports in 2027. In the ten reference years (up to 2027), the gains will come mostly from a further expansion of the dairy herd (1.3% per year) and productivity (1% per year).

Finally, according to the OECD / FAO (2018), in developed countries milk production will mainly be transformed into butter, cheese, skimmed milk powder and whole milk powder. The developed countries in the reference period will produce 87% of the skimmed milk powder, 79% of the cheese, 43% of the butter, and 46% of the whole milk powder. The share of developed countries in the production of butter and whole milk powder will decline slightly by 2027. The production of dry milk extract will increase by 10% in these countries, of which 37% will be used for cheese production, while 23% will go to the manufacture of skimmed milk powder, 20% for butter production, 10.5% will be processed in whole milk powder and 8.5% of other fresh milk products. In developing countries, in 2027, of the 33% of milk production, 85% will go to produce fresh products, 7% for butter, 4% for whole milk powder, 3% for cheeses and only 0.6 % for skimmed milk powder.

In this context, Brazilian milk production reaches 2017 with a total of 35.1 billion litres. The same level of 2014, although it has quadrupled between 1974 and 2014. This stagnation is due to the strong economic recession experienced by Brazil between the middle of 2014 and the end of 2016. The country is still an importer of milk and derivatives, having reached a total of 169 thousand tons of dairy products in 2017, of which 103 thousand tons were made of milk powder and another thousand tons of UHT milk. The largest suppliers of milk and dairy products to Brazil are the Mercosur neighbours. In 2017, Brazil bought US $ 264 million of these products from Argentina, which represented one-third of the country’s total dairy exports. Uruguay purchased US $ 236 million, representing 40% of total Uruguayan dairy sales in the same year. In turn, Rio Grande do Sul, in 2017, produced 4.6 billion litres of milk, ranking third in the country, behind Minas Gerais (8.8 billion) and Paraná (4.8 billion). This year, Rio Grande do Sul registered 1.4 million milked cows. In 2017, the industrial sector of the state of Rio Grande do Sul purchased a total of 3.4 billion litres of raw milk for processing. In this regard, Rio Grande do Sul is the second Brazilian state, behind only Minas Gerais with 6 billion litres purchased by the industry. Among the 14 largest Brazilian dairy industries, Rio Grande do Sul has one in seventh place. It is the Cooperativa Central Gaúcha Ltda (CCGL), established in the municipality of Cruz Alta, which in 2017 received 439 million litres of milk, an increase of 23.2% over the total collected in 2016. In 2017, The central cooperative worked with 4,302 milk producers, who provided 278 litres/ day of milk per day (an increase of 35.3% over the average of 2016). The average productivity of Rio Grande do Sul per milked cow, in 2017, reached 3,240 litres/ year, constituting the highest average productivity in the country (EMBRAPA, 2018).

2.2 Box & Jenkins Models

Box & Jenkins (1970) mathematical models can capture the serial correlation or autocorrelation behaviour among time series values, thus making future predictions (Souza, 2016). The modelling process consists of steps of model identification, parameter estimation and diagnostic analysis. Its objective is to find models that adequately describe data with fewer parameters, and the accuracy of forecasts obtained, thus obeying the principle of parsimony (Box & Jennkins,1970; Tibulo 2013; Souza 2016). The strategy for constructing the models according to Box & Jenkins methodology shall follow the iterative cycle (Figure 1).

Figure 1 – Box and Jenkins iterative cycle flowchart

Source: Souza (2006)

The flowchart refers to the systematic method of model identification, parameter estimation, model adequacy and forecasting. According to Marchezan and Souza (2010), in a first step, the process of filter identification and stationarity analysis of the series under study is carried out. In the estimation of the parameters, it is analysed whether the stationary series is generated by an autoregressive (AR), moving average (MA) and/or mixed (ARMA) process. Next, the identified and estimated model is checked to see if it is suitable for forecasting. This verification is carried out by analysing the waste and evaluating the order of the model. If the model is inadequate, the step of the estimation of new parameters is returned. After selecting the most suitable model comes the last part of the methodology Box & Jenkins, the realisation of forecasts (Souza, 2016).

The ARIMA models are the result of a combination of three filters: the autoregressive component, the integration filter (I) and the moving averages component, represented by the letters p, d and q, respectively (FAVA, 2000). For Souza (2016), the process of identifying the model consists of determining which of the AR, MA, ARMA or ARIMA filters make up the series. According to Box, Jenkins & Reinsel (1994), a process is stationary when it oscillates around a constant mean, with variance also constant (AKA homoscedastic) and defined as stochastic once it is controlled by probabilistic laws. ARIMA models capture and model short memory processes. Such feature is observed through the auto correlation function (FAC), which rapidly decreases to zero. A process {xt,t E Z} follows an ARIMA process (p, d, q) if Δdxt follows an ARMA model (p, q), that is:

(1)

where B is a retroactive operator, or a lag operator, d is positive integer and represents the difference, according to the order of integration; φ and θ are the parameters of the autoregressive and moving average order processes p and q, respectively and et~RB(0,σ2). Another factor to be observed in the modelling is the periodicity which data are presented. It occurs within a period of one year and it is called seasonality; it is represented by an S repetition period. Periodicity can be incorporated into the model by using a seasonal ARIMA model, as represented in equation 2:

(2)

The seasonal component (St) presents oscillations, occurred in subperiods such as, for instance, the years. They may occur monthly, quarterly, or even daily.

In SARIMA models, not only non-seasonal components (p, d, q), but also seasonal parameters (P, D, and Q) are added. These models are widely employed when the series suffers climate-related impacts or some specific annual events (WANKE; BARBOSA, 2010). Ferraz, Sáfadi and Lage (2009) define SARIMA(p,d,q)(P,D,Q)St as models that contemplate series which present seasonal autocorrelation.

2.3 Unit root test

The investigation of stationarity in time series modelling is essential to avoid spurious results and to define the appropriate model to be used. The basic idea of stationarity claims that laws of probability, acting through the process, do not change over time. In actual data such characteristic is rarely observed, since those series usually present some kind of linear trend, whether positive or negative.

Souza (2016) defines a process as strictly stationary when the probability density function, the set of observations Z1, Z2, ..., ZN, does not change with the origin of the time axis. However, the process remains in “statistical equilibrium” “. Whenever series present time-dependent mean and variance, it is not stationary.

The most used unit root tests are Dickey-Fuller Increased (ADF) and Kwiatkowski-Phillips-Schmidt-Shin (KPSS); they both validate ADF test (1984). Such test has as a null hypothesis that the series is non-stationary or integrated of order d (d>0). In most cases of series that are not stationary, when applying the first difference it will become stationary, that is, I (0). ADF test is based on the following regression for the variable Zt:

(3)

As a form of validation of the ADF test the KPSS (1992) test is used, where the hypothesis is the inverse of the ADF test, i.e., H0 postulates that the series is I (0) against the alternative that the series is I (1). The KPSS test is a Lagrange Multipliers (LM) test being represented by the model:

	(4)

Where: µt represents the deterministic trend, at a random walk and et a residual stationary variable. It is given et ∼ iid (0, σ2e). Since the residual variable (et) is stationary, thus Zt is a trend stationary process. Therefore, as the null hypothesis σ2e = 0 and the residuals et ∼ iid (0, σ2e), then et must be null. As the null hypotheses of the ADF and KPSS tests are distinct, it is expected that the two signal the same decision regarding the stationarity of the series (SOUZA; SOUZA, 2010).

2.6 The Johansen method

For the analysis of causality and cointegration, the econometric model used was Johansen cointegration, Granger causality tests and impulse-response function. The Johansen method is based on an unrestricted VAR model, represented in terms of the levels of the variables relevant to the analysis. To illustrate the process, we consider the following vector equation:

(5)

In this case, yt represents a vector of k non-stationary endogenous variables, Ai represents matrices k × k of parameters and ut represents a vector of residues iid, with zero mean and matrix of contemporaneous variances and covariances Ω. The matrix Ω is positive definite, so that the residues are not correlated in series, but may be contemporaneously correlated. System 5 is in reduced form, where each variable in yt depends on its out-of-date values, the out-of-phase values of the other system variables and the constant µ.

According to Sims (1980), this type of model has the particularity of allowing the modelling of dynamic relations between endogenous joint variables without imposing strong a priori constraints to the system, such as particular structural relations or the exogeneity of some of the variables. The VAR model also permits the historical decomposition of the variance of prediction errors, k periods ahead, in percentages to be attributed to each component variable of the system, analyzing the importance of each shock in each variable of the model, occurred in the past in the explanation of the deviations from the observed values of the variables in relation to their forecast at the beginning of the period considered.

The concept of cointegration was first introduced by Engle and Granger (1987). The economic interpretation of integration is that if two (or more) variables have a long-run equilibrium relationship, even if the series may contain stochastic trends they will move together over time and the ratio between them will be stable, that is, stationary. In short, the concept of cointegration indicates the existence of long-run equilibrium, for which the economic system converges in time (HARRIS, 1995).

Cointegration tests are vital for those who work with time series in economics since they make it possible to study and analyse structural relations between economic variables. More precisely, these tests allow us to determine if the variables have a long-term equilibrium relationship between the variables (MARGARIDO, 2004). Johansen’s methodology (JOHANSEN; JUSELIUS, 1990) considers a VAR model of order p, where it is necessary to know the number of lags used in the VAR model in an earlier step. For the determination of the number of integration vectors, the methodology determines two tests denominated trace statistic (λtraço) and maximum eigenvalue (λmax). In the trace statistic (λtraço) the null hypothesis (H0) is that there are at least r cointegration vectors. The maximum eigenvalue test (λmax) has the following premises: (H0) the number of cointegration vectors is equal to r vectors and (H1) the number of cointegration vectors is equal to r + 1 vectors. If the calculated values are higher than the critical values, we reject the null hypothesis of non-cointegration, the critical values of the trace test and the maximum eigenvalue test.

Briefly, Johansen’s cointegration methodology generally involves the following steps: (1) testing the order of integration of the model variables using the ADF and KPSS test; (2) choose the number of lags of the VAR model to be included in the cointegrating space, so that the residuals are white noise. The identification of the existence of cointegration between the studied variables and that a VAR model with error correction can be written, it is identified which are the causal relations existing between the series under study.

Causal relationships are used to identify causality relationships. The causality, according to Granger (1969), is when a variable X causes another variable Y, in the sense of Granger, if the observation of X in the present or the past helps to predict the future values of Y for some time horizon. Causality can exist in one of the directions, or X causes Y or Y causes X, or still, bi-causal. This causality means that oscillations made in X can be felt in Y or vice versa. Hall, Anderson and Granger (1992) note that cointegration between two variables is a sufficient (but not necessary) condition for the presence of causality in at least one direction.

Finally, the effect of a shock from one variable to another is verified, in order to verify the extent of the shocks from one to another without distinction, and for this, the impulse-response function is used. When applied a “shock” (an innovation of a standard deviation), in the error term of one of the variables of the system of equations and if this shock does not produce an effect on the predictions of variance of the error of the other variables, it can be said that this variable is exogenous, that is, it is independent of the others. If we consider a variable Xt, Yt to verify the existence of a cointegration and causality relationship between them, we observe that the effect of a shock on Yt not only changes immediately the values of the variable Xt, but also the future values of Xt and Yt, provided that the lagged values appear in the two equations.

The impulse-response function primarily represents the behaviour of a variable when another system variable, or itself, suffers a shock (impulse) at that particular time t, which is transferred to the future period, at t + 1, t + 2, and so on. In other words, the results presented in the impulse-response function allow to adequately evaluate shock results in any of the system variables.

The Cholesky decomposition will be used to identify the vector et. A shock in the variable does not only directly affect the ith of the model, but is also transmitted to the other endogenous variables through the dynamic structure of the VAR.

An impulse-response function traces the effect of a shock or an innovation on the current and future values of the model’s endogenous variables. The construction of the impulse-response function starts from the representation of the VAR model, which are expressed in current and past values of the shocks eti. Brooks (2003) points out that the impulse-response function makes it possible to determine if the changes of a given variable have positive or negative effects on the other variables of the system, as well as to determine the time required for that effect of being gathered. Thus, for each variable of the equation, a unitary shock is submitted to the disturbance. The effects on the system are represented graphically. If there are g variables in the system, a total of g2 shocks may be generated.

Bliska (1990) reports that one of the main advantages of orthogonalised innovations over others is that they are uncorrelated. However, there is a different decomposition for each ordering of the variables, and the direction of the effect is derived from the arbitrary selection of the order of the variables in the analysed vector. Therefore, the lower the contemporary covariance (the lower the correlation between the residues), the smaller the importance of the selected order. Therefore, even if there is no sense of causality between two variables, there may still be an effect of a shock on one of them over the other due to the presence of covariance between their respective errors.

3 Metodology

Data collected on a monthly basis refers to the price of milk paid to the producer and milk production in Rio Grande do Sul, from January 1995 to January 2017. These data were obtained from IBGE – Instituto Brasileiro de Geografia e EstatísticaI. The methodology is composed of three primary phases. The first phase deals with the series investigation of stationarity. The second phase consists of modelling the forecast, using the methodology proposed by Box and Jenkins (1976), where the milk price and production series were modelled through ARIMA models.

According to Morretin & Toloi (1987), the construction of Box-Jenkins models is based on an iterative cycle, in which the choice of the model is made based on its own data. According to Box & Jenkins (1976), three are the steps to build the model: 1) Identification: consists of discovering which of the various versions of the Box-Jenkins models, whether seasonal or not, describe the behaviour of the series. The identification of any estimated model occurs by the behaviour of autocorrelation functions (ACF) and by partial autocorrelation functions (PACF). Further details on how to obtain those functions and what behaviours represent the models can be researched in Makridakis et al. (1998). 2) Estimation: consists in estimating both parameters φ and Φ of the auto-regressive component, both parameters θ and Θ of the moving average component and the variance of εt. 3) Verification: consists of assessing whether the estimated model is adequate to describe data behaviour. If the model is not suitable, the cycle is repeated, returning to the identification phase. A widely used procedure is to identify not only a single model, but some models that will then be estimated and verified. When a satisfactory model is obtained, one moves on to the last stage of the Box-Jenkins methodology, which is also its main objective: to make forecasts. After the analysis of the residues of the competing models that explain the phenomenon under study, the model selected will be the one that meets the criteria of CTA and BIC; the lower value of CTA or BIC would signal the most likely model to make forecasts. It is worth remembering to always look for models that provide white noise, that is, with zero mean, constant variance, and which are not self-correlated. The information criteria CTA and KBS are calculated based on the estimated variance. As such criteria take into account the number of parameters used in modelling, they are called penalty criteria (MADDALA, 1992). Algebraically, they are defined as:

	(6)
	(7)

Where: T is the sample size; 1 is the value of the likelihood log function and k the number of estimated parameters.

These criteria will also be used to determine the number of lags to be incorporated in the VAR model. Besides explaining the phenomenon under study, the chosen model shall present white noise in the waste analysis (waste must present zero mean and constant), not auto-correlated variance and better statistics in the AIC penalty criteria (Souza & Lopes, 2009). The selected model will have its residues investigated in relation to homoscedasticity through the LM-ARCH test. These residues will be modelled whenever the presence of residual volatility in the variable under study is verified. The class of ARCH models will be employed, and the adequacy of the model will be investigated again, through analysis of residues. Market integration analysis is important for a country’s economic development; in view of such relevance, plus the role played by the dairy sector on national and global scale, Johansen’s is the method chosen for the second phase of data modeling.

Johansen method consists in the study of the cointegrating characteristic (r) of the VAR system. For this purpose, VAR represented in (8) can be written in the form of a vector error correction model (VECM), which was first introduced by Sargan (1964):

(8)

Where e . This specification of the system contains information on the parameters of adjustment of the model in the short and long term through estimates of ^Γk and ^Π. The use of the VEC model in the context of cointegration is ensured by the Granger Representation Theorem, which postulates that “if there is a linear dynamic model with stationary perturbations and data is I (1), then variables are cointegrated in order CI (1, 1)”. A more general version of the cointegration concept says that “if there is a matrix β such that the waste generated by the regression of the data is I(d−b), with data I(d), with d ≥ b, so that variables are to be cointegrated of CI(d,b)order”.

In terms of economic interpretation, the concept of cointegration tells us that if a chronological succession vector is bonded through a long-term equilibrium relationship, then involved variables move together over time, regardless of the existence of stochastic trends in each of them, so that their difference is stable. The vector error correction model has a great advantage – it separates both short- and long-term information involved in the adjustment process. If xt ∼ I(1), then ∆xt∼I(0) e Γk ∆xt-k∼I(0). Πxt-1 is a linear combination of variables I (1), while this one I (0) is based on the assumptions made in relation to the random disturbances of the long-term model. Matrix Π can be factored as Π=α β′, where α stands for the speed of balance adjustment and β is the matrix of long-term coefficients, i.e. the cointegrating vectors. This is valid when there are r cointegrating vectors with 0<r<I that represent the error correction mechanism in the VAR system. Once the number of cointegrating relationships in the system has been determined and matrices estimated α and β, VAR is estimated incorporating these cointegrating relationships. When the cointegrating characteristic is equal to the number of endogenous variables in the system (r=i), variables in levels become stationary and the usual methods for VAR estimation can be applied. If r=0, then Π=0. In this case, there is no cointegrating relationship between the variables in the system and it is only possible to apply VAR in the first differences with no long-term elements. When 0<r<i,, there is cointegration among system variables. Determining the cointegrating characteristic of the system is equivalent to establishing how many cointegrating vectors there are in β, that is, how many columns of α are null. It means determining the number of linearly independent lines in matrix Π. Johansen (1988, 1991) proposed two tests for the null hypothesis: the trace and the maximum eigenvalue tests. In the first case, the alternative hypothesis is that the characteristic is i,, while in the second case, it is r+1. Two important characteristics of cointegration motivate its use: the cointegration test, which estimates a coefficient capable of capturing the long-term relationship between the variables; and the short-term dynamics, which is estimated for any long-term relationship found through the inclusion of the error correction term. Soon after, to identify the causality relationships among variables, the Granger causality test is employed, followed by the impulse-response function analysis, that checks if the change of a single variable has positive or negative effects on the other.

4 Results and Discussions

Descriptive analysis is important to identify variables behaviour, for this purpose we analyse means, standard deviation, coefficient of variation, asymmetry, kurtosis, and the normality of data through Jarque-Bera’s test. According to Souza (2016), Jarque-Bera statistic is based on the asymmetry and kurtosis coefficients. If p-value is greater than the 5% significance level (whose critical value is 5.99), null hypothesis of normality will not be violated. Table 1 presents some descriptive statistics on the price and production of the paid milk to rural producers in Rio Grande do Sul from January 1995 to January 2017.

Table 1 – Descriptive measures of milk prices paid to producers in Rio Grande do Sul

	Mean	Std. Dev.	CV	Asymmetry	Kutosis	J-B	p-value
Price	0.556	0.310	0.557	0.787	3.178	27.709	0.000
Production	189774	75132	0.396	0.259	1.328	15.629	0.000

Source: Data formulated by the author

Note: Production is measured in litres

By analysing the variation coefficient for the price of milk paid to producers, we can interpret that there is a large dispersion about the mean, with deviations reaching 55% of the mean value. As for production, we observe that deviations reach 39% of the mean value. Percentages show the weight of Standard Deviation over the series.

When analysing the value of asymmetry and kurtosis, we realize asymmetry is different from zero, and both series have a positive asymmetry: one in the price (0.777) and the other in the production (0.259). Kurtosis (3,178 in price and 1,328 in production) infers that the distribution may have a heavy burden the price variable. Values of asymmetry and kurtosis are indicative there may be signs of volatility in the price and production series. In this study results do not point at it, though.

The next section brings the study of the temporal evolution to the stationarity and non-stationarity analysis of chronological successions (price paid to the producer and milk production in Rio Grande do Sul, from January 1995 to January 2017). ADF and KPSS tests are then performed.

In Figures 2 and 3 one can analyse the monthly averages of milk price and production in Rio Grande do Sul, where both series are growing over time. It can be observed that price and production follow an increasing trend over time. Graphically, the series have non-stationarity characteristics.

Figure 2 – Monthly means for milk price

Source: Authors (2020)

Figure 3 – Monthly means for milk production

Source: Authors (2020)

After the graphical analysis it is necessary to investigate the stationarity of the series as a preliminary analysis for the modelling of time series, based on Lütkepohl and Kratzig (2004), because it is necessary to consider such characteristics in the modelling of the process generating the data that compose a system of variables potentially related in time.

The stationarity check was performed using the ADF hypothesis test and as a form of validation of the ADF test the KPSS test is used, where the hypotheses presented are inverse of the ADF tests, i.e., H0 postulates that the series is I (0 ) against the alternative that the series is I (1), according to Souza and Souza (2010).

Table 2 – ADF and KPSS unit root tests at level and first differences

	ADF	∆(ADF)	KPSS	∆ (KPSS)
Price	-0.100	-10.829*	0.034*	0.135*
Production	-0.392	-7.055*	2.071	0.143*

Source: Data processed by the authors.

Notes: ADF: * Significant at the 5% level. Critical values of MacKinnon (1996): -3,501 (1%) and -2,892 (5%). KPSS values: Critical values of MacKinnon (1992, table 1): 0.739 (1%) and 0.463 (5%).

When the level series are analysed by the ADF hypothesis test, according to Table 2, it can be affirmed that the series are non-stationary and these results are corroborated by the KPSS test, becoming stationary in first difference.

In Figures 4 and 5 it is observed that the series in level do not have a steady behaviour, that is, they do not distribute around a mean, and after the differentiation the series become stationary.

Figure 4 – Price of milk in level and differentiated

Source: Authors (2020)

Figure 5 – Milk production in level and differentiated

Source: Authors (2020)

The next step is to analyse the autocorrelation (FAC) and partial autocorrelation (FACP) functions. The FAC and FACP signal information as the type of filter that will be used by the behaviour of these functions (AR, MA, ARMA or ARIMA) and also suggests the order of the model by the Box and Jenkins methodology.

Figure 6 – Autocorrelation (FAC) and partial autocorrelation (FACP) functions of the original milk price series, in sequence

Source: Authors (2020)

Figure 7 – Autocorrelation (FAC) and partial autocorrelation (FACP) functions of the original milk production series, in sequence

Source: Authors (2020)

Figures 6 and 7 presents respectively the FAC and FACP of the price of milk paid to the producer and of milk production in Rio Grande do Sul. It is noted that the FAC of the price and the FAC of the production show a slow decay to zero, therefore have non-stationarity characteristics.

The next methodological procedure, concerning the forecast models, is the identification of the significant models that represent the behaviour of the price and milk production series in Rio Grande do Sul. Through the analysis of Tables 3 and 4 it is possible to observe the competing models for price and milk production, respectively, which presented white noise and all significant parameters.

Table 3 – Models competing for the price of milk in Rio Grande do Sul

Models	Parâmetros	AIC	SBC
ARIMA (1, 1, 0)	Φ = 0.497	-1165.125	-1115.062
ARIMA (1, 1, 1)	Φ = 0.253 Θ = -0.377	-1181.472	-1127.832
ARIMA (0, 1, 1)	Θ = 0.535	-1176.529	-1126.466

Source: Authors (2020)

Among these models, it is necessary to choose the best one, in this way, we used the AIC and SBC penalty criteria, where the model that presented the lowest value for AIC and SBC was the ARIMA model (1, 1, 1) in the price series.

Table 4 – Competing models for milk production in Rio Grande do Sul

Models	Parameters	AIC	SBC
ARIMA (1, 1, 1)	Φ = 0.139 Θ = -0.029	5902.528	5916.832
SARIMA (1, 1, 1) (1, 1, 1)6	φ = 0.553 θ = -0.311 Φ = -0.752 Θ = -0.577	5285.477	5303.125
SARIMA (1, 1, 0) (1, 1, 1)6	φ = 0.574 Φ = -0.607 Θ = -0.781	5324.271	5341.918

Source: Authors (2020)

Residual series were tested for conditional heteroscedasticity, by employing ARCH LM test; statistic F resulted in 5.5438, not pointing any evidence for conditioned heteroscedasticity. Thus, ARIMA model is a feasible option to model the series under study and make predictions.

In the same way, the milk production series model was evaluated, where the model that presented the lowest AIC and SBC was the SARIMA model (1, 1, 1) (1, 1, 1)6. In Figures 8a and 8b, the FAC and FACP of the residues of the ARIMA (1, 1, 1) model chosen for the price variable are shown.

Figure 8 – (a) Autocorrelation functions (FAC) and (b) partial autocorrelation (FACP) of the residues of the milk price series model, ARIMA (1, 1, 1)

Source: Authors (2020)

Figures 9a and 9b show the FAC and FACP of the residues of the SARIMA model (1, 1, 1) (1, 1, 1)6 for the variable milk production.

Figure 9 – (a) Autocorrelation functions (FAC) and (b) partial autocorrelation (FACP) of the residues of the milk production series model, SARIMA (1, 1, 1) (1, 1, 1) 6

Source: Authors (2020)

The residual FACs of the price series models and the production series show that autocorrelations of residues are within 95% of the confidence limits, as occurs in FACPs. Therefore, these wastes are not autocorrelated. Thus, to confirm whether the residues are autocorrelated or not, the Ljung-Box Q * test is used, the residuals are considered as non-autocorrelated when the p-value is higher than the level of significance adopted.

After analysing the price series (p-value: 0.275) and the production series (p-value: 0.073), it was verified that in both series the p-values were higher than the level of significance adopted considering non-autocorrelated residues. For the models found, one can see in Table 5 the monthly forecasts, from February to July 2017 and, in Table 6, the monthly forecasts, for the same period, for the milk production series.

Table 5 – Milk price forecasts in Rio Grande do Sul in the year 2017

Months	Forecast	Real value	Limit lower than 95%	Limit higher than 95%	Standard error	VR-VP
Feb./17	1.28	1.32	1.23	1.33	0.03	0.04
Mar./17	1.29	1.34	1.20	1.39	0.05	0.05
Apr./17	1.32	1.35	1.19	1.45	0.07	0.03
May/17	1.34	1.38	1.19	1.50	0.08	0.04
Jun./17	1.36	1.36	1.18	1.54	0.09	0.00
Jul./17	1.38	1.32	1.18	1.58	0.10	-0.06

Source: Authors (2020)

The period of analysis was from January 1995 to January 2017 and a six months prognosis was made, corresponding to the period from February to July 2017. Analysing Table 5 it is observed that the model captured with a minimum error, obtaining very accurate forecasts, proving to be an excellent tool, allowing the rural producer financial planning that will also allow the development of his business. The applied methodology was able to predict, quite accurately, the values analysed in the worked period.

Table 6 – Forecast of milk production in Rio Grande do Sul in 2017

Months	Forecast	Real value	Limit lower than 95%	Limit higher than 95%	Standard error	VR-VP
Feb./17	246305	247111	229982	262628	8328	806
Mar./17	258253	256719	237349	279158	10665	-1534
Apr./17	233132	236390	209392	256872	12112	3258
May/17	239773	254331	213915	265631	13193	14558
Jun./17	256289	267698	228673	283905	14089	11409
Jul./17	278578	304528	249411	307746	14881	25950

Source: Authors (2020)

Analyzing Table 6, which presents the prediction of milk production in Rio Grande do Sul for six months, corresponding to the period from February to July 2017, it is noticed that in April there is a decrease in production, but already in May it is recovering again, and the milk intake in that region continues a growing tendency for the following months.

In reporting the empirical results of the forecast models, the second phase of the proposed methodology is the hypothesis of cointegration and causality between the prices paid to the consumer and the milk production. The results of the cointegration tests between the price and the production are reported in Table 7. The order of the number of lags inserted in the models was made based on the lag order selection criteria (the test calculates several criteria for selecting the delay order of an unrestricted VAR model. The table indicates the selected discrepancy of each column criterion by an asterisk (*). All criteria are discussed in Lütkepohl (1991, Section 4.3), and it was considered the smallest value of the SBC criterion for adjusting the model with a lag for each variable.

Table 7 – Results of the test lag order selection criteria between milk production and price paid to the producer

Lag	AIC	SC	HQ
1	17.96975	18.05307	18.00327
2	17.84477	17.98364	17.90063
3	17.74935	17.94377*	17.82755
4	17.70260	17.95257	17.80315
5	17.68656	17.99208	17.80945
6	17.65421	18.01528	17.79944
7	17.66100	18.07762	17.82858
8	17.50808*	17.98024	17.69800*

Source: Authors (2020)

After the number of lags is defined, we can evaluate the presence of cointegration vectors by the Johansen method, in which the cointegration relations of the variables were defined for the three lags model. The cointegration term is known as the error correction term because the long-term equilibrium deviation is gradually corrected through a series of short-term partial adjustments.

The cointegration test indicates cointegration between the price of milk paid to the producer and the production itself; that is, there is a long-term equilibrium relation between the variables. In this regard, it can be said that the price paid to the producer match the milk production.

When the behaviour of the series is analysed over time, it is verified that when production increases, then the price decreases, as a consequence of the law of supply and demand: when the supply is higher than the demand, the prices fall as a result.

Table 8 – Cointegration value

Hip.	Eigenvalue	λTrace	0,05 Critical value	Prob.**	λmax	0,05 Critical value	Prob.**
0*	0.364751	210.0985	20.26184	0.00	117.9719	15.89210	0.00
Until 1*	0.298359	92.12657	9.164546	0.00	92.12657	9.164546	0.00

Source: Authors (2020)

Note: * Denotes rejection of the hypothesis at the 0.05 level. ** Mackinnon, Haug and Michelis (1999).

From the modelling, the Granger causality test was performed among the variables under study. Analysing Table 9 it can be seen that the results presented are all significant at 5% significance. The results point to a bidirectional causality, that is, the prices paid to the producer influence milk production and vice versa. Shocks or fluctuations in prices are transferred between price and production.

Table 9 – Granger causality test for milk production and price paid to the producer

Null hypothesis	Direction	Lags*	Teste F	p-value	Decision	Causality
Milk Production doesn’t cause PPP	-->	3	7.64266	0,00	Rejects H0	Bidirectional
PPP doesn’t cause milk production	<--		4.51367		Rejects H0

Source: Authors (2020)

Note: Price Paid to the Producer (PPP).

After the causality analysis, the Cholesky impulse-response function is performed, which shows the long-term effects in the time series. The analysis of the impulse-response function is particularly useful for the observation of the direction, duration and reaction pattern of the variable response of interest to the impulses of a standard deviation and future in the endogenous variables of the model.

Graphical analysis of the impulse-response function is a practical way of visualising the behaviour of the dependent variable in response to a series of shocks. It is important to emphasise that the ordering of the variables has a significant influence on the analysis of the impulse-response function. The importance of the classification depends on the magnitude of the correlation coefficient between the error terms.

Figure 10 shows the impulse-response functions of the price paid to the producer and the milk production to the innovations of a standard deviation. It is found that a shock of standard deviation applied to the price paid to the producer influences the milk production. Thus, in the first month, the price decreases until the fourth month, but from the second period, it rises again until the third month and stabilises in the short term. On the other hand, milk production, when price changes occur, decreases up to the fourth month, increasing between the third and fourth month, returning to stability between the seventh and eighth month.

Figure 10 – Impulse response function analysis

Source: Authors (2020)

Note: Price Paid to the Producer (PPP).

The first graph shows a response to price variation and a shock in itself. It is observed that the shock produces a relatively high-intensity response, with a relatively long negative duration, ceasing the impact from the eight-month after an innovation. The second graph shows a response in the price paid to the producer to a shock in the production variation. It is verified that the shock produces a relatively low-intensity response and also of negative medium-duration, ceasing the effect starting from the sixth month after the innovation.

The third graph shows the milk production response to a shock in the price variation paid to the producer. It is found that the shock produces a low-intensity response and also a short one. The fourth graph shows the response of the milk production variation to a shock in itself. It is observed that the shock produces a response of relatively medium intensity, with a relatively long negative duration, ceasing the impact from the seventh month after the innovation.

The analysis of the impulse-response function showed that the shocks present significant impacts between milk production and the price paid to the producer, both regarding duration and intensity. In general, shocks showed six to eight-month duration.

5 Conclusions

The research was carried out through the analysis of time series, using data from milk prices and production practised in Rio Grande do Sul from 1995 to 2017. In developing the study on price forecast and milk production, producers were provided with a component in order to aid decision-making in production management.

Thus, management measures can be taken by knowing the behaviour of this variable in advance. The predictions found, as expected, show a minimal difference between actual and predicted values, proving the efficiency of the Box and Jenkins (1976) methodology.

Subsequently, the Johansen cointegration method was used. The test indicates cointegration between the price of milk paid to the producer and the production; that is, there is a long-term equilibrium relation between the variables. Therefore, it can be said that the price paid to the producer is synchronised with milk production.

Granger’s causality test points to a two-way causality, i.e. producer prices influence milk production and vice versa. Shocks or fluctuations in prices are transmitted between price and production. The analysis of the impulse-response function showed that the shocks present significant impacts between milk production and the price paid to the producer, both in terms of duration and intensity. In general, shocks showed durations of six to eight months.

The proposed methodology provides information and subsidies based on robust analyses. Therefore, it offers rural producers and/or decision makers pre-management measures for subsequent decision-making and planning, although it is important to consider that, in practice, producers take the prices as they are instead of setting the prices themselves, as the industry is the one that decides the amount to pay per litre of milk. Finally, it is recommended as significant points for future work the use of other Brazilian states. In order to carry out international comparisons, it is possible to use countries that are part of the Mercosur in order to evaluate the cointegration of the market in the sector.

Acknowledgment

We want to thank Universidade Regional Integrada do Alto Uruguai e das Missões – URI for funding the research translation.

Thanks to reviewers and collaborators.

References

ANUALPEC. Anuário da Pecuária Brasileira. Pecuária de leite. São Paulo: IEG/FNP; 2017.

BLISKA FMM. Formação de preços de carne bovina: uma aplicação do modelo de auto-regressão vetorial.  Agricultura em São Paulo. 1990;37(3):41-59.

BOX, G. E. P.; JENKINS, G. M. Time Series Analysis: Forecasting and Control. San Francisco: Holden-Day, 1970.

BOX, G. E. P.; JENKINS, G. M.; REINSEL, G. C. Time series analysis: forecasting and control. 3. ed. San Francisco: Holden-Day, 1994.

BOX GEP, JENKINS GM. Times series analysis: forecasting and control. Holden-day, San Francisco. 1976

BROOKS C.  Introductory econometrics for finance. Cambrige: University Press; 2003.

CAMELO HN et al. Modelagem box – jenkins aplicada a previsão de velocidade do vento em regiões do nordeste brasileiro para fins de geração eólica. Ciência e Natura. 2018;40(30):1-10.

CEPEA – Centro de Estudo Avançados em Economia Aplicada. Boletim Informativo. São Paulo; 2017. Available from: https://www.cepea.esalq.usp.br/br

DALZOTTO-ARTUZO et al. Relação entre os custos de produção e o preço de mercado do milho. Custos e @gronegócio on line. 2017; 13(2):448.

EMBRAPA - Empresa Brasileira de Pesquisa Agropecuária. Gado de leite: importância econômica. São Paulo; 2017. Available from: https://sistemasdeproducao.cnptia.embrapa.br/FontesHTML/Leite/LeiteCerrado/importancia.html

EMBRAPA - Empresa Brasileira de Pesquisa Agropecuária. Anuário Leite 2018: indicadores, tendências e oportunidades para quem vive no setor leiteiro. São Paulo; 2018. Available from:

https://www.embrapa.br/busca-de-publicacoes/-/publicacao/1094149/anuario-leite-2018-indicadores-tendencias-e-oportunidades-para-quem-vive-no-setor-leiteiro

ENGLE RF. GRANGER CWJ. Co-integration and error correction: representation, estimation, and testing. Econometrica. 1987:55(2):251-276.

FAO - Food and Agriculture Organization of the United Nations. Dairy production and products: milk production. 2016. Available from: http://www.fao.org/agriculture/dairygateway

FAVA VL. Metodologia de Box-Jenkins para modelos univariados. In: Vasconcellos MAS, Alves D, editors. Manual de econometria: nível intermediário. São Paulo: Atlas; 2000.

FERRAZ, M SÁFADI T, LAGE G. Uso de modelos de séries temporais na previsão de séries de precipitação pluviais no município de Lavras – MG. Revista Brasileira de Agrometeorologia. 1999;7(2):259-267.

GRANGER CWJ. Investigating causal relations bu econometric models and cross – spectral models. Econometrica. 1969;37(3):424-438.

HALL AD, ANDERSON HM, GRANGER CWJ. A cointegration analysis of Treasyry Bill Yields. The Review of Economics and Statistics. 1992;74(1):116-126.

HARRIS RID. Cointegration analysis in econometric modelling. London: Prentice Hall, 1995.

IBGE - Instituto Brasileiro de Geografia e Estatística. Rio de Janeiro; 2017. Contas nacionais trimestrais. Available from: https://www.ibge.gov.br/

IBGE - Instituto Brasileiro de Geografia e Estatística. Rio de Janeiro; 2017. Pesquisa da pecuária municipal 2014. Available from: https://www.ibge.gov.br/

JOHANSEN S, JUSELIUS K. Maximum likelihood estimation and inference on cointegration: with applications to the demand for money. Oxford Bulletin of Economics and Statistics; 1990.

KOVAČEVIĆ V et al. Causality between corn production cost and cash corn price. Custos e @gronegócio on line. 2017;13(4):1-15.

LEAL JÚNIOR BV et al. Modelo híbrido de previsão de séries temporais para possíveis aplicações no setor de geração eólica. Ciência e Natura. 2018;40(Edição Especial: X Workshop Brasileiro de Micrometeorologia):1-5.

LÜTKEPOHL H. Introduction to multiple time series analysis. Berlin: Springer; 1991.

LÜTKEPOHL H, KRATZIG M (eds.). (2004) Applied Time Series Econometrics. Cambridge: Cambridge University Press; 2004.

MACHADO SLD, SILVA CR, ARAÚJO AA. Descrição temporal do comportamento do Cerrado sensu strictu usando séries temporais. Ciência e Natura. 2018;40(30):1-10.

MACKINNON JG. Model specification tests and artificial regressions. Journal of Economic Literature. 1992;30(1):102-146.

MACKINNON JG. Numerical Distribution Functions for Unit Root and Cointegration Tests. Journal of Applied Econometrics. 1996;11(6):601-618.

MACKINNON JG, HAUG AA, MICHELIS L. Numerical distribution functions of likelihood ratio tests for cointegration. Journal of Applied Econometrics. 1999;14(5):563–577.

MADDALA GS. Introduction to econometrics. New Jersey: Prentice-Hall; 1992.

MARCHEZAN A, SOUZA A. Previsão do preço dos principais grãos produzidos no Rio Grande do Sul. Ciência Rural. 2010; 40(11):2368-2374.

MARGARIDO MA Teste de co-integração de Johansen utilizando o SAS. Revista Agrícola. 2004;51(1):87-101.

MAPA - Ministry of Agriculture, Livestock and Supply. Brasília; 2017. Available from: http://www.agricultura.gov.br/

MORETTIN, P. A. e TOLOI, C. M. C., Modelos de Função de Transferência, Rio de Janeiro, 3ª Escola e séries Temporais Econometria, 1989.

NOMELINI QSS et al. Uso de modelagem univariada e multivariada com séries temporais como ferramenta de gestão do agronegócio na cultura de soja do Brasil. Revista Espacios. 2017;38(8).

OECD - Organisation for Economic Co-operation and Development. Projections agricoles 2018-2027. Paris; 2018. Available from: https://www.oecd-ilibrary.org/agriculture-and-food/data/statistiques-agricoles-de-l-ocde/perspectives-agricoles-de-l-ocde-et-de-la-fao-edition-2018_371b829bfr?parentId=http%3A%2F%2Finstance.metastore.ingenta.com%2Fcontent%2Fcollection%2Fagr-data-fr

OECD/FAO. Perspectives agricoles de l’OCDE e de la FAO. Estatísticas agrícolas da OCDE (base de dados), Paris-França. Paris; 2018. Available from: http://dx.doi.org/10.1787/agr-data-fr

PINHEIRO CAO, SENNA V. Previsão de preços através de redes neurais e análise espectral: evidências para o mercado futuro das commodities açúcar e soja. Custos e @gronegócio on line. 2017;13(4):1-26.

RETES-MANTILLA RF, TORRES-MANCERA MT, HERNÁNDEZ MLC Estimación del precio de compra de la pulpa de café en México para su aprovechamiento en la obtención de productos de alto valor agregado. Custos e @gronegócio on line. 2017; 13(3):101-119.

SANTANDER. A importância do agronegócio para o Brasil. 2017. Available from: https://www.santandernegocioseempresas.com.br

SIMS C. Macroeconomics and Reality. Econometrica. 1980;48(1):1-48.

SOUZA FM et al. (2010). Previsão do consumo de cimento no Estado do Rio Grande do Sul. Pesquisa Operacional para o Desenvolvimento. 2010; 2(1):1-86.

SOUZA FM. Modelos Box & Jenkins aplicados à previsão de demanda de leitos hospitalares [monography]. Santa Maria: Programa de Pós-Graduação em Estatística e Modelagem Quantitativos/UFSM; 2006.

SOUZA FM. Modelos de previsão: aplicações à energia elétrica – ARIMA – ARCH – AI e ACP. Curitiba: Appris; 2016.

SOUZA FM, SOUZA AM. Previsão para o fornecimento de energia elétrica para o RS por meio de combinação de previsão utilizando ponderação por autovalores. X Semana de Engenharia de Produção Sul Americana - SEPROSUL; 2010; Santiago; Chile; USACH – Universidade de Santiago do Chile, 1:1-10; 2010.

SOUZA, F. M.; LOPES, L. F. D. Previsão de demanda de leitos hospitalares por meio de análise de Séries Temporais. Revista Ciência e Natura. v. 31, n. 1. jun. 2009. p. 33-47.

TERRA VIVA. Perspectivas agrícolas da OCDE e da FAO 2017-2026 – Produção. 2017. Available from: http://www.terraviva.com.br/site/index.php?option=com_k2&view=item&id=12467

TIBULO C, CARLI V. Previsão do preço do milho através de séries temporais. Scientia Plena. 2014;10(10).

TIBULO C, TIBULO V. Modelos de séries temporais para a previsão do preço médio mensal da soja no Rio Grande do Sul e análise da evolução da cultura no cenário nacional e regional. 2013. Available from: http://www.ucs.br

WANKE P, BARBOSA M. Combinação ou competição de previsões: um estudo de caso nos fretes do agronegócio. 2010. Available from: https://www.ilos.com.br/web/combinacao-ou-competicao-de-previsoes-um-estudo-de-caso-nos-fretes-do-agronegocio-parte-1/

WERNER L, RIBEIRO JLD. Previsão de demanda: uma aplicação dos modelos Box-Jenkins na área de assistência técnica de computadores pessoais. Gestão e Produção. 2003;10(1):47-67.

Authorship Contribution

1 - Francisca Mendonça Souza

Contribuições: Curadoria de Dados, Análise Formal

2 - Brunna Aver de Paula

Contribuições: Curadoria de Dados, Investigação, Análise Formal, Validação

3 - Daniel Knebel Baggio

Contribuições: Curadoria de Dados, Análise Formal

4 - Argemiro Luis Brum

Contribuições: Investigação, Análise Formal

5 - Juliana da Fonseca Capssa Lima Sausen

Contribuições: Metodologia, Análise Formal

---