Analytical Pricing of Commodity Futures with Correlated Jumps and Seasonal Effects: An Empirical Study of Thailand’s Natural Rubber Market

Abstract

This paper presents a novel multivariate mean-reverting jump-diffusion model that incorporates correlated jumps and seasonal effects to capture the complex dynamics of commodity prices. The model also accounts for the interplay between price volatility and convenience yield, offering a comprehensive framework for commodity futures pricing. By leveraging the Feynman–Kac theorem, we derive a partial integro-differential equation for the conditional moment generating function of the log price, enabling an analytical solution for pricing commodity futures. This solution is validated against Monte Carlo simulations, demonstrating high accuracy and computational efficiency. The model is empirically applied to historical futures prices of natural rubber from the Thailand Futures Exchange. Key parameters—including commodity price dynamics, convenience yields, and seasonal factors—are estimated, revealing the critical role of jumps and seasonality in influencing market behavior. Notably, our findings show that convenience yields are negative, reflecting higher inventory costs, and tend to increase with rising spot prices. These results provide actionable insights for traders, risk managers, and policymakers in commodity markets, emphasizing the importance of correlated jumps and seasonal patterns in pricing and risk assessment.

1. Introduction

The modeling of commodity prices is a cornerstone in financial economics, underpinning the valuation of derivatives, risk management strategies, and investment decisions. Commodity prices exhibit unique characteristics such as mean reversion, seasonality, stochastic convenience yields, and abrupt price jumps reflective of market shocks. Traditional stochastic models that assume smooth price paths often overlook these complexities, introducing estimation biases that can distort risk assessments and pricing accuracy. In response, this paper introduces the multivariate mean-reverting jump-diffusion (MMRJD) model, which incorporates both correlated jumps and seasonal effects to better capture the true dynamics of commodity markets.

The theory of storage provides a fundamental framework for understanding commodity price behavior by emphasizing the interplay between spot and futures prices, convenience yields, and inventory levels. Pioneering works by Kaldor [1], Working [2], and Brennan [3] laid the foundation for this theory, linking convenience yields to market conditions such as backwardation and contango. Subsequent stochastic models, beginning with the geometric Brownian motion approach of Samuelson [4] and later advanced by Schwartz [5], have been instrumental in jointly describing commodity prices and convenience yields. However, these models often fall short in addressing abrupt market events and cyclical behaviors.

Seasonality is a critical feature in many commodity markets, particularly in sectors like agriculture and energy, where supply and demand are influenced by weather patterns, harvest cycles, or production schedules. Researchers such as Cartea and Figueroa [6] and Borovkova and Geman [7] have incorporated seasonal components into pricing frameworks, while more recent studies like Frau and Fanelli [8] highlight its role in the long-run equilibrium of mean-reverting processes. Ignoring seasonal effects can lead to systematic misestimations of volatility and disrupt the accurate forecasting of price dynamics, ultimately impairing both pricing and risk management.

Abrupt market events—ranging from geopolitical crises to natural disasters—can induce sudden price jumps that traditional diffusion models fail to capture. Merton [9] introduced the jump-diffusion model to address discontinuities in price trajectories, and later works by Hilliard and Reis [10] and Cartea and Figueroa [6] extended these ideas to commodity markets. More recent studies by Fileccia and Sgarra [11], Nguyen and Prokopczuk [12], and Brignone et al. [13] emphasize the importance of modeling correlated jumps between commodity prices and convenience yields, noting that neglecting these phenomena can lead to significant estimation biases and ineffective hedging strategies during periods of market stress.

Building on these advancements, the proposed MMRJD model extends classical two-factor frameworks by integrating correlated jumps and a sinusoidal seasonal component into the joint dynamics of commodity prices and convenience yields. This dual approach effectively captures both the abrupt changes driven by rare market events and the periodic fluctuations inherent to seasonal patterns. By explicitly addressing these aspects, the model mitigates estimation biases in key parameters and provides a more comprehensive depiction of market behavior, enhancing its utility for forecasting and risk management.

The significance of the MMRJD model lies in its strong technical and economic justifications. Technically, by incorporating both seasonal effects and correlated jumps, the model overcomes the limitations of earlier frameworks and reduces the risk of parameter estimation bias. Economically, a more accurate representation of commodity price dynamics leads to better pricing of derivatives, more effective hedging strategies, and improved overall risk management. This is particularly relevant in markets such as Thailand’s natural rubber market, where the interplay of sudden shocks and pronounced seasonal trends demands a sophisticated modeling approach. In essence, the innovations presented in this paper offer both practical relevance and enhanced performance in real-world applications, paving the way for more robust financial strategies and informed decision making.

The primary contributions of this paper are as follows:

1. Model Development: The MMRJD model advances classical frameworks by incorporating correlated jumps and seasonal effects, capturing critical market dynamics such as abrupt price changes and periodic fluctuations.

2. Closed-Form Solution: By applying the Feynman–Kac theorem, a partial integro-differential equation for the conditional moment generating function (CMGF) of the log price is derived. This enables an analytical solution for pricing commodity futures, significantly enhancing computational efficiency compared to simulation-based methods.

3. Empirical Validation: The model is validated using historical futures data from Thailand’s natural rubber market (3 January 2019 to 27 December 2023), demonstrating its effectiveness in capturing market dynamics and outperforming existing models.

The remainder of this paper is organized as follows. Section 2 introduces the MMRJD model, emphasizing its integration of correlated jumps and seasonal effects to model the complex dynamics of commodity prices. Section 3 derives closed-form pricing formulas for commodity futures, leveraging the analytical framework of the model. Section 4 validates the model through numerical simulations and an empirical analysis of Thailand’s natural rubber market. Finally, Section 5 concludes the paper.

2. Multivariate Mean-Reverting Jump-Diffusion Model: Correlated Jumps and Seasonal Effects

This section introduces the MMRJD model, which extends the classical two-factor framework proposed by Schwartz [5]. The MMRJD model incorporates two critical features essential for capturing the dynamics of commodity markets: correlated jumps [9,10,11,12,13,14,15] and seasonal effects [6,7,8,16]. By integrating these elements, the model provides a comprehensive framework for describing the joint evolution of the commodity price $S_t$ and the instantaneous convenience yield $Y_t$. These enhancements enable the model to capture key market features, such as simultaneous jumps in prices and convenience yields, as well as periodic fluctuations driven by seasonality.

Under the risk-neutral measure $\mathbb{Q}$ with a filtration ${\mathcal{F}}_t$, the dynamics of $S_t$ and $Y_t$ are governed by the following stochastic differential equations:

(1)$\begin{array}{c}\begin{array}{cc}\hfill {\displaystyle \frac{d{S}_t}{S_t}}& =\left(r-\lambda \left(e^{\alpha +{\displaystyle \frac{{\beta }^2}{2}}}-1\right)-Y_t\right)dt+\sigma d{W}_t^{\left(1\right)}+\left(e^{J_t^{\left(X\right)}}-1\right)d{N}_t,\hfill \\ \hfill d{Y}_t& =\kappa \left(\overline{y}\left(t\right)-{\displaystyle \frac{\eta +\lambda \gamma }{\kappa }}-Y_t\right)dt+\xi d{W}_t^{\left(2\right)}+J_t^{\left(Y\right)}d{N}_t.\hfill \end{array}\end{array}$

Here, $W_t^{\left(1\right)}$ and $W_t^{\left(2\right)}$ are correlated Wiener processes with a constant correlation coefficient $\rho \in (-1,1)$, and $N_t$ is a Poisson jump process with intensity rate $\lambda >0$. The parameters $r>0$, $\sigma >0$, and $\xi >0$ represent the risk-free interest rate, the volatility of the commodity price, and the volatility of the convenience yield, respectively. The convenience yield reverts to its long-run equilibrium $\overline{y}\left(t\right)$ at a rate $\kappa >0$, while $\eta$ denotes the market price of convenience yield risk. The inclusion of the jump term introduces rare, abrupt changes in both the price and convenience yield, allowing the model to capture market shocks effectively.

The long-run equilibrium of the convenience yield, $\overline{y}\left(t\right)$, is modeled as a sinusoidal function to capture periodic fluctuations observed in many commodity markets:

(2)$\overline{y}\left(t\right)={\overline{y}}_0+{\overline{y}}_{1}sin(2\pi t+\varphi ),$

