Investigating seasonality, policy intervention and forecasting in the Indian gold futures market: a comparison based on modeling non-constant variance using two different methods

Introduction A futures contract is an agreement between two parties to transfer ownership of an underlying asset at a fixed time in the future. It is a standardized contract, so the trade takes place through an exchange. The underlying quantity is physical gold in the case of a gold futures contract. This type of financial contract acts as an insurance to the investor as it allows them to benefit from favorable price movements and acts as a hedge against systemic financial risk (Kou et al. 2019) during unfavorable future price movements (Choudhry et al. 2015; Iqbal 2017). There is evidence that metals play an important role in diversifying risk (Beckmann et al. 2015; Lean and Wong 2015; Daskalaki et al. 2017; Alkhazali and Zoubi 2020). Abstract

Forecasting futures prices is an integral component of a profitable futures trading strategy. In this respect, this paper uses one of the most popular methods for forecasting and econometric analysis, namely, ARIMA models (Box and Jenkins 1970;Al-Shaib 2006;Box et al. 2016;Guha and Bandyopadhyay 2016;Challa et al. 2018;Mallikarjuna and Rao 2019;Challa et al. 2020;Dong et al. 2020). The effects of seasonality are prevalent in the gold market (Ball et al. 1982;Ma 1986;Lucey and Tully 2006;Wang et al. 2019;Xiao and Maillebuau 2020) and are attributed to the Halloween effect (Bouman and Jacobsen 2002;Zhang and Jacobsen 2013;Pariyaprasert and Boonchuaymetta 2018;Schmidbauer and Rosch 2018), the calendar effect (Jain 2019), the festival season in India (Source: World Gold Council report), hours of daylight (Kamstra et al. 2003), and the autumn effect (Baur 2013).
Behavioural issues such as investor sentiment play an important role in gold prices (Aggarwal and Lucey 2007). Investment in gold has always been of great sentimental value to Indians due to cultural relevance and its objectification as a secure financial investment. Consequently, India ranks second largest in private gold holdings (Source: NITI Aayog report). The position of a major consumer is strengthened when the domestic economic need is linked with the international market for efficient price discovery and risk management. The presence of an established spot market provides members in the Indian value chain an opportunity to offset price risk through the futures market. The gold futures market allows investors to take positions by paying small initial margins, yielding higher returns, and the difference can be invested in risk-free government bonds. This additional incentive to invest in returns with zero risk, is possible when holding a portfolio containing the gold futures contract, unlike with equity investments and gold Exchange-Traded Funds (ETFs). The role of gold as a stabilizing force in the financial system during market shocks (Baur andMcDermott 2010, 2016) highlights the need to investigate the volatility properties of the asset. Understanding and explaining the behavior of commodity volatility is imperative as it is useful in designing optimal hedging strategies for derivatives such as futures and options (French et al. 1987;Chou 1988;Kocaarslan et al. 2017;Shakil et al. 2018;Mo et al. 2018). Analyzing the residuals after fitting an ARIMA model may suggest the need for volatility models (Ping et al. 2013;Todorova 2017) and volatility forecasting to aid investment decisions (Samouilhan and Shannon 2008;Thupayagale 2010;Kumar 2014;Mukherjee and Goswami 2017;Emenogu et al. 2020).
In India, trading in gold futures began in 2003 and the derivative contract was initially listed in the Multi Commodity Exchange of India (MCX) and launched under the commodity derivatives segment in the National Stock Exchange (NSE) and Bombay Stock Exchange (BSE) in 2018. Therefore, MCX accounts for 95% of the trade. The spot price on the MCX for the gold contracts is arrived at from polling a panel of representatives from the value chain of the physical market, with ex-polling conducted twice in a day (12.15 pm to 12.45 pm and 4.00 pm to 4.30 pm). The price is ex-Ahmedabad, which is based on the location of delivery, and on the lower tax structure, although the latter has become irrelevant after the imposition of the Goods and Service Tax (GST) in 2017. The price excludes GST and any other additional tax.
The trade volumes in gold futures fell to half and have remained stagnant in the last five years. This is due to the imposition of the Commodities Transaction Tax (CTT) of 0.01% in the 2013-2014 Union Budget. While the imposition of a transaction tax on the commodities futures market reduces speculative trading, it may not be able to reduce price volatility in futures or spot markets (Edwards 1993). In the Indian context, it is observed that the transaction tax will be a burden on the operators in the commodity market and act as a deterrent of entry to hedgers who may choose to shift to illegal platforms (Pavaskar and Ghosh 2008). The liquidity in the Indian commodity derivatives market started recovering in the first half of 2018-2019 from a steady decline during 2015-2018. A report published in the last quarter of FY 2019-2020 showed that returns from Indian gold funds were up 27% in the last one year, and up 11.93% in the last three years (Source: The Economic Times) This evidence suggests that investment in Indian gold investment products is a good portfolio diversification option and investment analysts and fund managers may be tempted to allocate more than the prudential 10%. Therefore, our study is aimed to assist fund managers by modelling the volatility using time series forecasting models used in the estimation of Value at Risk (VaR), which is a measure of risk of loss for investments.
From the above introduction about the Indian gold market, we see the importance of finding suitable methods of analyzing this time series and the impact of seasonality, intervention, volatility, and forecasting. Thus, Sect. 2 will introduce seasonal time series models, an intervention test, and methods for describing non-constant variances. Section 3 discusses forecast performance of different time series methods, and Sect. 4 provides concluding remarks.

