In [1] , Wets and Rios introduced a novel methodological approach to account for the distinct features of short term regimes via-a-vis the long term regimes. More specifically, a drift term was taken into consideration in case of the former, while recognizing that in case of the latter, the process is assumed to be drift-less, due to the volatile nature of the process of commodity prices in the long term, resulting in a mean-reverting process being ruled out for modeling of long term dynamics. Further, it was identified that the drivers of the short term process, in addition to historical data, also takes into consideration the aspects of state of the market as a whole. The proposed setup was then empirically analyzed in the context of copper prices, while recognizing the feasibility of the same in case of other commodities. The authors, recognizing that the transition from the short term to long term process, leading upto the "global" setup is non-trivial, carried out the same through empirical data analysis. Further, they cited existing literature to surmise that the long-term dynamics is triggered typically between the third and fourth year of the period under consideration. Accordingly, the global process was modeled using a linear combination of the short term and the long term processes. In our work, we adopt a mean reversal model for the long term process (after [1] ).
In [2] , Schwartz and Smith, developed a two factor model for commodity prices to capture two things, namely, the mean-reversion behavior of the short-term dynamics, and the uncertainly in the equilibrium levels of the prices, resulting from mean-reversion. These two factors were then collated to model the dynamics for the spot prices of the commodities. While estimation of the equilibrium price levels can be estimated from price movement of commodity futures with long maturity, the short-term behavior of the prices can be inferred from the price difference between short term and long term futures, and the model was validated using the prices of oil futures. In a detailed study [3] , Hilliard and Reis, investigated the impact of stochastic yields and interest rates, along with jumps in spot price, in the pricing of commodity derivatives. In particular, the authors modeled the spot prices as a jump-diffusion process, and also an equilibrium convenience yield diffusion. A two factor model was used for valuation of futures and options. This was then extended for three factors, as well as three factor jump-diffusion valuation model. In the context of our work, one of the assumptions in [3] , is that the net marginal rate of yield of the commodity follows a mean-reverting process. In [4] , the authors presented logarithmic mean reversion models for commodity pricing. Schwartz [5] , presented a comparative analysis for three stochastic models for commodity prices, by considering a mean-reversion process, in the paradigm of predictive modeling of futures prices, as well as their valuation vis-a-vis, other derivatives and assets. The comprehensive analysis undertaken in the paper, gave credence to the premise that there exists a strong mean reversion factor when it comes to pricing of commodities. Hart et al. [6] attributed the lack of existence of options on agricultural commodities, with maturity exceeding a year, to the lack of option pricing mechanism, in case of commodities, whose spot prices can be encapsulated by a mean-reversion process. Accordingly, they proposed and analyzed a model to account for the mean-reversion in the spot prices, and the model also corrected for seasonal influences, motivated by the empirical observance of both mean-reversion and seasonal influence, in case of price dynamics in soybean markets. In [7] , Bernard et al. carried out a detailed comparative analysis of several approaches, namely, GARCH, jump-diffusion and mean reversion, to forecast commodity prices. The empirical analysis carried out making use of spot prices of aluminum and futures price series, produced mixed results, that is, no clear uniform dominance of any model, could be established.
It is important to take note of the fact that all the above models cited from literature are Markovian, wherein the instantaneous drift term is a function of the instantaneous or spot price. However, in the context of our discussion is would be more realistic to account for the fact that, in the case of metals or other heavy industries, the industry takes time to adjust to the market information, namely, demand and supply. This motivated the consideration of a time lag factor, i.e., the instantaneous drift term will be a function of the price of the commodity at some earlier time point. Accordingly, we propose our Stochastic Differential Equation (SDE), to model the instantaneous price x(t) of the commodity, at time t, as,
is a Wiener process. The two key aspects of our model is that unlike the other models, this SDE can neither be analytically solved, nor is it Markovian. In fact, it is a Higher Order Markovian (HOM) model of order τ [8] .
Note that the Markovian analogue of equation (2.1), is given [1] by
Economic theory suggests that, if there are changes in the prices of a commodity due to the volatility of the markets, the industry will respond to that accordingly. For instance if there's a sharp rise in prices, the profitability of the commodity, increases, which gives the industry an incentive to increase the supply, which subsequently reduces the prices. Similarly, if there is a fall in the prices, the profitability reduces and it is likely that the industry will reduce supply, and thus prices would increase. However for most of the industries, like that of base metals, it is unlikely that producers will bring changes in production immediately, as a response to the market information. It is more likely that the industry will respond with a delay to the market information. This response time (τ ) serves as the order of our HOM process. Another significant advantage of a HOM model is that the speed of inversion is significantly lower than a Markovian model. This implies that the price can drop significantly below the mean. In fact, it also allows for negative prices of commodities (an example of which is, when the WTI May contract was in the negative zone, on 20th April, 2020).
We determine the long term mean of the SDE described in equation (2.1) . Let x(t) be the mean of x(t) as t → ∞ and let ∆t be a small interval, in which the price of the commodity increases by ∆x, in that interval. Therefore,
But, since t >> ∆t, therefore,
Taking the expectation of equation (2.1), we obtain,
At this point, we note that this result is useful while estimating parameters.
For the numerical analysis of the Copper Mini prices, we assume a(t), b(t) and σ(t) to be constants. However, over longer periods, more realistic estimates can be made by treating them as functions of time, like accounting for inflation in b. The parameter estimation was executed using the Euler-Maruyama Maximum Likelihood Analysis [9, 10] . We can use Maximum-Likelihood estimation as a HOM Process, which is the solution of the Fokker-Planck equation, and is stated as the following Theorem, the Proof of which is given in Appendix A.
Theorem 1. For a Higher Order Markov process z, of order τ , that takes initial values of the function f defined over an interval [t 0 , t 0 + τ ], z satisfies the following equation
In order to estimate the parameters of the model, we first exclude a fraction of the historical data to accommodate the time lag proposed in our model. After this, we set aside the initial 80% of the remaining data points for the training part of the implementation, with the final 20% being used in the validation part of the implementation. Initially, an approximate estimation of the model parameters a and τ is obtained using a trial and error approach. Further, the estimation of b is taken to be the mean of the spot prices of the training data set, from equation (2.4) . Finally, in order to estimate the value of σ, we start with an initial guess of σ = Std.Dev (ln(X)), where the random variable X takes the values of the spot prices of the training data set. After we have obtained an approximate estimate of the four parameters (as explained in the preceding discussion), we then worked on obtaining a more accurate estimation of the parameters, through the implementation of a coordinate descent algorithm. At this point, we note that instead of a coordinate descent of the loss function, we employ a coordinate ascent of the log-likelihood function, which is simply a maximization algorithm [11] . In the final step, we then forecast the prices for the remaining 20% of the data points, time by simulations (using Monte-Carlo simulation) and then average them out to obtain a mean path, and then compare the results obtained using our model, to that of the Markovian models. The full algorithm for parameter estimation and forecasting is given in Appendix B.
The primary focus of this article is the Copper Mini prices in India, in addition of a few other base metals, namely, Aluminum, Lead, Nickel and Zinc. The emphasis on Copper Mini prices was due to the availability of a reasonably sized data set. For the remaining metals, the choice was based on the liquidity and the availability of data, albeit to a much lesser level than Copper Mini.
For the purpose of empirical analysis reported in this paper, we made use of the available spot prices from the website of MCX-India [12] . We focus our analysis on Copper Mini for the following tabulated period ( Table  1) . As already noted, the approach to estimation of data, first involved the exclusion of a fraction of the historical 
The errors in the forecasting of Copper Mini prices for both the HOM and the Markovian models in the paradigm of the enumerated approaches for the estimation of errors is presented in Table 3 , wherein, it can be observed that, in all the five approaches, the errors in estimation is consistently and significantly lower in case of the HOM models. We also plot graphs to compare the actual price of the commodity with the mean forecast of the HOM and the Markovian models, as can be seen in Figure 1 .
We also analyze the probability distribution of the spot prices at a given time point. Accordingly, we test the null hypothesis that the probability distribution of the spot prices at a given time is log-normal. We use the Kolmogorov-Smirnoff test and the Anderson Darling test on an ensemble of 2000 simulations of the Copper Mini spot price using the HOM model for different time steps, which are presented in Table 4 and 5, respectively. We can see that the null hypothesis cannot be rejected for the Kolmogorov-Smirnoff test or the Anderson-Darling test for any of the time intervals. We also plot the log-distributions at the t = 210 days (the choice being for illustrative purpose) in Figure 2 . 
We extend our analysis to four other base metals traded on MCX-India [12] , namely, Aluminium, Lead, Nickel and Zinc, on their spot prices in the following ranges of dates, as detailed in Table 6 . Since the size
Date range Number of Data Points Aluminum 01-03-2019 to 07-08-2020 652 Lead 03-06-2019 to 07-08-2020 517 Nickel 03-06-2019 to 07-08-2020 532 Zinc 01-04-2019 to 07-08-2020 614 Table 6 : Details of data of other base metals.
of the available data for these metals is significantly smaller than that of Copper Mini, we leave alone 75 data points for historical information in the parameter estimation. The rest of the parameter estimation and forecasting methodology remains the same. The estimated values of τ are given in For the estimation of errors, we have excluded the period wherein the prices of the commodities were frozen as the market was not operational due to the lock-down. We can see that the HOM model gives better estimations of the mean price of the path, than the Markov model, in all cases, except that of Aluminum. The error in Aluminum may have been caused by the increased volatility in Asian markets due to the COVID-19 pandemic that could not be accounted for as our parameters were estimated in the pre COVID-19 period [14] . It is also likely that the estimation would improve if we could have access to more data and estimate a better response time as that would allow for the estimation of a bigger time, which we have been unable to do due to unavailability of data. However, it is interesting to note that the models make accurate predictions about the price of the commodity by August 2020, by the time the market had reopened. The rise of the price of Aluminum above the model prediction is also probably a consequence of the lock-down as well, and can be explained by pent up demand in the market following the lock-down.
The main contribution of this work is the presentation of our approach to the modeling of commodity prices, using the Higher Order Markovian method, as an improved alternative to the Markovian model. We show, using historical data of spot prices, that one achieves a better prediction of the spot prices of base metals in the Indian market (in all cases studied except that of Aluminium), by factoring in the response time of the industry to respond to market information, which is the order of the Higher Order Markov Process. Our analysis revealed three advantages, namely, encapsulating the delay in the response of various base metal industries to the market information, improving the predicted rate of long-term mean reversion, and reducing the error of estimation, therefore providing better predicting power as compared to the Markovian model. Proof.
For a τ -th order Markov process,
where y and z are the state of the system at time T and t respectively. Further, f is a function defined on R over (t 0 , t 0 + τ ) such that t 0 + τ < t < T and P (b, a) is the conditional probability of event b occurring, given a has occurred. Assuming time homogeneity we will omit the notation of time with the conditioning variable i.e., we will write P (y, t + ∆t|f, t) = P (y, ∆t|f ) where ∆t is the difference in time of y and the upper limit over which the interval f , is defined over. The Fokker-Planck equation can be derived using this (after the derivation in [8] ).
where V (z) = D (1) (z) and D(z) = D (2) (z) with D (n) (z) = 1 n! lim 