where ${\overline{y}}_0$ is the stationary level, ${\overline{y}}_1>0$ governs the amplitude of the seasonal pattern, and $\varphi \in (0,2\pi )$ determines the phase shift of the cycle. This functional form reflects the seasonality inherent in markets such as agriculture and energy, where supply and demand exhibit periodic behavior driven by factors like weather or production schedules.

The MMRJD model (1) also accounts for correlated jumps in commodity prices and convenience yields. The magnitudes of these jumps, denoted $J_t^{\left(X\right)}$ and $J_t^{\left(Y\right)}$, are assumed to follow a bivariate normal distribution:

(3)$\begin{array}{c}g(\acute{x},\acute{y})={\displaystyle \frac{1}{2\pi \beta \delta \sqrt{1-{\nu }^2}}}exp\left(-{\displaystyle \frac{1}{2(1-{\nu }^2)}}\left[{\displaystyle \frac{{(\acute{x}-\alpha )}^2}{{\beta }^2}}+{\displaystyle \frac{{(\acute{y}-\gamma )}^2}{{\delta }^2}}-2\nu {\displaystyle \frac{(\acute{x}-\alpha )(\acute{y}-\gamma )}{\beta \delta }}\right]\right),\end{array}$

where $\alpha$ and $\gamma$ are the means, $\beta >0$ and $\delta >0$ are the standard deviations, and $\nu \in (-1,1)$ is the correlation coefficient between $J_t^{\left(X\right)}$ and $J_t^{\left(Y\right)}$. This structure captures the simultaneous occurrence of abrupt changes in prices and convenience yields, which is a common feature during periods of market stress.

To facilitate analytical tractability, the commodity price process is transformed into its logarithmic form, $X_t=ln{S}_t$. The resulting system can be expressed as

(4)$\begin{array}{c}\begin{array}{cc}\hfill d{X}_t& =\left(\mu -Y_t\right)dt+\sigma d{W}_t^{\left(1\right)}+J_t^{\left(X\right)}d{N}_t,\hfill \\ \hfill d{Y}_t& =\kappa \left(\vartheta \left(t\right)-Y_t\right)dt+\xi d{W}_t^{\left(2\right)}+J_t^{\left(Y\right)}d{N}_t,\hfill \end{array}\end{array}$

where the adjusted drift term is given by $\mu =r-{\displaystyle \frac{{\sigma }^2}{2}}-\lambda \left(e^{\alpha +{\displaystyle \frac{{\beta }^2}{2}}}-1\right)$, and the modified long-run equilibrium is

(5)$\vartheta \left(t\right)={\overline{y}}_0+{\overline{y}}_{1}sin(2\pi t+\varphi )-{\displaystyle \frac{\eta +\lambda \gamma }{\kappa }}.$

This representation retains the jump-diffusion structure of the original model while providing a clearer framework for deriving closed-form solutions for derivative pricing.

The MMRJD model (1) generalizes the classical two-factor framework by Schwartz [5]. When the jump intensity $\lambda$ and the seasonal amplitude ${\overline{y}}_1$ approach zero, the MMRJD model (1) reduces to the simpler two-factor setup. This flexibility allows the model to address a wide range of market conditions from purely mean-reverting dynamics to markets characterized by significant seasonal patterns and abrupt jumps. By integrating correlated jumps and seasonality, the MMRJD model (1) captures critical features of commodity markets, making it a powerful tool for pricing derivatives, evaluating risk, and analyzing market dynamics.

3. Analytical Pricing of Commodity Futures

The MMRJD model presented in (4) effectively captures both continuous and discontinuous dynamics in commodity prices and convenience yields, making it suitable for commodities subject to seasonal supply–demand fluctuations and rare events like market shocks. This formulation also simplifies the computation of various statistical properties such as the CMGF of the log price process, which is crucial for pricing derivatives and managing risk in commodity markets [17].

Given $T>0$, the CMGF of the log price process $X_T$ given $X_t=x$ and $Y_t=y$ can be written as

(6)$M(t,x,y;u,T)={\mathbb{E}}^{\mathbb{Q}}\left[e^{u{X}_T}|X_t=x,Y_t=y\right]$

where ${\mathbb{E}}^{\mathbb{Q}}\left[e^{u{X}_T}|X_t=x,Y_t=y\right]<\infty$ for $(t,x,y)\in [0,T]\times {\mathbb{R}}^2$. We denote ${\mathbb{E}}^{\mathbb{Q}}$ as the condition expectation under the risk-neutral measure $\mathbb{Q}$.

3.1. A Closed-Form Formula for the CMGF of the Log Price Process

We now present the closed-form formula for the CMGF (6), which decomposes the log price process into its diffusion, jump, and seasonal components. This decomposition enhances interpretability and aids in understanding the distinct roles played by continuous fluctuations, discontinuous jumps, and periodic trends in determining the overall behavior of commodity prices.

Theorem 1. 

Let $X_T$ denote a stochastic process driven by the MMRJD model (4). The CMGF of $X_T$, denoted by $M(t,x,y;u,T)$, can be decomposed into three distinct components: a diffusion component, a jump component, and a seasonal component. The CMGF can be expressed as

(7)$M(t,x,y;u,T)=e^{ux-u\left({\displaystyle \frac{1-e^{-\kappa (T-t)}}{\kappa }}\right)y+\mathrm{\Psi }(t;u,T)+\mathrm{{\rm Y}}(t;u,T){\overline{y}}_1+\mathrm{\Phi }(t;u,T)\lambda }$

for $(t,x,y)\in [0,T]\times {\mathbb{R}}^2$. The function $\mathrm{\Psi }(t;u,T)$ represents the diffusion component and is given by

(8)$\begin{array}{cc}\hfill \mathrm{\Psi }(t;u,T)& =u\left(r-{\displaystyle \frac{{\sigma }^2}{2}}-{\overline{y}}_0+{\displaystyle \frac{\eta }{\kappa }}\right)(T-t)\hfill \\ & +u\left({\overline{y}}_0-{\displaystyle \frac{\eta }{\kappa }}\right)\left({\displaystyle \frac{1-e^{-\kappa (T-t)}}{\kappa }}\right)+u^2\left({\displaystyle \frac{{\sigma }^2}{2}}-{\displaystyle \frac{\rho \sigma \xi }{\kappa }}+{\displaystyle \frac{{\xi }^2}{2{\kappa }^2}}\right)(T-t)\hfill \\ & +u^2\left({\displaystyle \frac{\rho \sigma \xi }{\kappa }}-{\displaystyle \frac{{\xi }^2}{2{\kappa }^2}}\right)\left({\displaystyle \frac{1-e^{-\kappa (T-t)}}{\kappa }}\right)-{\displaystyle \frac{u^2{\xi }^2}{4\kappa }}{\left({\displaystyle \frac{1-e^{-\kappa (T-t)}}{\kappa }}\right)}^2.\hfill \end{array}$

The function $\mathrm{{\rm Y}}(t;u,T)$ describes the seasonal component and is given by

(9)$\begin{array}{c}\mathrm{{\rm Y}}(t;u,T)=u\left({\displaystyle \frac{cos(2\pi T+\varphi )-cos(2\pi t+\varphi )}{2\pi }}\right)\hfill \\ -2\pi u\left({\displaystyle \frac{cos(2\pi T+\varphi )-cos(2\pi t+\varphi )e^{-\kappa (T-t)}}{4{\pi }^2+{\kappa }^2}}\right)\\ \hfill +\kappa u\left({\displaystyle \frac{sin(2\pi T+\varphi )-sin(2\pi t+\varphi )e^{-\kappa (T-t)}}{4{\pi }^2+{\kappa }^2}}\right).\end{array}$

The jump component, represented by $\mathrm{\Phi }(t;u,T)$, is given by

(10)$\begin{array}{c}\mathrm{\Phi }(t;u,T)=-u\left(e^{\alpha +{\beta }^2/2}-1-{\displaystyle \frac{\gamma }{\kappa }}\right)(T-t)-{\displaystyle \frac{u\gamma }{\kappa }}\left({\displaystyle \frac{1-e^{-\kappa (T-t)}}{\kappa }}\right)\hfill \\ \hfill +{\int }_0^{T-t}{e}^{u\left(\alpha -{\displaystyle \frac{\gamma }{\kappa }}(1-e^{-\kappa s})\right)+u^2\left({\displaystyle \frac{{\beta }^2}{2}}-{\displaystyle \frac{\nu \beta \delta }{\kappa }}(1-e^{-\kappa s})+{\displaystyle \frac{{\delta }^2}{2{\kappa }^2}}{(1-e^{-\kappa s})}^2\right)}ds-T+t.\end{array}$