Modeling methodology
In this paper, we use a time series approach to build a statistical model on a nonhomogeneous, nonstationary time series after taking a proper degree of differencing and variance stabilizing transformations. In our preliminary analysis, we plot the time series, autocorrelation function and partial autocorrelation function for the monthly data to observe the presence of seasonality. Based on the results of our preliminary analysis, we construct models to accommodate the effects of the seasonality and intervention present in the series. However, for forecasting comparison and inference, we use only those models whose parameters are statistically significant. Further, the error distribution for the proposed model with non-constant variance is also examined and the results are tabulated.
The objective of this paper is to compare methods used to forecast future prices and to suggest the optimal method which is useful for practitioners. To obtain meaningful results through prediction, modelling the historical evolution of the series through time is imperative. Therefore, we use a time series approach to incorporate the distinctive features of this data series into model building.
Let Z t be a sequence denoting the monthly time series under analysis. The data used is the closing price of the derivative contract titled Multi Commodity Exchange of India Gold Commodity Future Continuation 1 (RIC: MAUc1) and is obtained from Thomson Reuters Eikon platform for 14 years from January 2005 to January 2019. It is a time series of closing price considering the rollover of the contracts on the fifth day of the contract expiry month. The trading unit is 1 kg, and the quotation value is 10 g. The data has been converted into a monthly series. So, the time unit t, used in this analysis is months. The currency denomination used to represent the closing price is the Indian Rupee. The whole data set is from January 2005 (t = 1) to January 2019 (t = 169) as shown in Fig. 1. We will use the first 156 observations for modeling analysis and the remaining observations for forecast comparison.
From Fig. 1, we see that the monthly data has an increasing trend. However, we are interested in analyzing if the seasonal pattern observed in the gold spot market is evident in the gold futures market, and in non-stationarity of the price series. We plot the autocorrelation and partial autocorrelation functions of the seasonally differenced monthly series to confirm our claim on seasonality. The effect of seasonality is evident with seasonal period, s = 12, in the plots shown below in Fig. 2.
Subsequently, we carry out the well-known augmented Dickey-Fuller (ADF) test (Dickey et al. 1984) to investigate the presence of a unit root in a seasonal time series. The result is given in Table 1, with p-values less than 0.05 indicating both regular differencing and seasonal differencing are needed.
A good time series modelling method is the well-known Box-Jenkins multiplicative seasonal ARIMA model (p, d, q) × (P, D, Q) s where (i) φ p (B) and θ q (B) are the regular autoregressive and moving average polynomials of order p and q respectively, P B s and Q B s are the seasonal autoregressive  where µ is the mean of the stationary Z t ; and (iii) a t is assumed to be the sequence of independent normally distributed random variables with zero mean and constant variance σ 2 a . A simple way to find whether a series has a constant variance is to check the plot of residual squares after a preliminary model fitting. The following Fig. 3 of the residual squares after a preliminary ARIMA model fitting clearly shows the variance is not constant.
There are two approaches to solve the non-constant variance problem: The first approach is the power transformation.
introduced by Box and Cox (1964), where in practice the value of is often chosen from the set {1, 0.5, 0, − 0.5, − 1} as the one that yields a model with the minimum Akaike Information Criterion (AIC) and/or Bayesian Information Criterion (BIC).
The other approach is the one first proposed by (Bollersev 1986), where Eq. (1) is rewritten as and we have the more general error process where e t are i.i.d random variables with zero mean and variance one, independent of past realizations of n t−i , and the roots of (1 − φ 1 B − . . . − φ r B r ) = 0 are outside the unit circle. To guarantee σ 2 t > 0, we assume that θ 0 > 0 and that φ i and θ j are nonnegative. The model in (4) with the  (5) is known as the generalized autoregressive conditional heteroscedasticity (GARCH) model of order (r, s) and is denoted by GARCH (r, s) (Wei 2006) (Chapter 15). We will now illustrate the two approaches in the following subsections.

