Modelling Risk for Commodities in Brazil: An Application for Live Cattle Spot and Futures Prices

Abstract

This study analyses a series of live cattle spot and futures prices from the Boi Gordo Index (BGI) in Brazil. The objective is to develop a model that best portrays this commodity’s behaviour to estimate futures prices more accurately. The database created contains 2010 daily entries in which trade in futures contracts occurs, as well as BGI spot sales in the market, from 1 December 2006 to 30 April 2015. One of the most important reasons why this type of risk needs to be measured is to set loss limits. To identify patterns in price behaviour in order to improve future transaction results, investors must analyse fluctuations in asset values for longer periods. Bibliographic research reveals that no other study has conducted a comprehensive analysis of this commodity using this approach. Cattle ranching is big business in Brazil given that in 2021, this sector moved BRL 913.14 billion (USD 169.29 billion). In that year, agribusiness contributed 26.6% of Brazil’s total gross domestic product. Using the proposed risk modelling technique, economic agents can make the best decision about which options within these investors’ reach produce more effective risk management. The methodology is based on Holt–Winters exponential smoothing algorithm, autoregressive integrated moving-average (ARIMA), ARIMA with exogenous inputs, generalised autoregressive conditionally heteroskedastic and generalised autoregressive moving-average (GARMA) models. More specifically, five different methods are applied that allow a comparison of 12 different models as ways to portray and predict the BGI commodity behaviours. The results show that GARMA with order c(2,1) and without intercept is the best model. Investors equipped with such precise modelling insights stand at an advantageous position in the market, promoting informed investment decisions and optimising returns.

1. Introduction

Cattle ranching is a major industry in Brazil. According to the report from [1], in 2021, this sector accounted for BRL 913.14 billion (approximately USD 169.29 billion). Ref. [2] highlighted that, in the same year, agribusiness contributed 26.6% to Brazil’s total gross domestic product (GDP). Drawing from past data, in 2017, livestock profits represented 31% of agribusiness contributions to the GDP. Beef exports have been crucial in maintaining Brazil’s positive trade balance, with this commodity accounting for 3.2% of all exports in 2017. There was a notable growth that year, with an increase of 9.6% in volume and 13.9% in sales compared to the previous year (see [3]).

The term ‘risk’ in this context is associated with the likelihood of unexpected events occurring in business ventures, financial investments or any other situations in which large financial losses may occur. For example, auto insurance companies need to estimate the occurrence of claims in order to price insurance premiums offered in the market. Equity investment must also take into account the potential for sharp declines in local and international markets.

Similarly, stock market risk management’s goal is to calculate capital requirements as accurately as possible by minimising potential losses [4,5]. An unrealised loss can result from the decision to keep financial assets in portfolios even when the assets’ values have dropped significantly. In this case, investors prefer to maintain their position with a devalued asset, as they expect its value will recover, reaching a break-even point (i.e., no loss or purchase price). A question that arises in this situation is how to quantify the value of potential overly high losses [5].

A decrease in share value reduces prospective gains or, conversely, creates an increase in unrealised losses. When confronted with this situation, investors must choose between two contrasting decisions: liquidate the share to conserve the capital or keep it in hopes of recovering the amount invested [6]. Unrealised losses may thus be recovered by retaining portfolio assets rather than selling and accepting the consequent realised losses. However, if assets continue to depreciate, the unrealised losses can become quite significant. Successful decisions reflect investors’ ability to perceive fluctuations in their position’s value and minimise the likelihood of making wrong decisions whether by buying, selling or maintaining their financial position [7].

One of the main reasons why this type of risk needs to be measured is to facilitate setting loss limits. To identify patterns in price behaviour and improve future transactions’ results, investors must analyse fluctuations in assets’ values for longer periods. Investors can thus define the maximum potential loss that they are willing to accept during financial transactions [6].

According to [8], using value at risk (VaR) to manage an asset portfolio allows investors to position their portfolio properly based on risk aversion and typical price fluctuations. VaR also allows investors to incorporate their expectations in terms of future volatilities, enabling more accurate, improved adjustments to their stop-loss limits. VaR was conceptualised by [9] and [10] to optimise economic agents’ choices when seeking to minimise the probability of high losses. VaR is obtained by estimating loss distributions and tail quantile values that measure distributions’ tail weight.