Thus, the CMGF $M(t,x,y;u,T)$ provides a comprehensive representation of the contributions from diffusion, seasonal variations, and jumps within the model.

Proof. 

See Appendix A.1. □

Theorem 1 presents a comprehensive formula for the CMGF of the log price process (6) under the MMRJD model (4). The decomposition into diffusion, seasonal, and jump components enhances the interpretability of the log price process by separating continuous price fluctuations, rare jumps, and periodic trends. This decomposition allows for a more detailed understanding of how these elements interact to influence commodity prices. The explicit form of CMGF is instrumental in applications such as pricing derivative contracts, managing risk, and analyzing market shocks. Its practical significance lies in its ability to provide clear insights into the dynamics of commodities subject to both predictable and unpredictable forces.

3.2. Analytical Pricing of Commodity Futures

Commodity futures are financial instruments designed to standardize the terms of a transaction between two parties who agree to buy or sell a commodity at a predetermined price on a specified date in the future [18]. These contracts serve two main purposes: hedging and speculation. Hedging involves reducing exposure to price fluctuations, offering producers or consumers of commodities a way to lock in future prices, thereby stabilizing cash flows and managing risk. On the other hand, speculation allows traders to take positions in commodity markets, anticipating price changes to make a profit.

3.2.1. A Closed-Form Formula for Pricing Futures Contracts

Suppose that a futures contract is signed at inception time $t\in [0,T]$ with the transaction and delivery occurring at expiration time T. Following the no-arbitrage assumption, the futures price of a commodity is defined as the conditional expectation of the commodity price at expiration time T given the current information at time t [19]. Using the logarithm of the commodity price, $X_t=ln{S}_t$, the futures price $F(t,s,y;T)$ can be computed directly from the CMGF of the log price process:

(11)$\begin{array}{c}\hfill F(t,s,y;T)={\mathbb{E}}^{\mathbb{Q}}\left[S_T\mid S_t=s,Y_t=y\right]=M(t,lns,y;1,T),\end{array}$

where $M(t,lns,y;u,T)$ denotes the CMGF evaluated at $u=1$. This relationship establishes the futures price as a direct function of the CMGF, providing a practical framework for pricing.

Theorem 1 establishes the explicit form of the CMGF, which plays a crucial role in deriving closed-form expressions for futures prices. Proposition 1 demonstrates the analytical computation of the futures price using the CMGF and highlights its dependence on underlying parameters and model components.

Proposition 1. 

Given $T>0$, the futures price of a commodity can be expressed as

(12) $F(t,s,y;T)=s{e}^{f_0(t;T)+f_1(t;T)y+f_2(t;T){\overline{y}}_1+f_3(t;T)\lambda },$

where the functions $f_0(t;T)$, $f_1(t;T)$, $f_2(t;T)$, and $f_3(t;T)$ are defined as follows:

(13) $\begin{array}{cc}\hfill f_0(t;T)& =\left(r-{\overline{y}}_0+{\displaystyle \frac{\eta }{\kappa }}-{\displaystyle \frac{\rho \sigma \xi }{\kappa }}+{\displaystyle \frac{{\xi }^2}{2{\kappa }^2}}\right)(T-t)\hfill \\ \hfill & +\left({\overline{y}}_0-{\displaystyle \frac{\eta }{\kappa }}+{\displaystyle \frac{\rho \sigma \xi }{\kappa }}-{\displaystyle \frac{{\xi }^2}{2{\kappa }^2}}\right){\displaystyle \frac{1-e^{-\kappa (T-t)}}{\kappa }}-{\displaystyle \frac{{\xi }^2}{4\kappa }}{\left({\displaystyle \frac{1-e^{-\kappa (T-t)}}{\kappa }}\right)}^2,\hfill \end{array}$

(14) $\begin{array}{cc}\hfill f_1(t;T)& =-{\displaystyle \frac{1-e^{-\kappa (T-t)}}{\kappa }},\hfill \end{array}$

(15) $\begin{array}{cc}\hfill f_2(t;T)& ={\displaystyle \frac{cos(2\pi T+\varphi )-cos(2\pi t+\varphi )}{2\pi }}\hfill \\ \hfill & -2\pi {\displaystyle \frac{cos(2\pi T+\varphi )-cos(2\pi t+\varphi )e^{-\kappa (T-t)}}{4{\pi }^2+{\kappa }^2}}+\kappa {\displaystyle \frac{sin(2\pi T+\varphi )-sin(2\pi t+\varphi )e^{-\kappa (T-t)}}{4{\pi }^2+{\kappa }^2}},\hfill \end{array}$

(16) $\begin{array}{cc}\hfill f_3(t;T)& =-\left(e^{\alpha +{\displaystyle \frac{{\beta }^2}{2}}}-{\displaystyle \frac{\gamma }{\kappa }}\right)(T-t)-{\displaystyle \frac{\gamma }{\kappa }}{\displaystyle \frac{1-e^{-\kappa (T-t)}}{\kappa }}\hfill \\ \hfill & +{\int }_0^{T-t}{e}^{\alpha +{\displaystyle \frac{{\beta }^2}{2}}-{\displaystyle \frac{\gamma +\nu \beta \delta }{\kappa }}(1-e^{-\kappa v})+{\displaystyle \frac{{\delta }^2}{2{\kappa }^2}}{(1-e^{-\kappa v})}^2}dv.\hfill \end{array}$

Proof. 

See Appendix A.2. □

The futures price $F(t,s,y;T)$ reflects the value of the underlying commodity at the future expiration date T. Its derivative with respect to the current commodity price s is always positive, indicating direct proportionality. Conversely, the futures price is inversely related to the convenience yield $Y_t$, as $f_1(t;T)<0$. This inverse relationship aligns with the theory of storage, where higher convenience yields reduce the reliance on futures contracts due to the benefits of physical ownership, leading to lower futures prices.

3.2.2. Extraction of Commodity Price and Convenience Yield

On a given trading day t, the futures price depends on the unobservable variables $S_t$ and $Y_t$. Proposition 2 provides explicit formulas for extracting these state variables using observed futures prices from two contracts with different maturities.

Proposition 2. 

Given two futures contracts with maturities $T_1$ and $T_2$, where $0<T_1<T_2$, and their corresponding prices $F_t^{T_1}$ and $F_t^{T_2}$ observed on day t, the log price $X_t$ and convenience yield $Y_t$ can be expressed as

(17) $\begin{array}{cc}\hfill X_t& =ln{F}_t^{T_2}+{\displaystyle \frac{ln\left(\frac{F_t^{T_1}}{F_t^{T_2}}\right)+A(t;T_2)-A(t;T_1)}{e^{-\kappa (T_2-t)}+e^{-\kappa (T_1-t)}}}\left(1-e^{-\kappa (T_2-t)}\right)-A(t;T_2),\hfill \end{array}$

(18) $\begin{array}{cc}\hfill Y_t& =\kappa {\displaystyle \frac{ln\left(\frac{F_t^{T_1}}{F_t^{T_2}}\right)+A(t;T_2)-A(t;T_1)}{e^{-\kappa (T_2-t)}+e^{-\kappa (T_1-t)}}},\hfill \end{array}$

where

(19) $\begin{array}{cc}\hfill A(t;T_i)& =f_0(t;T_i)+f_2(t;T_i){\overline{y}}_1+f_3(t;T_i)\lambda ,\hfill \end{array}$

for $i=1,2$.

Proof. 

See Appendix A.3. □

Proposition 2 enables the practical extraction of $S_t$ and $Y_t$ from observed futures prices, bridging theoretical pricing models with empirical applications. This facilitates parameter estimation and a deeper analysis of market dynamics.

4. Numerical Results and Discussion

This section provides an in-depth analysis of the numerical performance and practical applications of the analytical pricing framework presented in Theorem 1. The analysis is divided into three main components. First, the accuracy and computational efficiency of the closed-form formula from Proposition 1 are evaluated by benchmarking its results against Monte Carlo (MC) simulations, using the Euler–Maruyama (EM) approximation for jump-diffusion processes, as detailed in Chapter 3 of Glasserman [20]. This validation demonstrates the reliability of the analytical derivation and highlights its advantages in terms of computational speed and precision.