The power transformation approach
In this subsection, we apply the power transformation function given in Eq.
(2) to stabilize the non-constant error variance. Based on the minimum AIC and/or BIC given in Table 2 obtained from a SAS program, it suggests a square root transformation of the series for a constant variance.

A seasonal ARIMA model for the square root transformed series
After regular and seasonal differencing, the autocorrelation and partial autocorrelation plots of the transformed series shown in Fig. 4 suggests a seasonal ARIMA model of order (3, 1). The model can then be written as follows where (1) t ranges from 1 to 156, (2) s = 12, 24, 36, (3) B is the backshift operator, (4) a t is assumed to be a sequence of independent normally distributed random variables with zero mean and constant variance σ 2 a . The parameter estimates of the coefficients to the above fitted model are obtained. After removing the non-significant seasonal autoregressive coefficients associated with B 12 and B 24 , the obtained model is:

Check intervention on square root transformed series
Intervention analysis is a useful technique to study the effect of an external event that causes a level shift in the time series (Box and Tiao 1975). In our study, we will check the effect of the tax known as CTT on the series, which is represented as a step function and included in the model. In this section, the impact of CTT on the transformed series is examined.
The model with the intervention is: where (1) ω o represents the initial impact of CTT which is also felt during the period after introduction, (2) the effect of intervention represented by I t is expected to produce a step change from March 2013, which is the time of announcement of the tax and is denoted by t = 100 in our dataset and (3) a t is assumed to be a sequence of independent normally distributed random variables with zero mean and variance one. The estimation result is given below: with all parameters being significant except the one pertaining to the intervention variable and coefficients of the seasonal autoregressive operators B 12 , B 24 and σ 2 a = 7.221. The residual autocorrelations and partial autocorrelations indicate the adequacy of the model. From the estimates of model fitting, we see that the effect from intervention is (9) Z t = −1.5916 (2.5279) Fig. 5 The plot of autocorrelation and partial autocorrelation of residuals from model in Eq. (7) negative but not statistically significant for a significance level of 5% or for any of the commonly used significance levels, since its p-value = P Z < −1.5916 2.5279 = 0.2645, and the commonly used level at the final conclusion is 5%. Therefore, we omit the variable for tax intervention I t , and use Eq. (7) as our final chosen model for the square root transformed series.

An ARIMA/GARCH modelling approach
In this section, after fitting a time series model for the gold future prices, we will fit an ARCH/GARCH model to model its error variance. We refer readers to existing literature that uses similar methods of analysis and discusses the reliability of results obtained through econometric software (McCullough and Vinod 1999;McCullough and Renfro 2000).