As an extension of Gaussian autoregressive moving-average (ARMA) models, the class of generalised autoregressive moving average (GARMA) was proposed by [11]. It allows time series modelling for non-Gaussian time series data. In this sense, we can highlight some recent works focusing on predicting commodity prices based on GARMA models (see, e.g., [12,13]).

The present study sought to analyse the behaviour of a specific commodity with a significant impact on Brazil’s GDP, namely, the Boi Gordo Index (BGI), in order to help investors who trade futures contracts. To this end, the concept of VaR was applied to examine a price series of BGI futures contracts using different methodologies and compare the respective forecasts. The research’s goal was to help economic agents make the best decision about the options within these investors’ reach to ensure more effective risk management.

The analysed database contained 2010 daily entries in which BGI futures and spot sales were traded on the market, from 1 December 2006 to 30 April 2015. Among the 2010 daily entries under study, 1994 were used to develop models, and the remaining 16 were kept apart for later comparison between actual data and model predictions produced by the different methodologies. Futures prices were obtained with GrapherOC software (https://www.quantsis.com.br/GrapherOC/compare.asp), and BGI spot prices were accessed on the Centro de Estudos Avançados em Economia Aplicada—Universidade São Paulo (Centre for Advanced Applied Economic Studies-University of São Paulo) (CEPEA-USP) website (see http://cepea.esalq.usp.br/boi/, first accessed on 2 May 2015). The modelling was performed using R software, after which the method that best addresses the research objectives and has the highest predictive capacity in terms of BGI futures prices’ VaR was selected.

The manuscript is organised as follows. Section two introduces the BGI with reference to commodities and futures prices, as well as relevant definitions. Section three first describes the data analysis process and then the exponential smoothing algorithms, autoregressive integrated moving-average (ARIMA), ARIMA with exogenous inputs (ARIMAX), generalised autoregressive conditionally heteroskedastic (GARCH) and generalised autoregressive moving-average (GARMA) models. The final section presents the conclusions.

2. BGI and Futures Prices

A commodity is merchandise (i.e., raw material or primary merchandise) that is subject to a few or even no change processes, which is produced in large quantities and priced according to international supply and demand. This type of goods’ quality can vary slightly, but commodities are largely uniform and independent of their provenance so they can be traded on international stock exchanges [14]. Various types of commodities can be found in international markets, such as natural gas, gold, soy or cattle.

Brazil is a major producer of various commodities including the one under analysis in the present study—the BGI—which is a male castrated bovine carcass that has a net weight of 16 arrobas (one arroba is approximately 15 kg) or more, with a maximum age of 42 months. Futures contracts for this meat may be traded on the Bolsa de Mercadorias & Futuros (BM&F, Stock Exchange of Commodities and Futures) market and other international markets such as the Chicago Mercantile Exchange or the Australian Securities Exchange. Notably, the Australian Securities Exchange and its dynamics have been extensively researched in recent studies, such as the work by [15] that examines the predictive performance of various financial, economic and sentiment variables for the Australian All Ordinaries index equity risk premium spanning the last 28 years. Commodities’ purchase and sale are made through ‘futures’ contracts that stipulate the quantities and quality of commodities to be traded. This research focused on estimating these futures contracts’ VaR and, more specifically, BGI futures, as described in the next section.

Commodity markets have a major impact on economies worldwide, as these markets can produce fluctuations in the price of goods essential to countries’ optimal functioning (e.g., food and energy). The trade in this type of merchandise may have an effect on many other goods’ prices that would otherwise simply be subject to the laws of supply and demand. Various authors have developed models of commodities and futures price behaviour in financial markets.

One frequently used model is the GARMA process, also referred to as the Gegenbauer autoregressive moving-average model, which was proposed by [16]. This model harnesses the properties of Gegenbauer polynomials’ generating function. Ref. Ref. [17] further introduced a $k$-factor extension while conducting research, in which a two-factor GARMA model was applied to Mauna Loa atmospheric CO₂ data. The results confirm that the model has a reasonable goodness of fit and produces excellent forecasts.

[11] paved the way for innovative time series modelling with the introduction of the generalised autoregressive moving-average (GARMA) class, designed specifically for non-Gaussian data. Their work set a foundation that would be built upon in subsequent years.

[18] also used a GARMA approach to model exchange rates’ conditional mean. The cited authors estimated a GARMA model with integrated GARCH using a wavelet-based maximum likelihood estimator. Subsequently, [19] proposed a class of observation-driven time series models known as generalised autoregressive score models.

Furthermore, exploring the nature of interaction between stock market returns in multiple countries, [20] employed a multivariate generalised autoregressive conditional heteroskedasticity (MGARCH) model. Their study encompassed stock market returns of four countries: Australia, Singapore, the UK and the US. Using data from January 1992 to December 2008, their findings underscored significant positive mean spillovers, especially from the US stock market returns, impacting primarily Australia and Singapore. This research also highlighted significant volatility and cross-volatility spillovers across all four markets, suggesting limited risk reduction benefits for investors solely diversifying their portfolios with stocks from these countries.

In addition, [21] examined the links between the crude oil market and stock markets of G-7 countries. The cited researchers used multivariate GARCH models and wavelet analysis to detect a significant volatility spillover between oil and stock markets. Ref. [22] estimated GARMA and seasonal autoregressive fractionally integrated moving-average (ARFIMA) (i.e., a special case of the multi-factor GARMA) models to analyse time series for precious metals volumes. These series exhibited periodic-stochastic behaviour and the presence of long-range dependence, so the cited scholars used a multi-factor GARMA model.

[23] developed GARMA and special ARFIMA subfamily models to deal with US commodity futures time series of counts. The cited study confirmed consistently long memory structures in these highly liquid and important financial instruments and, more specifically, revealed systematic forecastable patterns in trading behaviour and liquidity. Ref. [23] also report that GARMA models present excellent forecasts. Ref. [24] further adopted a generalised long-memory GARMA model to estimate the conditional mean of electricity spot price time series. The cited researchers applied VaR, conditional VaR and other models to assess electricity market exposure and concluded that electricity market commodities can be adopted to support diversification and hedging against stock market risks.

[25], in turn, used moving average of moving-average, exponential weighted moving-average and GARCH models to conduct parametric analyses of Indian equity mutual funds’ downside risk. The cited author applied a VaR-based approach and found that time series presented considerable downside risk during the selected period of 1999–2014. Ref. [26] also examined financial markets, specifically investigating the volatility transmission from developed markets to Asian emerging markets. The latter-cited researcher used a time series from 1996 to 2015 and worked with a heterogeneous autoregressive distributed lag model to study spillover effects.

In line with the increased adoption of the GARMA model for predicting commodity prices, recent endeavours such as those by [12] and a significant study by [13] have further showcased its utility in this domain.

Similar to previous studies, the present research also applied a statistical approach using a GARMA model to analyse data from 1997 to 2015 and relied on the GARMA model’s predictive power. The VaR was extracted and then compared to the real data to see whether the proposed approach provides a reliable VaR and accurate forecasts.

3. Empirical Research on BGI

The data consisted of BGI futures prices acquired with GrapherOC software and BGI spot prices taken from the CEPEA-USP website. The futures price data referred to the period from 11 February 2005 to 30 December 2014, whereas the spot prices came from trades between 23 July 1997 and 30 April 2015.

When a future contract is due in March and another in April, for instance, both contracts will be available to buy and will have a price months before March. In other words, futures contracts’ terms overlap, and a common market strategy is also to operate using futures with short maturities. If a contract’s last month is chosen, the value in the last days will converge on the spot value, and, after launching the series, the next contract will diverge again, because it will be due in 30 days at the end of the futures contract’s term. Thus, the data included the next-to-last month of each futures contract along with the spot prices for the same dates. The resulting period used in the data modelling was from 1 December 2006 to 30 April 2015.

Figure 1 presents a graph of the daily data under study. The black, blue and red lines represent the futures price, spot price without Funrural and spot price with Funrural, respectively. ‘Funrural’ is a social security contribution tax in Brazil that is levied on gross revenue from sales of rural products. The graph shows a clear positive trend over time.

Complementing this visual portrayal,

Table 1. Descriptive statistics for futures and spot prices (with and without FUNRURAL).

Variable	Min	First Qu	Median	Mean	Third Qu	Max	SD
Future price	50.20	79.50	95.25	94.02	102.94	151.55	21.79
Spot price without FUNRURAL	51.31	77.28	93.95	93.09	102.11	150.65	22.00
Spot price with FUNRURAL	52.52	79.11	96.17	95.28	104.52	154.20	22.52

provides a statistical overview, displaying the minimum, first quartile, median, mean, third quartile, maximum and standard deviation. As anticipated, the spot price inclusive of the Funrural tax consistently registers higher values than its counterpart without the tax. This difference aligns with the inherent nature of the Funrural as an added tax. Meanwhile, the futures price demonstrates its own unique distribution and fluctuations, reflective of market dynamics and anticipations.

For the futures price, the mean stands at 94.02, with a standard deviation of 21.7942, ranging from a minimum of 50.20 to a maximum of 151.55. This indicates that there has been an overall upward trend in the ’Future Price’ from its start to its end, though the actual movement has intermediate ’ups’ and ’downs’. The seasonal components are approximately −0.04887 and 0.05295, respectively. This shows that the seasonal fluctuations (or oscillations) around the trend have a fairly small magnitude. The seasonal patterns go both above and below the trend. Finally, the random component range goes from about −3.279720 to 3.254873. This shows that the fluctuations or irregularities are relatively larger than the seasonal component. The spot price without Funrural exhibits a mean value of 93.09 and a standard deviation of 21.99908 and spans from 51.31 to 150.65. When considering the Funrural contribution, the mean spot price slightly rises to 95.28, with a standard deviation of 22.51686 and a range between 52.52 and 154.20.

Figure 2 is a scatter plot of futures prices and spot prices. This graph reveals a strong positive correlation between the two prices.

Figure 3 displays the histogram of BGI futures prices. Most prices are concentrated in the distribution’s centre, which at first suggested that they would fit translated t-Student and normal distributions, although this would entail some positive skewness. However, low prices were a point of concern. This figure shows that the lowest price has a heavier tail, so, when the GARMA models were estimated (see Section 3.4), a t-Student distribution was selected as having the best fit. Formal tests were conducted to validate this choice, using a significance threshold of p-value = 0.05.

3.1. Exponential Smoothing Algorithms

This research relied on simple, Holt and Holt–Winters exponential smoothing models. At first glance, the simple exponential smoothing method resembles the moving-average technique in that both approaches extract random behaviour from time series observations by smoothing historical data (see [27], chapter 10). The innovation introduced by the simple exponential smoothing model is that this method assigns different weights to each observation in the series. In the moving-average technique, the observations used to formulate future value forecasts contribute equally to calculating the predicted values. Simple exponential smoothing, in contrast, uses the latest information by applying a factor that determines its importance (see [28], chapter 4). According to [29], the argument in favour of the latter treatment of time series observations is that the latest observations are assumed to contain more information about futures and thus are more relevant to the forecast process.

The simple exponential smoothing model is represented as follows. If $\alpha \in \left(0,1\right),$ then Equations (1) and (2) can be written:

$$N_t=\alpha y_t+\alpha \left(1-\alpha \right)y_{t-1}+\alpha {(1-\alpha )}^2{y}_{t-2}+\alpha {(1-\alpha )}^3{y}_{t-3}+\dots$$

(1)

$$N_t=\alpha y_t+\left(1-\alpha \right)N_{t-1},$$

(2)

in which $N_t$ is the forecast at time $t,$ and $\alpha$ is the weight assigned to observation $y_t$. In the present model’s fit, $\alpha \mathrm{equals} 0.9999223,$ and the forecast coefficient is $150.6$. Figure 4 presents the fitted values of the simple exponential smoothing algorithm applied.

Figure 5 shows a ‘zoom in’ image of the data, with only the last 31 daily entries. This facilitated a comparison of the real dataset not used in the data modelling with the model’s predictions in blue. The dashed blue line is the confidence interval. The mean squared error is $4.982926$. This figure highlights the model’s forecasting inadequacies.

When the simple exponential smoothing method is applied to predicting time series that involve trends in past observations, the predicted values overestimate or underestimate the actual values [27]. Thus, the forecast’s accuracy is impaired. The linear exponential smoothing method was developed to avoid this systematic error, in an effort to reflect the trends present in data series [29].

The latter method obtains forecast values by applying Equation (3):

$${\widehat{y}}_t\left(h\right)=N_t+h{T}_t, h=1,2,\dots$$

(3)

in which $N_t$ corresponds to the forecast at time $t$. $T_t$ represents the trend component obtained with Equation (4), and $h$ is the forecast horizon:

$$N_t=\alpha y_t+\left(1-\alpha \right)\left(N_{t-1}+T_{t-1}\right)T_t=\beta \left(N_t-N_{t-1}\right)+\left(1-\beta \right)T_{t-1}$$

(4)

in which $\alpha$ is the weight attributed to the observation $y_t,0<\alpha <1$, and $\beta$ is the smoothing coefficient. The present study’s Holt exponential smoothing model has the values $\alpha \approx 1^{-}$ and $\beta =0.01184473$. The coefficients for the forecast are 150.6 and 0.1450552.

Figure 6 shows, in red, the fitted values for the Holt exponential smoothing model. This graph looked similar to Figure 4 above, so the next step was to zoom in on the data’s last part. Figure 7 shows the results with the simple exponential smoothing model’s values in blue and the Holt exponential smoothing model’s values in red. This graph highlights the differences between the latter two models.

Figure 8 presents the daily entries from 1980 to 2015 and, in red, the fitted values with the posterior predictions. The dashed lines are the confidence interval, and the mean squared error for the predictions is $14.36051$. This figure again shows the Holt model’s forecasting inadequacy, as the values are clearly overestimated.

The third exponential smoothing algorithm considered was the Holt–Winters approach. In addition to the unobservable components present in the previous model, this algorithm also incorporates a third—seasonality [30]. Since the analysed data did not present seasonality, the current study only used the simple and Holt exponential smoothing models. Of these two models, the best fit was the simple exponential smoothing option. Despite only giving fixed-value predictions, this model presented a lower mean squared error in the forecasts.

3.2. ARIMA

This section discusses the fit obtained for the futures prices using an ARIMA approach. The ARIMA model was proposed by [31]. The present research’s first model was the traditional ARIMA(1,1,0) and ARIMA(2,1,1), after which a regression variable was used, namely, an ARIMAX model. A seasonal ARIMA model was not included because, as mentioned previously, the data did not have a seasonal component.

The futures price series were differentiated once, and the autocorrelation function (ACF) was obtained as shown in Figure 9. This graph confirmed that one differentiation was enough, so $d=1$.

The first line in Figure 9 above represents the ARIMA(1,1,0) results, with an autoregressive coefficient equal to $0.0894$ and an Akaike information criterion (AIC)$\mathrm{of} 4830.15$. Regarding the diagnostic analysis of this model, the standardised residuals behaved as a white noise sequence, and they did not present either an autoregressive or moving-average component.

The second model shown in Figure 9 above is ARIMA(2,1,1), with autoregressive coefficients of $0.6272,0.0037$ and a moving-average coefficient of $-0.5445$. This model has an AIC $\mathrm{equal} \mathrm{to} 4826.42$, which is smaller than the previous model’s criterion. Regarding the diagnostic analysis of the ARIMA(2,1,1) model, the standardised residuals again behaved as a white noise sequence, and they did not present either an autoregressive or moving-average component.

In both models, the p-values for Ljung–Box Q-statistics are higher than 0.05. The futures prices were predicted for 16 days. Figure 10 shows in blue and red the ARIMA(1,1,0) and ARIMA(2,1,1) models’ results, respectively. Figure 11 zooms in on only these two models’ predictions and presents the futures prices’ real value. Neither prediction clearly captured the fall in values happening in the real data. The sum of the predictions’ squared errors is $4.761213$ and $5.205089$ for the ARIMA(1,1,0) and ARIMA(2,1,1) models, respectively.

The results support the conclusion that neither prediction and, consequently, neither model was a good fit for this time series. However, more information was available to use in the modelling process, that is, the spot prices. Thus, the next step was to estimate the ARIMAX model with the spot prices as a regression variable. Orders c(1,1,0) and c(2,1,1) were used, as previously, for the purposes of comparison, but this time—after using the auto.arima function in R—order c(0,0,5) was also included.

The ARIMAX(1,1,0) model has an autoregressive coefficient of $0.0048$ and a regression coefficient of $0.4416$. The ARIMAX(2,1,1) model produced autoregressive coefficients of 0.0681 and 0.0116, a moving-average coefficient of $-0.0630$ and a regression coefficient of 0.4388. The ARIMAX(0,0,5) model’s coefficients are $1.1211, 1.1261, 0.9690, 0.6602 \mathrm{and} 0.3271$ for the moving-average component, and $3.2594$ was obtained for the intercept and $0.9750$ for the regression.

The models’ AICs are 4711.76, 4715.48 and 5606.42, respectively. Figure 12 presents the predictions generated by the ARIMAX models, with a sum of squared errors of $6.090724$, $6.06713$ and 18.64785, respectively. Although the predictions looked somewhat better, the ARIMAX models presented a higher sum of squared errors compared to the ARIMA approach.

Overall, the ARIMA(1,1,0) model presented the lowest sum of squared errors for its predictions, that is, $4.761213$. This figure is not only the lowest as opposed to the ARIMA or ARIMAX models’ results but also compared with the exponential smoothing algorithms and ARIMAX(0,0,5). The results obtained by using the auto.arima function in R include the highest sum of squared errors for the predictions: $18.64785$. Another class of models was estimated to facilitate a retrospective comparison.

3.3. GARCH

The GARCH approach was originally proposed by [32]. In the present study, this model was developed using the ‘rugarch package’ written in R. As described in the previous section, the futures prices series were differentiated for modelling purposes, and, after the predictions were generated, the data were transformed back into their original form to compare the forecasts with the real data. Given the dynamic conditional variance, the GARCH model estimated was configured as GARCH(1,1), the mean model was ARFIMA(1,0,1), and the conditional distribution was normal. All model parameters and distributions were statistically significant with a p-value of 0.05.

The optimal parameters are $\mu =$0.0538910, $ar1=0.684777$, $ma1=-0.635453$, $\omega =0.012223$, ${\alpha }_1=0.080113$ and ${\beta }_1=0.903628$, in which $ar1$ and $ma1$ represent the autoregressive and moving-average coefficients. The AIC is 2.2184. Figure 13 shows that the lines are close to each other, and the sum of squared errors is $0.6774444$—the lowest value up to this point.

After comparing the GARCH, ARIMA or ARIMAX and exponential smoothing models, the GARCH(1,1) version was found to be the best fit for the data. To finalise the modelling process, one more type of model, GARMA, was estimated for comparison purposes.

3.4. GARMA

To apply the GARMA method, the most appropriate family of distributions needed to be determined, as mentioned above in Section 3. The t-Student distribution was selected because the weight of the tails matches the financial data analysed. Formal tests were conducted to confirm the suitability of this distribution choice, using a significance threshold of 0.05. The orders that best fitted the real figures were $c\left(1,0\right)$ and $c\left(2,1\right)$. The models with and without intercept were considered to ascertain the best fit.

In this study, the VaR at 5% was always used, as it corresponds to the lower limit of a 90% confidence interval in terms of the predictions, that is, a t-Student distribution quantile. A confidence interval of 90% means 5% in each tail of the distributions, and only the 5% lower weight was of interest in this research. Regarding the distribution, the t-Student family was used. Ref. [33] report that the data binding function for $\mu$ should be identity and, for $\sigma$ and $\nu ,$ logarithmic.

The GARMA model was developed with order $c\left(1,0\right)$, that is, a GAR(1) model with intercept. Every estimation has a p-value lower than 1%. The resulting model is represented by Equations (5) and (6):

$${\mu }_t={\nu }_t=97.011+0.3787x_t+{\tau }_t,$$

(5)

with

$${\tau }_t=0.9987 \mathrm{A}\left(y_{t-1}, x_{t-1}, 97.011, 0.3787\right),$$

(6)

in which $A$ is the function that represents the autoregressive term, with an associated parameter $\varphi =0.9987$.

The GARMA model with order $c\left(1,0\right)$ and without intercept was also a generalised autoregressive GAR(1), except that it had no intercept. All estimates have a $p$-value less than 1%. The model developed is shown in Equations (7) and (8):

$${\mu }_t={\nu }_t=0.3961x_t+{\tau }_t,$$

(7)

with

$${\tau }_t=\mathrm{A}\left(y_{t-1}, x_{t-1}, 0.3961\right),$$

(8)

in which A is the function that represents the autoregressive term, with an associated parameter $\varphi =1.00$.

This model was not just autoregressive, because it also included a moving-average component. The GARMA model used order $c\left(2,1\right)$, namely, the autoregressive part had order 2, and the moving average component had order 1. As mentioned previously, all this model’s estimations have a $p$-value below 1%. This model is expressed in Equations (9) and (10):

$${\mu }_t={\nu }_t=94.5611+0.3873x_t+{\tau }_t,$$

(9)

with

$${\tau }_t=\sum _{j=1}^2 \mathrm{A}\left(y_{t-j}, x_{t-j}, 94.5611, 0.3873\right)+0.1495 \mathrm{M}\left(y_{t-j}, {\mu }_{t-j}\right),$$

(10)

in which $A$ is the function that represents the autoregressive term, with the associated parameter $\underset{\_}{\varphi }\prime =0.8353 \mathrm{and} 0.1631$. $M$ represents the moving-average component that has as an associated parameter $\theta =0.1495$.

These three models’ results are presented in Figure 14; Figure 15. All three are similar, so they are not shown separately. These graphs focus only on the data’s last part to be able to show the models by zooming in closely enough. The straight blue line separates the data used for modelling from the part utilised to compare the real data with the models’ prediction.

Figure 15 above presents the models’ forecasts. Of the four options analysed, the ones with order $c\left(1,0\right)$ with and without intercept and the one with order $c\left(2,1\right)$ with intercept also provided similar results in terms of forecasts. The predictions overlap regarding both goodness of fit (see Figure 14 above) and forecasts (see Figure 15 above). In this graph, the three models’ futures price values cross the estimated VaR line.

The last model estimated also had an order of $c\left(2,1\right)$, without intercept. As in the previous cases, all this model’s estimations present a $p$-value lower than 1%. The model is shown in Equations (11) and (12):

$${\mu }_t={\nu }_t=1.0086x_t+{\tau }_t,$$

(11)

with

$${\tau }_t=\sum _{j=1}^2 \mathrm{A}\left(y_{t-j}, x_{t-j}, 1.0086\right)+0.1573 \mathrm{M}\left(y_{t-j}, {\mu }_{t-j}\right) ,$$

(12)

in which A is the function that represents the autoregressive term, with an associated parameter $\underset{\_}{\varphi }\prime =0.7512 \mathrm{and} 0.2013$. In addition, M represents the moving-average component, with an associated parameter $\theta =0.1573$.

Figure 16 presents the real data in black, and the blue line marks the separation between data used in the final model and data separated for comparison with the forecasts. The fitted values are in orange. This graph zooms in to show only the last 60 days in order to display the overlapping lines correctly.

Figure 17 displays the futures price forecasts. This model with order $c\left(2,1\right)$ and without intercept stands out as having a better fit. Notably, the graph reveals that the model’s forecast curve closely follows the real figures. In addition, the lower confidence interval, that is, the estimated VaR using 95% confidence, is always below the actual future prices at all points, which supports the conclusion that the estimated VaR is reliable. The analyses also included finding the mean square deviation of the four GARMA models estimated. The results are presented in

Table 2. Mean square deviation of models’ forecast.

Model	Mean Square Deviation	AIC
Simple exponential smoothing	4.9829	–
Holt exponential smoothing	14.3605	–
ARIMA(1,1,0)	4.7612	4830.15
ARIMA(2,1,1)	5.2051	4826.42
ARIMAX(1,1,0)	6.0907	4711.76
ARIMAX(2,1,1)	6.0671	4715.48
ARIMAX(0,0,5)	18.6478	5606.42
GARCH(1,1)	0.6774	2.2184
GARMA c(1,0) with intercept	1.7036	3965.59
GARMA c(1,0) without intercept	1.7036	3963.59
GARMA c(2,1) with intercept	1.7862	3950.68
GARMA c(2,1) without intercept	0.1284	3957.93

.

Table 2 highlights in bold that the forecast generated by the GARMA model with order $c\left(2,1\right)$ is considerably more adequate, and this model’s mean square deviation is the smallest. Thus, this model was chosen among the GARMA options as the best way to obtain the VaR of futures prices, as shown in bold in Table 2.

The conclusions drawn from the models’ fits and Table 2 are that the two exponential smoothing models did not have as good a fit as the GARMA models, and that these two models presented a higher mean square deviation in the forecasts. Therefore, the best model to estimate futures prices’ VaR is the GARMA with order $c\left(2,1\right)$ and without intercept.

The body of work in the domain of forecasting using GARMA models, especially with commodity prices, has recently been enriched with studies exploring novel hybrid techniques. A pivotal study by [12] brought forward the GARMAWLLWNN process, a hybrid model that synthesises the k-factor GARMA process, empirical wavelet transform and the local linear wavelet neural network (LLWNN) methods. When evaluated using data from the Polish electricity markets, their hybrid model showed enhanced predicting performance, outdoing both the dual generalised long-memory k-factor GARMA-G-GARCH model and the individual WLLWNN.

On the other hand, [13] proposed the k-factor GARMA-EWLLWNN model, an inventive hybrid approach linking the dual long-memory process (GARMA and G-GARCH), the empirical wavelet transform (EWT) and the local linear wavelet neural network (LLWNN). Their model, tested using data from the same Polish electricity markets, demonstrated superior forecasting capabilities in comparison to the GARMA-G-GARCH process and the hybrid EWLLWNN.

In our research, though we leveraged the GARMA model to forecast commodity prices, we found results that parallel the performance of the models proposed in both studies. Our forecasts, derived using the GARMA model, align with the strong predictive capacities observed in both these studies. It is imperative to note, however, that though these studies combine several techniques to enhance the prediction accuracy, our research provides a testament to the robustness and efficacy of the GARMA model in its capacity.

Moreover, though both aforementioned studies focus on the Polish electricity markets, our analysis extends to a broader range of data spanning from 1997 to 2015. This provides a comprehensive insight into the model’s utility across varying datasets and time frames.

Our results contribute to the growing literature that attests to the utility of the GARMA model and its hybrids in forecasting, further cementing its position as an invaluable tool for time series forecasting, particularly in the context of commodity prices.

4. Conclusions

This study sought to apply the VaR concept in order to analyse the behaviour of BGI futures contract price series using different methodologies, so that investors can make better decisions than those offered by the options currently available for effective risk management. Based on empirical research, the model with the best fit in terms of VaR estimations for BGI futures prices is the GARMA with order $c\left(2,1\right)$ and without intercept, which was presented in Section 3.4. This model predicted the BGI futures prices accurately enough that the confidence interval includes the real values.

Since the t-Student distribution was used, the VaR could be estimated from the quantiles of the $t$ distribution. As shown in Section 3.4, the model estimated behaved in the desired manner, that is, disposing of all points below the actual futures price values. Thus, this model can protect futures investors by computing the VaR as the maximum loss. In future studies, it would be beneficial for researchers to update the database to encompass more recent economic changes and trends. Furthermore, they could examine and analyse different commodities to confirm whether the proposed model can also be utilised to help protect investors against losses in other futures contracts.