Second, a sensitivity analysis explores the effects of jumps and seasonal components on futures prices. This analysis offers valuable insights into the robustness of the model and its capacity to incorporate these essential features. Finally, an empirical study assesses the applicability of the proposed model (12) using real-world data, specifically futures prices for ribbed smoked rubber sheets traded on the Thailand Futures Exchange (TFEX). The empirical results underscore the model’s effectiveness in capturing market dynamics and estimating key parameters through nonlinear regression.

All numerical experiments were conducted using MATHEMATICA 13.0 on a notebook equipped with an Apple M1 processor (8 GB RAM) running macOS Sonoma 14.0.

4.1. Accuracy and Efficiency of the Closed-Form Formula for Futures Pricing

Example 1. 

This example evaluates the accuracy and computational efficiency of the closed-form Formula (12) for futures pricing under the MMRJD model (1), comparing it against MC simulations. The MC method employs the Euler–Maruyama (EM) approximation, with a time step of $\mathsf{\Delta }t=1/252$, representing daily intervals over a trading year. The initial time was set to $t_0=0$ with the initial states of the spot price and convenience yield given as $(s_0,y_0)=(50,-0.50)$.

The parameter vector $\overline{\mathbf{p}}$ was initialized as follows: for the diffusion component, the risk-free interest rate was $r=0.05$, and the volatilities of the log price and convenience yield were $\sigma =0.25$ and $\xi =0.30$, respectively. The mean-reversion rate and the market price of convenience yield risk were set to $\kappa =0.50$ and $\eta =0.20$. For the seasonal component, the stationary level of long-term equilibrium was ${\overline{y}}_0=0.50$, while the amplitude and phase shift of the sinusoidal function were specified as ${\overline{y}}_1=0.30$ and $\varphi =1.50$, respectively. For the Poisson jump process, the intensity rate was $\lambda =0.50$, and the means of the jump sizes for the log price and convenience yield were $\alpha =0.20$ and $\gamma =-0.15$ with standard deviations $\beta =0.15$ and $\delta =0.25$, respectively. The correlation between the two Wiener processes was set to $\rho =0.15$, while the correlation between the two jump components was set to $\nu =-0.25$.

To emphasize the distinct roles of jump and seasonal components in shaping futures price dynamics, parameters of interest $(\lambda ,\alpha ,\beta ,\gamma ,\delta ,\nu ,{\overline{y}}_1,\varphi )$ were varied within specified ranges and fixed $\tau =0.100,0.101,0.102$, as shown in Figure 1. Each subfigure highlights the effect of a single parameter, while others were fixed at their baseline values $\overline{\mathbf{p}}$.

Figure 1a depicts the effect of the jump intensity rate λ, varied within the range $0.25\le \lambda \le 0.75$, while Figure 1b illustrates the impact of the mean of log price jumps α, which varied within $-0.5\le \alpha \le 0.5$. Figure 1c,d examine the influence of the standard deviation of log price jumps β and the mean of convenience yield jumps γ, which are varied within $0.25\le \beta \le 0.75$ and $-0.5\le \gamma \le 0.5$, respectively. Figure 1e demonstrates the effect of the standard deviation of convenience yield jumps δ, varied within $0.25\le \delta \le 0.75$, while Figure 1f shows the correlation between jumps ν, varied within $-0.5\le \nu \le 0.5$. Finally, Figure 1g,h present the effects of the seasonal component, specifically the amplitude ${\overline{y}}_1$ and the phase shift ϕ, varied within $0.25\le {\overline{y}}_1\le 0.75$ and $0\le \varphi \le 2\pi$, respectively.

The futures price is an increasing function of λ, α, β, and δ, indicating that higher jump intensity, mean and the variance in jumps in spot price, as well as variance in jumps in convenience yield contribute to higher futures prices. Conversely, the futures price decreases with γ, ν, and ${\overline{y}}_1$, suggesting that larger convenience yield jumps, stronger jump correlation, and seasonal amplitude reduce futures prices. Additionally, the futures price exhibits a periodic pattern with respect to ϕ, highlighting the influence of seasonality in the model.

To demonstrate the convergence of MC simulations to the closed-form Formula (12), we fixed the parameter of interest as the highest value in Figure 1 with time to maturity $\tau =1$ and varied the number of sample paths $n_p$ from $2\times {10}^5$ to $2\times {10}^7$. In Figure 2, panel (a) considers the jump intensity rate $\lambda =0.75$, while panel (b) examines the mean of log price jumps $\alpha =0.5$. Panel (c) presents results for the standard deviation of log price jumps $\beta =0.75$, and panel (d) evaluates the mean of convenience yield jumps $\gamma =0.5$. Similarly, panels (e) and (f) illustrate the convergence behavior for the standard deviation of convenience yield jumps $\delta =0.75$ and the correlation between jumps $\nu =0.5$. Finally, panels (g) and (h) depict the convergence results for the amplitude of the seasonal component ${\overline{y}}_1=0.75$ and the phase shift $\varphi =\pi$, respectively. In each case, as $n_p$ increases, the MC simulation results converge to the futures prices computed using the closed-form formula (blue lines), confirming the accuracy and reliability of the formula.

Next, we examine the computational efficiency of the closed-form Formula (12) in comparison to MC simulations. Table 1 summarizes the effect of the number of sample paths, $n_p$, on computational performance using the baseline parameter set $\overline{\mathbf{p}}$. For $T=1$, the sample paths were varied as $n_p=1\times {10}^3$, $5\times {10}^3$, $1\times {10}^4$, $5\times {10}^4$, and $1\times {10}^5$. The computational time for MC simulations increased proportionally with the number of sample paths, peaking at $2.8366$ s for $n_p=1\times {10}^5$. In contrast, the closed-form formula consistently required only $0.0035$ s regardless of $n_p$. This results in a substantial computational advantage with the closed-form formula being up to 811-fold faster than the MC simulation for $n_p=1\times {10}^5$.

In conclusion, this analysis validates the closed-form formula’s accuracy and computational efficiency for futures pricing under the MMRJD model (1), showcasing its robust performance across a wide range of parameter settings.

4.2. Impact of Jumps and Seasonality on Futures Pricing

Example 2. 

This example explores the sensitivity of futures prices to changes in model parameters through a systematic sensitivity analysis. Sensitivity analysis quantifies how variations in input parameters influence model outputs, offering insights into the robustness of the model under different assumptions [21]. The analysis uses the baseline parameter set $\overline{\mathbf{p}}$ with the initial state of the system defined as $(s_0,y_0)=(50,-0.50)$. Each parameter in the model is individually increased by $10\%$ from its baseline value, and the resulting percentage absolute change in the futures price is computed using the closed-form Formula (12) [22].

Figure 3 presents a bar chart illustrating the percentage absolute changes in futures prices resulting from a 10% increase in each parameter. Parameters represented by taller bars exhibit greater sensitivity, reflecting a stronger influence on the futures price. Among the parameters, the phase shift ϕ shows the highest sensitivity within the seasonal component with a percentage change of approximately $0.2800\%$. Similarly, the intensity rate λ demonstrates the highest sensitivity within the jump component with a percentage change of approximately $0.1687\%$. These findings underscore the critical roles of ϕ and λ in shaping the dynamics of futures prices under the proposed model.

4.3. Empirical Study on Ribbed Smoked Rubber Sheet Traded in Thailand

Example 3. 

This example presents an empirical study to evaluate the applicability and effectiveness of the MMRJD model (1) in capturing the dynamics of commodity futures prices. The primary objective is to assess the model’s performance in explaining observed market prices and identifying key parameters influencing the term structure of futures prices.