A seasonal ARIMA model fitting
As indicated earlier, in order to determine whether the series has a constant variance, we preliminarily fitted an ARIMA (0, 1, 0) × (3,1,1) 12 model to the data set, which is based on its plot of autocorrelation and partial autocorrelation of the regularly and seasonally differenced series. The exact model form is given below: where (1) t ranges from 1 to 156, (2) B is the backshift operator, (3) n t is assumed to be a sequence of independent normally distributed random variables with zero mean and non-constant variance, and both regular and seasonal differencing is order 1 with seasonal period of 12.
The parameter estimates of the coefficients to the above fitted model are obtained and after removing the non-significant seasonal autoregressive coefficients associated with B 12 and B 24 , the obtained model is: The time dependent error variance will be discussed later.

A seasonal ARIMA model with intervention
In this subsection, we will check the effect of a CTT on the regular and seasonally differenced series. The model with the intervention is: where n t is assumed to be a sequence of independent normally distributed random variables with zero mean and non-constant variance and the definitions of other terms are the same as those in (8).
The parameter estimates of the coefficients are obtained and it is observed that the seasonal autoregressive coefficients associated with B 12 and B 24 are not significant and (10) hence they are omitted, and the parameters are re-estimated. The model in (12) with reestimated parameters is: with all parameters except the one pertaining to intervention variable being significant and the time dependent error variance will be discussed later. The residual autocorrelations and partial autocorrelations indicate the adequacy of the model.
From the estimates of model fitting, it is evident that the negative signs of the intervention parameter estimates show that the introduction of CTT reduced the gold futures prices but not quite significantly. The estimates of the intervention variable are not statistically significant, since its p value = P Z < −437.8838 737.7356 = 0.2764, for the commonly used 5% significance level. Therefore, we omit the intervention and simply use Eq. (11) in the further analysis.

Volatility models for nt
Here, we apply the general framework discussed at the beginning of this section to the squared residuals obtained from fitting the model given in Eq. (11). The autocorrelation and partial autocorrelation plots of the squared residuals of the model in Eq. (11) as shown in Fig. 6 suggest fitting a seasonal ARCH (2) model, with seasonal period s = 12.
Based on the highest significant spike at lag 24 in the plot of the partial autocorrelation function, a seasonal ARCH (2) model with seasonal period 12 is fitted on the squared residuals of the model in Eq. (11).
The parameter estimates of the seasonal ARCH (2) model are obtained and after removing the non-significant seasonal autoregressive coefficient associated with B 12 the obtained model is: where e t is assumed to be a sequence of independent normally distributed random variables with zero mean and variance one. Fig. 6 The plot of autocorrelation and partial autocorrelation of the squared residuals from model in Eq. (11) Nargunam et al. Financ Innov (2021) 7:62 So, for the gold futures prices, we have a seasonal ARIMA-GARCH model with its estimation results given below.

Forecasts from price models and their performance
The models from Eqs. (7) and (15) are used to forecast future values. The l-step ahead forecasts are calculated from the time origin t = 156 as follows.
The l-step ahead forecast equation for the model in Eq. (7) is where Y t = √ Z t and a t is the white noise series of independent N (0, σ 2 a ) random variables.
The l-step ahead forecast equation for the model in Eq. (15) is For comparison, the l-step ahead forecasts from Eq. (16) are converted to the degree of the original series by raising the values to the power 2. This conversion also raises the value of the standard deviation by the power 2 and is used in calculating forecast error variances and the associated forecast limits to be discussed in Sect. 3.2 below.
The predicted future price values are compared with the actual values of the closing price of the derivative and the forecast errors are calculated as the difference between the actual values and the forecast values as follows.
where n is the forecast origin, and l is the lead time of the forecasts from the same origin. Further, the forecast measures are calculated to ascertain the model which provides better forecasts. The comparison is usually based on the following summary statistics, mean square error (MSE) and mean absolute percentage error (MAPE), which are the most popular measures for forecast accuracy and consequently used in model selection. MSE is scale dependent and MAPE is scale independent.