The study utilizes a data set comprising 61 futures contracts for Ribbed Smoked Rubber Sheet (RSS) traded on the TFEX (https://www.setsmart.com, accessed on 28 March 2024) between 3 January 2019 and 27 December 2023. The futures prices are categorized into two groups based on their maturity: the first closest maturity ($F1$) and the second closest maturity ($F2$). These are depicted in Figure 4 with $F1$ represented by the red line and $F2$ by the blue line. This classification provides a clear representation of the term structure of RSS futures prices, offering insights into the temporal dynamics of the market.

In this study, we calibrate the MMRJD model (1) using futures contracts observed over K = 1138 days, which are denoted as $0=t_0<t_1<t_2<\cdots <t_K=5$. For each day $t_i$, there are two different maturities of futures contracts available [23]. The observed futures prices can be decomposed into expected and unexpected components, representing the no-arbitrage futures price and random speculation, respectively, as follows:

(20)$ln{F}_{i,j}=lnF(t_i,s_i,y_i;T_{i,j},\mathbf{p})+{\mathrm{e}}_{i,j},$

for $i=1,2,\dots ,K$ and $j=1,2$. Here, $lnF(t_i,s_i,y_i;T_{i,j},\mathbf{p})$ represents the logarithm of the no-arbitrage futures price based on the MMRJD model (1), and ${\mathrm{e}}_{i,j}$ denotes the residual noise, which captures random speculation in the futures market. The residual ${\mathrm{e}}_{i,j}$ is assumed to follow a normal distribution with zero mean and variance ${ϵ}^2$, as detailed in Section 3 of Gonzato and Sgarra [17].

For the policy interest rate r, the risk-free interest rate is assumed to be constant with an average value of $0.0105$ over the study period (https://www.bot.or.th/th/our-roles/monetary-policy/mpc-publication/policy-interest-rate.html, accessed on 23 May 2024). The parameter vector $\mathbf{p}=(\sigma ,\kappa ,{\overline{y}}_0,{\overline{y}}_1,\varphi ,\eta ,\xi ,\rho ,\lambda ,\alpha ,\beta ,\gamma ,\delta ,\nu )$ is estimated using a nonlinear regression approach, which minimizes the discrepancy between the observed and modeled logarithm of futures prices. The objective function for the parameter estimation is given by

(21) ${\mathbf{p}}^{*}=arg\underset{\mathbf{p}}{min}\sum _{i=1}^K\sum _{j=1}^2{\left(ln{F}_{i,j}-lnF(t_i,s_i,y_i;T_{i,j},\mathbf{p})\right)}^2,$

subject to the following parameter constraints:

$\begin{array}{c}\sigma ,\kappa ,{\overline{y}}_1,\xi ,\lambda ,\beta ,\delta \in {\mathbb{R}}_{+},\\ {\overline{y}}_0,\eta ,\alpha ,\gamma \in \mathbb{R},\\ \rho ,\nu \in (-1,1).\end{array}$

The optimization was performed using commercial software equipped with numerical methods, such as the “NMinimize” command in Wolfram MATHEMATICA version 13. Table 2 provides a detailed summary of the estimated parameters, including their standard errors and $95\%$ confidence intervals, as described in Chapter 7 of Ross [24].

Figure 5 illustrates the distinct components of the closed-form formula for futures pricing under the MMRJD model (1), which was calibrated using the parameter vector ${\mathbf{p}}^{*}$ obtained from empirical analysis. The figure consists of six subfigures, each highlighting a unique aspect of the model’s structure.

Figure 5a depicts the diffusion component $f_0(t;T)$, which accounts for continuous price fluctuations and evolves over time to maturity τ. Figure 5b presents the coefficient of convenience yield $f_1(t;T)$, demonstrating a diminishing influence of the convenience yield as the maturity increases. Figure 5c illustrates the seasonal component $f_2(t;T)$, which captures periodic price effects driven by the sinusoidal function. Figure 5d highlights the jump component $f_3(t;T)$, showing the increasing impact of discrete jumps as the time to maturity extends.

Figure 5e visualizes the 3D probability density function $g(\acute{x},\acute{y})$, representing the joint distribution of log price jumps and convenience yield jumps, while Figure 5f provides a contour plot of the same distribution. Together, these subfigures offer a detailed decomposition of the key drivers of futures prices, showcasing the interactions between continuous fluctuations, seasonal trends, and discrete market shocks. This comprehensive visualization provides a deeper understanding of price dynamics under the MMRJD model (1).

Utilizing the parameter estimates ${\mathbf{p}}^{*}$ from Table 2, we apply the extraction formulas (17) and (18) provided in Proposition 2 to compute the extracted values of $S_t$ and $Y_t$ from observed futures prices over the period 3 January 2019 to 30 December 2020. This timeframe notably includes the period of significant price jumps in RSS futures, which are attributed to the market disruptions caused by the COVID-19 pandemic. The extracted results are plotted in Figure 6.

Figure 6a displays the extracted values of $S_t$ against the observed spot prices of RSS. The plot reveals a strong positive correlation between the two sequences, confirming the accuracy of the extraction process in capturing the dynamics of spot prices from observed futures prices. This correlation underscores the effectiveness of the MMRJD model (1) in reflecting real-world market behaviors, particularly in volatile periods.

Figure 6b illustrates the extracted values of $Y_t$ (convenience yields) against the modified long-run equilibrium $\vartheta \left(t\right)$ derived in (5). The extracted convenience yields exhibit fluctuations around the long-run equilibrium, which is consistent with the mean-reverting behavior expected in commodity markets. However, during periods of significant price jumps, the extracted convenience yields deviate markedly from the long-run equilibrium and are predominantly negative. This pattern aligns with the negative value of ${\nu }^{*}$ (the correlation coefficient between $J_t^{\left(X\right)}$ and $J_t^{\left(Y\right)}$) shown in Table 2, indicating an inverse relationship between price jumps and convenience yield jumps. Negative convenience yields imply higher inventory costs and are commonly observed during periods of supply disruptions or heightened market stress.

These results highlight the ability of the MMRJD model (1) to capture the nuanced interplay between spot prices, convenience yields, and market shocks. The extracted values of $S_t$ and $Y_t$ not only validate the model’s predictive accuracy but also provide deeper insights into market behavior during volatile periods, making it a valuable tool for both empirical analysis and practical applications in risk management and derivative pricing.

Finally, we adopt the methodology outlined in [25] to compare the MMRJD model (1) with Schwartz’s two-factor model [5] that does not incorporate jumps and seasonal effects. To evaluate the predictive performance of the MMRJD model (1), we define the percentage absolute errors (${\mathrm{E}}_{i,j}$) as

${\mathrm{E}}_{i,j}={\displaystyle \frac{|{\mathrm{e}}_{i,j}|}{lnF(t_i,s_i,y_i;T_{i,j},\mathbf{p})}}\times 100\%,$

where $lnF(t_i,s_i,y_i;T_{i,j},\mathbf{p})$ represents the logarithm of the futures price derived from (12). For comparison, Schwartz’s two-factor model is considered, which excludes seasonal and jump components by setting ${\overline{y}}_1=\lambda =0$. The logarithm of the futures price in Schwartz’s model is denoted by $ln\tilde{F}(t_i,s_i,y_i;T_{i,j},\tilde{\mathbf{p}})$, where the parameter vector is reduced to $\tilde{\mathbf{p}}=(\sigma ,\kappa ,{\overline{y}}_0,\eta ,\xi ,\rho )$.

The percentage absolute errors for Schwartz’s two-factor model (${\tilde{\mathrm{E}}}_{i,j}$) are defined as

${\tilde{\mathrm{E}}}_{i,j}={\displaystyle \frac{|{\tilde{\mathrm{e}}}_{i,j}|}{ln\tilde{F}(t_i,{\tilde{s}}_i,{\tilde{y}}_i;T_{i,j},\tilde{\mathbf{p}})}}\times 100\%.$

It is important to note that based on the above definitions of percentage absolute errors, we assume the log-futures prices used in this study are nonzero to ensure that the computations are well defined.

As shown in Figure 7, the errors ${\mathrm{E}}_{i,j}$ associated with the MMRJD model (1) demonstrate both a smaller mean and variance compared to ${\tilde{\mathrm{E}}}_{i,j}$ from Schwartz’s two-factor model. This highlights the superior performance of the MMRJD model (1) in capturing the intricate patterns present in the RSS futures price data, significantly outperforming Schwartz’s simpler two-factor model.

5. Conclusions

This study has introduced the MMRJD model (1), effectively capturing key characteristics of commodity prices, including correlated jumps, seasonal effects, and mean-reversion dynamics in both prices and convenience yields. By extending Schwartz’s two-factor model, the MMRJD framework (1) has provided a comprehensive approach for modeling market shocks, periodic fluctuations, and stochastic convenience yields. The derived closed-form formula for the CMGF, leveraging the Feynman–Kac theorem, has demonstrated significant computational efficiency and high accuracy compared to MC simulations. Empirical application to historical RSS futures prices traded on TFEX has validated the model’s robustness, highlighting the critical roles of jumps and seasonality in shaping market behavior. The findings show that convenience yields are negative, reflecting higher inventory costs, and increase with rising spot prices. Moreover, the MMRJD model (1) has outperformed the Schwartz two-factor model in capturing the observed dynamics of spot prices and convenience yields, offering actionable insights for traders, risk managers, and policymakers in commodity markets.

Author Contributions

A.S., S.R. and B.D. contributed equally to the analysis and writing of this paper. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are contained within the article.

Acknowledgments

This research was supported by the National Research Council of Thailand (NRCT) under grant number NRCT5-RGJ63016-150 and by Walailak University Partial Funding (Contract No. 01/2562) for the first and second authors. We extend our deepest gratitude to Willi Jä ger and Peter Bastian from the Interdisciplinary Center for Scientific Computing (IWR) at Heidelberg University for their invaluable discussions, which provided crucial insights for solving the PIDEs in this paper. We are also grateful for the insightful feedback from the anonymous reviewers, which has greatly enhanced the quality of this paper.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:

MMRJD	multivariate mean-reverting jump-diffusion
CMGF	conditional moment generating function
PIDE	partial integro-differential equation
ODE	ordinary differential equation
MC	Monte Carlo
EM	Euler–Maruyama
TFEX	Thailand Futures Exchange

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Comparative study on the variation in futures prices derived from the closed-form Formula (12) across diverse parameter ranges [Forumla omitted. See PDF.] and time to maturity [Forumla omitted. See PDF.]. The plots highlight the effects of seasonal and jump components, as detailed in Example 1.

Figure 2. Convergence analysis of MC simulations compared to the closed-form Formula (12), illustrating accuracy across varying sample sizes and selected parameter values [Forumla omitted. See PDF.], as detailed in Example 1.

Figure 3. Impact of jumps and seasonality on futures prices under the MMRJD model (1). The figure illustrates the percentage change in futures prices resulting from a 10% increase in each parameter of the baseline set [Forumla omitted. See PDF.], with initial conditions [Forumla omitted. See PDF.], as detailed in Example 2.

Figure 4. Futures prices for the first and second closest maturities of RSS traded on the TFEX from 3 January 2019 to 27 December 2023. These data are used for parameter estimation, as discussed in Example 3.

Figure 5. The characteristic of specific components in our closed-form Formula (12) corresponding to the estimated parameters [Forumla omitted. See PDF.], as detailed in Example 3.

Figure 6. Extraction of spot prices and convenience yields for RSS futures using the Formulas (17) and (18) and the parameter estimates [Forumla omitted. See PDF.] from Table 2, as detailed in Example 3.

Figure 7. Distributions of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], illustrating the superior performance of the MMRJD model (1) in capturing the nuanced patterns in RSS futures price data compared to Schwartz’s two-factor model. This comparison highlights the MMRJD model’s enhanced predictive accuracy, as detailed in Example 3.

Table 1

Comparison of the computational efficiency of the closed-form Formula (12) and MC simulations. The second column reports the absolute percentage error ($\varsigma =100\%\times |\mathrm{model}-\mathrm{simulation}|/\mathrm{model}$) between the closed-form formula and MC results. The computation times for MC simulations ($t_{MC}$) and the closed-form formula ($t_{CF}$) are shown in the third and fourth columns, respectively. The final column provides the approximate fold reduction in computational time achieved by the closed-form formula, as illustrated in Example 1.

${\mathit{n}}_{\mathit{p}}$	$\mathit{\varsigma }$	${\mathit{t}}_{\mathbf{MC}}$	${\mathit{t}}_{\mathbf{CF}}$	Reduction
(Samples)	(%)	(Seconds)	(Seconds)	(Folds)
$1\times {10}^3$	$0.6695$	$0.0870$	$0.0035$	25
$5\times {10}^3$	$0.5314$	$0.2715$	$0.0035$	78
$1\times {10}^4$	$0.4228$	$0.6013$	$0.0035$	172
$5\times {10}^4$	$0.1139$	$1.8431$	$0.0035$	527
$1\times {10}^5$	$0.0214$	$2.8366$	$0.0035$	811

Table 2

Parameter estimates of the MMRJD model (1), as detailed in Example 3, including standard errors and $95\%$ confidence intervals. These estimates are obtained via least squares using two sequences of daily futures prices ($F1$ and $F2$), each consisting of K = 1138 observation dates.

Model	Parameter	Estimated	Standard Error	Confidence Interval
MMRJD	${\sigma }^{*}$	$1.5515$	$0.0054$	$(1.5407,1.5622)$
	${\kappa }^{*}$	$7.6128$	$0.8064$	$(6.0322,9.1935)$
	${\overline{y}}_0^{*}$	$0.3687$	$0.1720$	$(0.0314,0.7059)$
	${\overline{y}}_1^{*}$	$0.4882$	$0.1945$	$(0.1070,0.8694)$
	${\varphi }^{*}$	$3.5941$	$0.0174$	$(3.5599,3.6283)$
	${\eta }^{*}$	$0.1339$	$0.0531$	$(0.0297,0.2381)$
	${\xi }^{*}$	$0.0125$	$0.0329$	$(0.0052,0.0771)$
	${\rho }^{*}$	$0.4556$	$0.0109$	$(0.4341,0.4770)$
	${\lambda }^{*}$	$1.1178$	$0.1147$	$(0.8929,1.3426)$
	${\alpha }^{*}$	$-0.4927$	$0.0422$	$(-0.5755,-0.4098)$
	${\beta }^{*}$	$0.6393$	$0.0630$	$(0.5156,0.7628)$
	${\gamma }^{*}$	$-0.6169$	$0.0414$	$(-0.6980,-0.5357)$
	${\delta }^{*}$	$1.9610$	$0.1719$	$(1.6240,2.2980)$
	${\nu }^{*}$	$-0.5528$	$0.0345$	$(-0.6205,-0.4850)$
Schwartz’s	${\sigma }^{*}$	$0.5924$	$0.0525$	$(0.4893,0.6954)$
Two-Factor	${\kappa }^{*}$	$0.0245$	$0.0297$	$(0.0037,0.0828)$
	${\overline{y}}_0^{*}$	$1.2993$	$0.1887$	$(0.9294,1.6692)$
	${\eta }^{*}$	$-0.4804$	$0.1402$	$(-0.7552,-0.2055)$
	${\xi }^{*}$	$2.4081$	$0.5521$	$(1.3259,3.4902)$
	${\rho }^{*}$	$0.9125$	$0.0324$	$(0.8488,0.9761)$

Appendix A

Appendix A.1. Proof of Theorem 1

Proof. 

For fixed $u\in \mathbb{R}$, the CMGF of $X_T$ given $X_t=x$ and $Y_t=y$ is defined as(A1)$$M(t,x,y;u,T)=M_{X_T}^{\left(u\right)}(t,x,y),$$ for $(t,x,y)\in [0,T]\times {\mathbb{R}}^2$. By applying the Feynman–Kac theorem [17], this function satisfies the following partial integro-differential equation (PIDE):(A2)$$\begin{array}{c}{\displaystyle \frac{\partial M}{\partial t}}+\left(\mu -y\right){\displaystyle \frac{\partial M}{\partial x}}+{\displaystyle \frac{{\sigma }^2}{2}}{\displaystyle \frac{{\partial }^{2}M}{\partial x^2}}+\kappa \left(\vartheta \left(t\right)-y\right){\displaystyle \frac{\partial M}{\partial y}}+{\displaystyle \frac{{\xi }^2}{2}}{\displaystyle \frac{{\partial }^{2}M}{\partial y^2}}+\rho \sigma \xi {\displaystyle \frac{{\partial }^{2}M}{\partial x\partial y}}\hfill \\ \hfill +\lambda {\int }_{\mathbb{R}}{\int }_{\mathbb{R}}\left(M(t,x+\acute{x},y+\acute{y};u,T)-M(t,x,y;u,T)\right)g(\acute{x},\acute{y})d\acute{x}d\acute{y}=0,\end{array}$$ for $(t,x,y)\in [0,T)\times {\mathbb{R}}^2$ with the terminal condition(A3)$$M(T,x,y;u,T)=e^{ux}.$$

To solve (A2) analytically, we perform a change of variables by setting $\tau =T-t$, noting that $\partial \tau =-\partial t$. The PIDE becomes(A4)$$\begin{array}{c}{\displaystyle \frac{\partial M}{\partial \tau }}=\left(\mu -y\right){\displaystyle \frac{\partial M}{\partial x}}+{\displaystyle \frac{{\sigma }^2}{2}}{\displaystyle \frac{{\partial }^{2}M}{\partial x^2}}+\kappa \left(\vartheta (T-\tau )-y\right){\displaystyle \frac{\partial M}{\partial y}}+{\displaystyle \frac{{\xi }^2}{2}}{\displaystyle \frac{{\partial }^{2}M}{\partial y^2}}+\rho \sigma \xi {\displaystyle \frac{{\partial }^{2}M}{\partial x\partial y}}\hfill \\ \hfill +\lambda {\int }_{\mathbb{R}}{\int }_{\mathbb{R}}\left(M(\tau ,x+\acute{x},y+\acute{y};u,T)-M(\tau ,x,y;u,T)\right)g(\acute{x},\acute{y})d\acute{x}d\acute{y},\end{array}$$ for $(\tau ,x,y)\in (0,T]\times {\mathbb{R}}^2$ with the initial condition(A5)$$M(0,x,y;u,T)=e^{ux}.$$

We propose a solution of the following form:(A6)$$M(\tau ,x,y;u,T)=e^{m_1(\tau ;u,T)x+m_2(\tau ;u,T)y+m_3(\tau ;u,T)},$$ where $m_1(\tau ;u,T)$, $m_2(\tau ;u,T)$, and $m_3(\tau ;u,T)$ are functions of $\tau$. Substituting (A6) into (A4), the partial derivatives and integral term are evaluated as in the original proof, leading to the following system of ODEs:(A7)$$\begin{array}{cc}\hfill {\displaystyle \frac{d{m}_1}{d\tau }}& =0,\hfill \end{array}$$(A8)$$\begin{array}{cc}\hfill {\displaystyle \frac{d{m}_2}{d\tau }}& =-\kappa m_2(\tau ;u,T)-m_1(\tau ;u,T),\hfill \end{array}$$(A9)$$\begin{array}{cc}\hfill {\displaystyle \frac{d{m}_3}{d\tau }}& =\mu m_1(\tau ;u,T)+{\displaystyle \frac{{\sigma }^2}{2}}{m}_1^2(\tau ;u,T)+\kappa \vartheta (T-\tau )m_2(\tau ;u,T)\hfill \\ \hfill & +{\displaystyle \frac{{\xi }^2}{2}}{m}_2^2(\tau ;u,T)+\rho \sigma \xi m_1(\tau ;u,T)m_2(\tau ;u,T)\hfill \\ \hfill & +\lambda \left(e^{\alpha m_1(\tau ;u,T)+{\displaystyle \frac{{\beta }^2}{2}}{m}_1^2(\tau ;u,T)+\gamma m_2(\tau ;u,T)+{\displaystyle \frac{{\delta }^2}{2}}{m}_2^2(\tau ;u,T)+\nu \beta \delta m_1(\tau ;u,T)m_2(\tau ;u,T)}-1\right).\hfill \end{array}$$

The initial conditions derived from (A5) are(A10)$$\begin{array}{c}\hfill m_1(0;u,T)=u,m_2(0;u,T)=0,m_3(0;u,T)=0.\end{array}$$

The solutions of the ODEs are(A11)$$\begin{array}{cc}\hfill m_1(\tau ;u,T)& =u,\hfill \end{array}$$(A12)$$\begin{array}{cc}\hfill m_2(\tau ;u,T)& =-u{\displaystyle \frac{1-e^{-\kappa \tau }}{\kappa }},\hfill \end{array}$$(A13)$$\begin{array}{cc}\hfill m_3(\tau ;u,T)& =u{\psi }_1(\tau ;T)+\lambda {\psi }_2(\tau ;u,T)+u{\overline{y}}_1{\psi }_3(\tau ;T)+u^2{\psi }_4(\tau ;T),\hfill \end{array}$$ where$$\begin{array}{cc}\hfill {\psi }_1(\tau ;T)& =\left(r-{\displaystyle \frac{{\sigma }^2}{2}}-{\overline{y}}_0+{\displaystyle \frac{\eta }{\kappa }}\right)\tau +\left({\overline{y}}_0-{\displaystyle \frac{\eta }{\kappa }}\right){\displaystyle \frac{1-e^{-\kappa \tau }}{\kappa }},\hfill \\ \hfill {\psi }_2(\tau ;u,T)& =-u\left(e^{\alpha +{\displaystyle \frac{{\beta }^2}{2}}}-1-{\displaystyle \frac{\gamma }{\kappa }}\right)\tau -{\displaystyle \frac{u\gamma }{\kappa }}{\displaystyle \frac{1-e^{-\kappa \tau }}{\kappa }}\hfill \\ \hfill & +{\int }_0^{\tau }{e}^{u\alpha -{\displaystyle \frac{u\gamma }{\kappa }}(1-e^{-\kappa s})+{\displaystyle \frac{u^2{\beta }^2}{2}}-{\displaystyle \frac{u^2\nu \beta \delta }{\kappa }}(1-e^{-\kappa s})+{\displaystyle \frac{u^2{\delta }^2}{2{\kappa }^2}}{(1-e^{-\kappa s})}^2}ds-\tau ,\hfill \\ \hfill {\psi }_3(\tau ;T)& ={\displaystyle \frac{cos(2\pi T+\varphi )-cos\left(2\pi \right(T-\tau )+\varphi )}{2\pi }}\hfill \\ \hfill & -2\pi {\displaystyle \frac{cos(2\pi T+\varphi )-cos(2\pi (T-\tau )+\varphi )e^{-\kappa \tau }}{4{\pi }^2+{\kappa }^2}}\hfill \\ \hfill & +\kappa {\displaystyle \frac{sin(2\pi T+\varphi )-sin(2\pi (T-\tau )+\varphi )e^{-\kappa \tau }}{4{\pi }^2+{\kappa }^2}},\hfill \\ \hfill {\psi }_4(\tau ;T)& =\left({\displaystyle \frac{{\sigma }^2}{2}}-{\displaystyle \frac{\rho \sigma \xi }{\kappa }}+{\displaystyle \frac{{\xi }^2}{2{\kappa }^2}}\right)\tau +\left({\displaystyle \frac{\rho \sigma \xi }{\kappa }}-{\displaystyle \frac{{\xi }^2}{2{\kappa }^2}}\right){\displaystyle \frac{1-e^{-\kappa \tau }}{\kappa }}-{\displaystyle \frac{{\xi }^2}{4\kappa }}{\left({\displaystyle \frac{1-e^{-\kappa \tau }}{\kappa }}\right)}^2.\hfill \end{array}$$

Substituting these solutions into (A6) completes the proof. □

Appendix A.2. Proof of Proposition 1

Proof. 

The futures price of the commodity under the MMRJD model (1) is obtained by evaluating the CMGF (7) at $u=1$. Using the definition of futures prices (11), the futures price is expressed as(A14)$$\begin{array}{c}\hfill F(t,s,y;T)=M(t,lns,y;1,T)=e^{lns-\left({\displaystyle \frac{1-e^{-\kappa (T-t)}}{\kappa }}\right)y+\mathrm{\Psi }(t;1,T)+\mathrm{{\rm Y}}(t;1,T){\overline{y}}_1+\mathrm{\Phi }(t;1,T)\lambda },\end{array}$$ for $(t,s,y)\in [0,T]\times {\mathbb{R}}_{+}\times \mathbb{R}$, where $\mathrm{\Psi }(t;1,T)$, $\mathrm{{\rm Y}}(t;1,T)$, and $\mathrm{\Phi }(t;1,T)$ correspond to the diffusion, seasonal, and jump components of the CMGF, respectively.

Simplifying (A14), the futures price can be rewritten as(A15)$$F(t,s,y;T)=s{e}^{f_0(t;T)+f_1(t;T)y+f_2(t;T){\overline{y}}_1+f_3(t;T)\lambda },$$ where the functions $f_0(t;T)$, $f_1(t;T)$, $f_2(t;T)$, and $f_3(t;T)$ are derived directly from the exponential terms in (A14) as follows:$$\begin{array}{cc}\hfill f_0(t;T)& =\mathrm{\Psi }(t;1,T),\hfill \\ \hfill f_1(t;T)& =-{\displaystyle \frac{1-e^{-\kappa (T-t)}}{\kappa }},\hfill \\ \hfill f_2(t;T)& =\mathrm{{\rm Y}}(t;1,T),\hfill \\ \hfill f_3(t;T)& =\mathrm{\Phi }(t;1,T).\hfill \end{array}$$

Thus, we have shown that the futures price can be expressed in the form of Equation (12), completing the proof. □

Appendix A.3. Proof of Proposition 2

Proof. 

Let the futures prices for the closest and second-closest maturities on a fixed trading day t be $F_t^{T_1}$ and $F_t^{T_2}$, where $T_1$ and $T_2$ denote the respective maturities. From Proposition 1, the logarithm of the futures price is expressed as(A16)$$\begin{array}{cc}\hfill ln{F}_t^{T_1}& =ln{S}_t+A(t;T_1)+f_1(t;T_1)Y_t,\hfill \end{array}$$(A17)$$\begin{array}{cc}\hfill ln{F}_t^{T_2}& =ln{S}_t+A(t;T_2)+f_1(t;T_2)Y_t,\hfill \end{array}$$ where $A(t;T_i)$ is defined as in (19) for $i=1,2$, and $f_1(t;T_i)$ is given in (14). Subtracting (A17) from (A16) yields(A18)$$ln{\displaystyle \frac{F_t^{T_1}}{F_t^{T_2}}}=A(t;T_1)-A(t;T_2)+\left(f_1(t;T_1)-f_1(t;T_2)\right)Y_t.$$

Substituting $f_1(t;T_i)=-{\displaystyle \frac{1-e^{-\kappa (T_i-t)}}{\kappa }}$ from (14) into (A18), we obtain$$ln{\displaystyle \frac{F_t^{T_1}}{F_t^{T_2}}}=A(t;T_1)-A(t;T_2)+{\displaystyle \frac{e^{-\kappa (T_2-t)}-e^{-\kappa (T_1-t)}}{\kappa }}{Y}_t.$$ Rearranging for $Y_t$ gives(A19)$$Y_t=\kappa {\displaystyle \frac{ln\frac{F_t^{T_1}}{F_t^{T_2}}+A(t;T_2)-A(t;T_1)}{e^{-\kappa (T_2-t)}-e^{-\kappa (T_1-t)}}}.$$

Finally, substituting (14) and (18) into (A17) allows us to express $S_t$ as(A20)$$S_t=F_t^{T_2}{e}^{\left({\displaystyle \frac{ln\frac{F_t^{T_1}}{F_t^{T_2}}+A(t;T_2)-A(t;T_1)}{e^{-\kappa (T_2-t)}-e^{-\kappa (T_1-t)}}}\right)\left(1-e^{-\kappa (T_2-t)}\right)-A(t;T_2)}.$$

Thus, both $S_t$ and $Y_t$ are expressed in terms of the two closest futures prices, $F_t^{T_1}$ and $F_t^{T_2}$, completing the proof. □

References

1. Kaldor, N. Speculation and Economic Stability. Rev. Econ. Stud.; 1939; 7, pp. 1-27. [DOI: https://dx.doi.org/10.2307/2967593]

2. Working, H. The Theory of the Price of Storage. Am. Econ. Rev.; 1949; 29, pp. 1254-1262.

3. Brennan, M.J. The Supply of Storage. Am. Econ. Rev.; 1958; 48, pp. 50-72.

4. Samuelson, P.A. Rational Theory of Warrant Pricing. Ind. Manag. Rev.; 1965; 6, pp. 13-31.

5. Schwartz, E. The Stochastic Behavior of Commodity Prices Implications for Valuation and Hedging. J. Financ.; 1997; 52, pp. 925-973. [DOI: https://dx.doi.org/10.1111/j.1540-6261.1997.tb02721.x]

6. Cartea, A.; Figueroa, M.G. Pricing in Electricity Markets: A Mean-Reverting Jump-Diffusion Model with Seasonality. Appl. Math. Financ.; 2005; 12, pp. 313-335. [DOI: https://dx.doi.org/10.1080/13504860500117503]

7. Borovkova, S.; Geman, H. Seasonal and stochastic effects in commodity forward curves. Rev. Deriv. Res.; 2006; 9, pp. 167-186. [DOI: https://dx.doi.org/10.1007/s11147-007-9008-4]

8. Frau, C.; Fanelli, V. Seasonality in commodity prices: New approaches for pricing plain vanilla options. Ann. Oper. Res.; 2024; 336, pp. 1089-1131. [DOI: https://dx.doi.org/10.1007/s10479-022-05128-x]

9. Merton, R.C. Option pricing when underlying stock returns are discontinuous. J. Financ. Econ.; 1976; 3, pp. 125-144. [DOI: https://dx.doi.org/10.1016/0304-405X(76)90022-2]

10. Hilliard, J.E.; Reis, J. Valuation of Commodity Futures and Options under Stochastic Convenience Yields, Interest Rates, and Jump Diffusions in the Spot. J. Financ. Quant. Anal.; 1998; 33, pp. 61-86. [DOI: https://dx.doi.org/10.2307/2331378]

11. Fileccia, G.; Sgarra, C. A particle filtering approach to oil futures price calibration and forecasting. J. Commod.; 2018; 9, pp. 21-34. [DOI: https://dx.doi.org/10.1016/j.jcomm.2017.12.003]

12. Nguyen, D.B.B.; Prokopczuk, M. Jumps in commodity markets. J. Commod. Mark.; 2019; 13, pp. 55-70. [DOI: https://dx.doi.org/10.1016/j.jcomm.2018.10.002]

13. Brignone, R.; Gonzato, L.; Sgarra, C. Commodity Asian option pricing and simulation in a 4-factor model with jump clusters. Ann. Oper. Res.; 2024; 336, pp. 275-306. [DOI: https://dx.doi.org/10.1007/s10479-022-05152-x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36643858]

14. Fileccia, G.; Sgarra, C. Historical and risk-neutral estimation in a two factors stochastic volatility model for oil markets. Int. J. Comput. Econ. Econom.; 2015; 5, pp. 451-479. [DOI: https://dx.doi.org/10.1504/IJCEE.2015.072289]

15. Mason, C.F.; Wilmot, N.A. Jumps in the convenience yield of crude oil. Resour. Energy Econ.; 2020; 60, 101163. [DOI: https://dx.doi.org/10.1016/j.reseneeco.2020.101163]

16. Mirantes, A.G.; Población, J.; Serna, G. The stochastic seasonal behavior of energy commodity convenience yields. Energy Econ.; 2013; 40, pp. 155-166. [DOI: https://dx.doi.org/10.1016/j.eneco.2013.06.011]

17. Gonzato, L.; Sgarra, C. Self-exciting jumps in the oil market: Bayesian estimation and dynamic hedging. Energy Econ.; 2021; 99, 105279. [DOI: https://dx.doi.org/10.1016/j.eneco.2021.105279]

18. Hull, J. Options Futures and Other Derivatives; Pearson Education: London, UK, 2018.

19. Geman, H. Commodities and Commodity Derivatives: Modeling and Pricing for Agriculturals, Metals and Energy; John Wiley & Sons: Hoboken, NJ, USA, 2005.

20. Glasserman, P. Monte Carlo Methods in Financial Engineering; Springer: New York, NY, USA, 2003.

21. Saltelli, A.; Tarantola, S.; Campolongo, F. Sensitivity Analysis as an Ingredient of Modeling. Stat. Sci.; 2000; 15, pp. 377-395.

22. Rujivan, S. Analytically pricing variance swaps in commodity derivative markets under stochastic convenience yields. Commun. Math. Sci.; 2021; 19, pp. 111-146. [DOI: https://dx.doi.org/10.4310/CMS.2021.v19.n1.a5]

23. Cortazar, G.; Schwartz, E.S. Implementing a stochastic model for oil futures prices. Energy Econ.; 2003; 25, pp. 215-238. [DOI: https://dx.doi.org/10.1016/S0140-9883(02)00096-8]

24. Ross, S.M. Introduction to Probability and Statistics for Engineers and Scientists; 4th ed. Elsevier: Amsterdam, The Netherlands, 2014.

25. Rujivan, S.; Thamrongrat, N.; Jantanon, P.; Djehiche, B. Analytical computation of conditional moments in the extended Cox–Ingersoll–Ross process with regime switching: Hybrid PDE system solutions with financial applications. Math. Comput. Simul.; 2025; 229, pp. 176-202. [DOI: https://dx.doi.org/10.1016/j.matcom.2024.09.032]