Comparisons of forecast values
The l-step ahead forecasts of Model (7) through Eq. (16) and Model (15) through Eq. (17) are computed. For comparison, forecast values and forecast errors are tabulated in Table 3 below.
It is evident that between the fitted models in Eqs. (7) and (15), the forecasts from the model in Eq. (7) are closer to the actual values than those from the model in Eq. (15) in eleven out of thirteen cases. The forecast performance measures for the fitted models can further be illustrated using the summary statistics in Table 4.
The forecast measures between both models are coherent. The values of errors are lower for the model in Eq. (7) than the errors for the model in Eq. (15). This indicates a higher forecasting accuracy for model in Eq. (7), implying that it slightly outperforms the model in Eq. (15) in terms of forecast accuracy.

Comparisons of forecast limits
In time series modeling, we check the efficacy of our proposed model and its relevance to theory by forecasting. This, therefore, is an important application of time series analysis. In our study, we also want to investigate their associated forecast limits for the forecasts from both models.
When models in Eqs. (7) and (15) are used to compute forecasts using time series software like SAS or R, we can also calculate their associated forecast limits. For comparison, we have calculated the 95% forecast limits and the associated lengths of intervals  (7) and (15) Table 4 Forecast performance measures between models (7) and (15) Forecast performance measure Model in ( between upper and lower limits. The plots of the forecast limits along with predicted future values are given in Fig. 7, and the values are presented in Table 5.

Concluding remarks
This paper employs time series models to examine the gold futures market in India taking into account the seasonality prevalent in this market and the implementation of a new tax scheme, namely CTT. The estimation results show that the parsimonious forecasting models without intervention variable provide a better fit to the gold futures data. It is observed that the coefficient of the intervention variable is negative, but the results from the analysis show that the parameter estimates of the intervention variable are not statistically significant. Therefore, the forecasts carried out in this study do not include the intervention variable. The decrease in price trend observed in Fig. 1 after 2013, shows the lack of financial literacy among the market participants with regard to the new tax regime that offers an income tax benefit and that the market was dominated by speculative traders. For the futures market to work, both hedgers and speculators are needed, and both play a crucial role. This highlights the need for stringent monitoring to avoid manipulation in futures markets. It also calls for the need to organize  awareness programs advocating investment alternatives and their benefits. A recent empirical study on an emerging market shows how speculative capital movement is monitored to improve efficiency in the financial markets . Market sentiments also play a major role in the Indian context, mainly due to the sentimental attachment to the underlying asset. Since India is a price taker due to the absence of domestic mining the fluctuations in gold spot price are to a greater extent influenced by macroeconomic variables. The world gold price is derived from the London price fix. China follows a similar benchmarking methodology to deciding the price, known as the Shanghai price fix or the Shanghai gold benchmark. In India, the Indian Bullion Jewelers Association (IBJA) publishes daily price for gold. However, there is a lack of transparency due to the absence of a regulatory framework to effectuate the Indian price fix or benchmark price for gold in the domestic currency (Indian rupees) for as low a denomination as one gram or ten grams of gold and calls to attention the need for an independent organization to oversee commercial gold trade.
In Table 3, we observe that according to the model in Eq. (7) and the model in Eq. (15), the price forecasts are increasing for first nine months, but are decreasing for next four months based on model in (15), however, the results from model in (7) is increasing in lead time = 10, 11 and subsequently decreases. The forecasts from the model in Eq. (15) suggests that the Indian gold futures market does not work in tandem with the gold spot market, suggesting the presence of "reverse autumn effect" in the gold futures market. This shows that the anomaly present in the spot market is not likely to be found in the derivative market. Therefore, gold futures is an effective hedge and an attractive investment option to investors who might want to include gold in their portfolio to hedge against the volatility in the spot market in case of a long position and also earn short sale profits due to prices moving in opposite directions in different markets. There may also be short windows for arbitrage due to the delay in convergence of basis (difference between spot and futures prices) due to market frictions. It is observed that the trade in the gold futures market has increased over the first quarter in 2018.
Since all actual values are contained in the 95% forecast limits, the forecasts from the fitted models suggest that the future values of the gold futures contract can be predicted from historical values with a certain degree of accuracy. Consequently, they can be used by fund managers and active investors to re-adjust their portfolio allocations.
The method used for determining the power transformation was carried out using SAS and the rest of the analysis in our study was carried out using R Studio software. The code for the same is available upon request to the authors.