Estimating the Competitive Storage Model with Stochastic Trends in Commodity Prices

We propose a state-space model (SSM) for commodity prices that combines the competitive storage model with a stochastic trend. This approach fits into the economic rationality of storage decisions, and adds to previous deterministic trend specifications of the storage model. Parameters are estimated using a particle Markov chain Monte Carlo procedure. Empirical application to four commodity markets shows that the stochastic trend SSM is favored over deterministic trend specifications. The stochastic trend SSM identifies structural parameters that differ from those for deterministic trend specifications. In particular, the estimated price elasticities of demand are significantly larger under the stochastic trend SSM.


Introduction
Economic theories are often developed in a stationary context. However, the real world does not always correspond to stationarity. This potential mismatch creates a challenge when attempting to relate theory to historical data. This is a well-known problem in empirical macroeconomics, where structural parameters of business cycle models are often estimated on data that have been filtered in order to remove variation at frequencies that the model is not intended to explain, such as low-frequency trend variations and seasonal fluctuations (DeJong and Dave, 2011;Sala, 2015). For an overview of alternatives to the use of pre-filtered data in order to address this general problem, see Canova (2014).
In the competitive storage model for commodity prices introduced by Gustafson (1958), the situation is similar to that of business cycle models. The rational expectations equilibrium implied by the solution of this model is only known to exist in a stationary market. Accordingly, it is a model for describing dynamic price adjustments towards an exogenously given fixed steady-state equilibrium. However, it cannot explain low-frequency price movements due to persistent shocks. This is problematic when attempting to estimate the structural parameters of the model using commodity price data, since time series of commodity prices typically display a strongly persistent behavior in the price level, so that non-stationarity cannot be rejected when using conventional statistical tests (Wang and Tomek, 2007;Gouel and Legrand, 2017). As a result, the estimates for the structural parameters, which determine quantities like the price elasticity of demand and storage costs, are likely to be biased. This issue was recognized by Deaton and Laroque (1995) in one of the earliest attempts to directly estimate the structural parameters of the storage model. This paper proposes an approach to estimate the structural parameters of the competitive commodity storage model using a state-space model (SSM) for commodity prices, which decomposes the observed price into a stationary component which is due to the storage model and a stochastic trend component included to capture low-frequency price variations the storage model is unable to explain. Using a stochastic trend specification to account for non-stationary price data, our empirical approach aims at fitting into the economic rationality of the stationary storage model so that it preserves theoretical coherence, promising meaningful estimates of the structural parameters. Such a fit results from the fact that a stochastic trend that scales equilibrium prices can be isolated in the storage model by assuming that the innovations to the trend do not interfere with the agents' equilibrium storage decisions. In the baseline storage model, unrestricted equilibrium storage decisions lead to an intertemporal pricing restriction of the form P t = βE t (P t+1 ), where E t (P t+1 ) is the rational period-t expectation of the commodity price P t+1 and β represents some discount factor. Thus, a stochastic price scaling K t will not impair the equilibrium storage decisions if K t P t = βE t (K t+1 P t+1 ). This generically identifies stochastic trends as shifts in the price levels that do not interfere with intertemporal stock allocations, allowing a coherent integration of the stationary rational expectations equilibrium into a non-stationary environment, thus providing the theoretical basis of our empirical SSM approach. The corresponding SSM, that jointly identifies the trend parameters and the structural parameters of the storage model, is non-linear in the latent states so that its likelihood function is not available in closed form. To overcome this difficulty, we propose to use a Bayesian posterior analysis based on a particle Markov chain Monte Carlo (PMCMC) procedure (Andrieu et al., 2010).
With our proposed approach we contribute to the literature concerned with the general problem of adapting stationary economic models to non-stationary data, and more specifically to the problem of estimating the structural parameters of the competitive storage model on non-stationary commodity price data. Legrand (2019) identifies reliable estimation as one of the main issues of structural models for commodity prices.
Early attempts of estimating the structural parameters revealed that fitted competitive storage models are not able to satisfactorily approximate the observed strong serial dependence in commodity price data, indicating misspecification of the empirical model and casting doubt on the reliability of the parameter estimates (Deaton and Laroque, 1995). Suggested solutions to this problem include ad-hoc enrichments of the dynamic structure of the storage model by including weakly dependent supply shocks (Deaton and Laroque, 1996;Kleppe and Oglend, 2017), or the tuning of the grid for the commodity stock state variable, used for approximating the policy function . Other approaches replace the estimation techniques applied in early empirical implementations of the storage model, like the pseudo maximum likelihood (ML) procedure of Deaton and Laroque (1996), by more sophisticated ones, such as the ML technique developed by Cafiero et al. (2015) or the particle filtering methods proposed in Kleppe and Oglend (2017).
Empirical approaches that, like ours, decompose the observed price into a component to be explained by the storage model and a trend component are those of , Bobenrieth et al. (2013), Guerra et al. (2015) and Gouel and Legrand (2017). The first three of these studies propose to account for the strong persistence in the price data that the storage model is not able to approximate, by detrending the prices using a deterministic log-linear trend prior to the estimation of the structural parameters. Gouel and Legrand (2017) improves upon this procedure by jointly estimating the structural and deterministic trend parameters using the ML-estimator of Cafiero et al. (2015). The trend specifications Gouel and Legrand (2017) consider in their empirical application include log-linear trends as well as more flexible trends specified as restricted cubic splines. One of their main findings is that empirical models accounting for a properly specified trend component in the observed commodity price yield more plausible estimates of the structural parameters than models without a trend. However, the deterministic trends used in those studies inherently imply well predictable capital gains in the storage model, and so question the economic logic of separating the trend from structural economic pricing components. Moreover, the appropriate functional form of the deterministic trend needs to be tailored to the specific commodity market and the sampling frequency for which the storage models are applied. In contrast, the stochastic trend as used in our SSM approach represents, in Bayesian terms, a hierarchical prior for the low-frequency price component, which is not only consistent with the rationality of the economic model, but also flexible in its design to account for variation that the storage model is not intended to explain. This makes our approach applicable to a broad range of commodity markets and different sampling frequencies. The strategy of scaling prices to address nonstationarity was also done by Routledge et al. (2000) in their equilibrium term structure model of crude oil futures. However, they did not do so in a rigorous estimation framework.
A stochastic trend as used in our storage SSM allows a potentially large fraction of the observed variation in commodity prices to be accounted for by the trend component. This risks miss-assigning price variation due to speculative storage to the trend component. Thus, if considered as an evaluation of the empirical relevance of the storage model, the use of a stochastic trend can be considered as a conservative test. To explore this issue further we perform a simulation experiment. The simulation results suggest that our proposed approach is able to accurately assign price variation to trend and model components. We further apply our storage SSM to monthly observations of nominal coffee, cotton, aluminum and natural gas prices.
The results show, not surprisingly, that most of the observed price variation is due to the stochastic trend component. In order to assess the empirical relevance of the competitive storage model, we compare the storage SSM to the nested model that results in the absence of storage. The comparison reveals that the storage model predicting non-linear price dynamics with episodes of isolated price spikes and increased volatility adds significantly to explaining the observed commodity price behavior. We also compare the stochastic trend SSM to the deterministic trend models of Gouel and Legrand (2017) by using the Bayes factor and a model residual analysis. Results show that the SSM with stochastic trend fits the price data much better than models with deterministic trends. The estimates for the price elasticity of demand obtained from the stochastic trend SSM are substantially larger than for the deterministic trend models. Also, the estimated storage costs vary considerably depending on the commodity. This highlights the importance of properly accounting for the trend behavior when evaluating the role of speculative storage in commodity markets.
The rest of this paper is structured as follows. In the next section, we present the storage model used in the paper and the assumed price representation. We then present the estimation methodology (Section 3), simulation results (Section 4) and empirical results for historical data (Section 5). We discuss the findings before we offer some concluding remarks (Section 6).

State-Space Formulation with Stochastic Trend
Our approach relies upon the commodity storage model of Oglend and Kleppe (2017). It extends the Deaton and Laroque (1992) model by including an upper limit of storage capacity, C ≥ 0, in addition to the conventional non-negativity constraint for stocks, so that the storage space is completely bounded. This upper limit takes into account possible congestion of the storage infrastructure, which can lead to negative price spikes in the event of substantial oversupply in the market. In addition, the assumption of a completely bounded storage space allows numerical solutions of the model that are more robust over a wider parameter range than those for a model without this assumption . This, in turn, simplifies estimation of the model parameters.
The economic model of commodity storage is a canonical dynamic stochastic partial equilibrium model in discrete time for a commodity market with risk neural storage agents and rational expectations. The rational expectations equilibrium is characterized by a price function, denoted by f (x), which maps the stocks x to commodity prices. For empirical implementation, we assume that the observed commodity price can be decomposed into a component to be explained by the commodity storage model and a stochastic trend component. The corresponding time series model that we propose for the commodity log-price p t , observed at time t (t = 1, . . . , T ), has the form where the available quantity of commodity stocks x t is treated as a latent state variable. Its dynamics are linear in the equilibrium storage policy σ(·), with stock depreciation rate δ and Gaussian supply shocks z t .
The latent trend component of the log-price k t is specified as a driftless Gaussian random walk, so that it is allowed to vary gradually over time. The innovations of this stochastic trend ε t and the supply shocks z t are assumed to be serially and mutually independent.
The rational expectations equilibrium price function f (x) satisfies for all x where D(p) represents a continuous and monotonically decreasing aggregate demand function in the market, P (x) is the corresponding inverse demand, and φ(z) is the probability density function of the supply shock z.
The storage cost discount factor is given by β = (1−δ)/(1+r), where r is a relevant interest rate. According to Equation (4), the equilibrium pricing function exhibits three different pricing regimes: (i) a stock-out pricing wheref (x) is the expected next period commodity price, and (iii) a full capacity pricing regime, where The stock-out regime is characterized by positive price spiking and high price volatility due to reduced shock buffering capabilities in the market. Under the no-arbitrage regime, prices evolve smoothly with a relatively low volatility. Full capacity pricing mirrors the stock-out regime but with negative price spikes. As the market transitions between regimes, prices move between periods of quiet and turmoil, generating non-linear dynamics in the price process. The rational expectations equilibrium f (x) is stationary, having an associated globally stationary price density .
Using k t = p t−1 − log f (x t−1 ) + ε t in the price equation, the model as given in Equations (1)-(3) can be written as This defines a non-linear Gaussian state-space model, with measurement equation (7) for the observed price and state-transition equation (8) for the latent stocks.

Stochastic Trends and Storage Decisions
Separating the trend from the storage model pricing component in a consistent way that does not compromise the rationality of storage agents in the market requires that the trend does not interfere with intertemporal allocation incentives. The martingale property of a trend component specified as a stochastic trend with innovations that are independent of supply shocks ensures that this requirement is met. By using a separable stochastic trend we are assuming that storage agents do not alter their storage decisions based on trend innovations. In other words, trend innovations are assumed perceived by agents as permanent scalings of price levels that do not warrant adjustments to storage allocations.
As an example, consider permanent shocks K to the inverse aggregate demand in the market, P * = KP .
The aggregate demand implied by P * is D * , and the resulting rational expectations equilibrium is given by Assume the scaling process is given by K = γK + , where is a random variable with density φ which is independent of the supply shock z. This scaling does not affect the optimal storage policy if f * (x) = Kf (x) solves the functional equation problem in Equation (9), where f (x) is the rational expectations equilibrium for the original non-scaled prices. Substituting the proposed solution for f * (x), we get Note that D * (Kf (x)) = D(f (x)) by the definition of D * as the inverse of the scaled inverse demand function . And so, which establishes f * (x) = Kf (x) as the solution to the inverse demand scaled rational expectations equilibrium. Consequently, any observed commodity price can be represented as P = Kf (x). Formally, the rational expectations equilibrium is linear homogeneous to the proportional scaling K of the inverse aggregate demand function when E(K ) = K and innovations to the trend are orthogonal to supply shocks.
Note that in our econometric model as given by Equations (1)-(3), it is the logarithm of the scaling process for the price levels k = log(K) and not K, for which we assume a stochastic trend. Hence, the martingale property for the scaling term process K will not apply exactly, and the assumed Gaussian process for log(K) implies that E(K ) = K exp v 2 /2 > K. By ignoring this bias in our econometric model we make the behavioral assumption that agents do not alter storage decision based on the capital gain due to the expected mark-up factor exp v 2 /2 > 1. We consider this a reasonable trade-off to allow us to empirically analyze the storage model within a log-linear state and measurement space, which is comparatively convenient for statistical inference. In addition, the bias is small when v 2 is small, that is, when the trend is fairly smooth. In fact, the estimates we obtain for v in our empirical application discussed below imply that the factor exp v 2 /2 varies in a range between 1.001 and 1.004 so that it is essentially negligible. Ignoring this factor is essentially equivalent to transforming the probability space to a setting where agents ignore information from trend innovations, similar to a risk-neutral valuation setting where v 2 defines a required risk premium term or a nominal inflation term.

Preliminaries
In our empirical application of the storage model with stochastic trend based on its state-space representation as given in Equations (7) and (8), we use monthly commodity spot prices and rely on a Bayesian Markov chain Monte Carlo (MCMC) posterior analysis. For this application, we follow Kleppe and Oglend (2017) and use P (x) = exp(−bx) as the inverse demand function, where the parameter b measures the semi-elasticity of the demand price. In line with Gouel and Legrand (2017), we fix the yearly interest rate at 5%, so that the monthly storage cost discount factor is given by β = (1 − δ)/(1 + r), with r = 1.05 1/12 − 1. The set of parameters then consists of the structural parameters (δ, b, C) and the trend parameter v.
In initial experiments to estimate the parameters, we found that the capacity limit C is empirically not well identified separately from the remaining parameters. This appears to be mainly due to the fairly small sample size of our data, ranging from 264 to 360 monthly spot price observations. Therefore, we decided to fix C at a positive predetermined value. Since C determines the full capacity threshold for equilibrium storage, it bounds the space for the unit-free latent state variable x, and by fixing its value (together with normalizing the mean of z to zero) we pin down the range of this space. The values of the remaining parameters (δ, b, v) and their implications for the price dynamics are then to be interpreted relative to this scale of x. In our empirical application below we set the capacity limit C = 10. This ensures that it is a fairly rare event for the market to be in the full-capacity regime. Suppose, for example, that all realizations of the supply shocks z for a sequence of periods are equal to one standard deviation and that nothing is consumed in those periods, so that x t+1 = (1 − δ)x t + 1. Then a storage infrastructure with C = 10 and a notable depreciation rate of δ = 0.01 can store those unconsumed supplies for about 10 months before reaching the capacity limit 1 .
In order to solve the functional equation for the equilibrium price function f (x) as defined by Equations (4)-(6), we use a numerical algorithm which is based on the method of Kleppe and Oglend (2019), detailed in Appendix A.1. This algorithm takes advantage of the fact that the storage space is completely bounded by the non-negative constraint and the capacity limit C, thus providing numerically robust and computationally fast solutions. This is critical for a Bayesian MCMC posterior analysis because it requires a significant number of reruns of the algorithm to obtain a solution to the pricing function for each new parameter value.

Bayesian Inference Using Particle Markov Chain Monte Carlo
In the SSM model as given by Equations (7) and (8) the vector of parameters to be estimated is given by θ = (v, δ, b) and the vector of latent state variables is x 1:T , where the notation a s:s is used to denote (a s , a s+1 , . . . , a s ). The posterior of the parameters is π(θ|p 1:T ) ∝ π θ (p 1:T )π(θ), where π(θ) denotes the prior density assigned to θ and π θ (p 1:T ) represents the likelihood function, given by where N (·|µ, σ 2 ) denotes a normal density function with mean µ and variance σ 2 . For the joint density of the price and state in the initial period π θ (p 1 , x 1 ) we assume that it factorizes into a uniform density on , and a dirac measure for the price p 1 located at its actually observed value (effectively conditioning the likelihood on the first price observation).
Due to the non-linear nature of the pricing function f (x) and the storage function σ(x) entering the measurement and state transition density as given in Equation (14), the likelihood (and hence the resulting posterior for θ) are not available in closed form, so that a Bayesian and likelihood-based inference requires approximation techniques. Several Monte Carlo (MC) approximation approaches have been developed for statistical inference in non-linear SSMs with analytically intractable likelihood functions. However, only a few of them are suited to the model considered here due to the discontinuous derivatives of f (x). In particular, methods using MC estimators for the likelihood π θ (p 1:T ) based on approximations to the conditional posterior of the states π(x 1:T |θ, p 1:T ), including second order/Laplace approximations (Shephard and Pitt, 1997;Durbin and Koopman, 2012) or global approximations as used by the efficient importance sampler (Liesenfeld and Richard, 2003;Richard and Zhang, 2007), perform poorly in such a context. The same applies to the Gibbs approach targeting the joint posterior distribution of the states and parameters π(x 1:T , θ|p 1:T ) and alternately simulating from the conditional posteriors π(x 1:T |θ, p 1:T ) and π(θ|x 1:T , p 1:T ).
It is known that such a Gibbs procedure typically has problems in efficiently approximating the targeted joint posterior in non-linear SSMs due to a fairly slow mixing (Bos and Shephard, 2006). Moreover, the Gibbs procedure is also not very computationally attractive in the present context, since both the (joint) conditional posterior of all the states π(x 1:T |θ, p 1:T ) and the single-site conditional posterior of the individual states π(x t |x 1:t−1 , x t+1:T , θ, p 1:T ) are non-standard distributions.
Here, we propose to use the particle marginal Metropolis-Hastings (PMMH) approach as developed by Andrieu et al. (2010), which is well suited for a posterior analysis of our proposed storage SSM as it can cope with the discontinuity of the gradients of f (x) and is very easy to implement. The PMMH uses unbiased MC estimates of the likelihood π θ (p 1:T ) inside a standard Metropolis-Hastings (MH) algorithm targeting the posterior of the parameters π(θ|p 1:T ). The MC estimation error of the likelihood estimate does not affect the invariant distribution of the MH so that the PMMH allows for exact inference. The PMMH produces an MCMC sample {θ i } S i=1 from the target distribution by the following MH updating scheme: Given the previously sampled θ i−1 and the corresponding likelihood estimateπ θi−1 (p 1:T ), a candidate value θ * is drawn from a proposal density Q(θ|θ i−1 ), and the estimate of the associated likelihoodπ θ * (p 1:T ) is computed. Then the candidate θ * is accepted as the next simulated θ i with probability otherwise θ i is set equal to θ i−1 . Under weak regularity conditions, the resulting sequence {θ i } S i=1 converges to samples from the target density π(θ|p 1:T ) as S → ∞ (Andrieu et al., 2010, Theorem 4).
For the PMMH, we use a Gaussian random walk proposal density Q(θ|θ i−1 ) = N (θ|θ i−1 , Σ) and follow the approach of Haario et al. (2001) to adaptively set the proposal covariance matrix Σ during the burn-in period of the MCMC iterations. After dropping the draws from the burn-in period, we use the θ draws from the next M PMMH iterations to represent the posterior π(θ|x 1:T ). The posterior mean of the parameters, used as point estimates, is approximated by the sample mean over the M PMMH draws. For numerical stability of the PMMH computations, we reparameterize the likelihood function using the transformed parameters θ = log(v), arctanh(2δ − 1), log(b) so that the resulting parameter space is unconstrained.
The prior densities for the parameters are selected as follows: For log(b) we assume a N (0, 1) prior, and for v 2 an inverted chi-squared prior with v 2 ∼ 0.1/χ 2 (10) , where χ 2 (10) denotes a chi-squared distribution with 10 degrees of freedom. Under this prior for v 2 , the mean is given by 0.01 and the standard deviation by 0.007. The prior density assigned to δ is a Beta with δ ∼ B(2, 20) so that the mean and standard deviation are given by 0.09 and 0.05, respectively.

Particle Filter Likelihood Evaluation
In order to obtain unbiased MC estimates for the likelihood in Equation (13), required as an input of the PMMH, we follow Andrieu et al. (2010) and Flury and Shephard (2011) and use a simple sampling importance resampling (SIR) particle filter (PF). For given values of the parameters θ, it produces MC estimates for the sequence of period-t likelihood contributions π θ (p t |p 1:t−1 ) by sequentially sampling and resampling using an importance sampling (IS) density q(x t |x 1:t−1 ) for the states x t (see, Doucet and Johansen, 2009;Cappé et al., 2007, for a detailed treatment of PFs). For the implementation of the PF we use the state-transition density π θ (x t |x t−1 ) as IS density (Gordon et al., 1993), and rely on a dynamic resampling scheme in which the particles are resampled only when their effective sample size falls below one half of the number of particles (Doucet and Johansen, 2009). This simple version of the PF (also known as the bootstrap PF, BPF) for approximating the likelihood as given by Equations (13) and (14) consists of the following steps: . Compute the IS weights as and their normalized versions Then use the IS weights to obtain the period-t likelihood contribution asπ θ (p t |p 1:t−1 ) = ( N k=1 w k t )/N , and compute the effective particle sample size defined by If N e t < N/2, resample from the particles {x k 1:t } N k=1 with replacement according to their IS weights W k t to obtain the resampled particles {x k 1:t } N k=1 , and set their weights to W k t = 1/N .
The resulting BPF estimate for the likelihood (conditional on the first price observation) is given bŷ π θ (p 1:T ) = [ T t=2π θ (p t |p 1:t−1 )]. The measurement density π θ (p t |p t−1 , x t−1:t ) is not very informative about the states x t for empirically relevant parameter values, resulting in a low signal-to-noise ratio. Thus, the simple BPF yields fairly precise MC estimates of the likelihood with a modest number of particles N (Cappé et al., 2007). High precision likelihood estimates are a critical requirement for the PMMH to produce a well mixing MCMC sample from the posterior of the parameters π(θ|x 1:T ) (Flury and Shephard, 2011). In our applications below, we use N = 10, 000 particles. For a time series with T = 360, one BPF likelihood estimate requires approximately 2.5 seconds (on a computer with an Intel Core i5-6500 processor running at 3.20 GHz). The MC numerical standard deviation of the log-likelihood estimate logπ θ (p 1:T ), computed from reruns of the BPF for a fixed θ value under different seeds, is about 0.1 percent of the absolute value of the log-likelihood, illustrating the high accuracy of the BPF.

State Prediction and Diagnostics
The BPF outlined in the previous section, and used for the PMMH implementation, can also be used to produce MC estimates for the predicted values of the latent state vector x 1:t+1 and functions thereof, given the prices observed up to period t, p 1:t . MC estimates of such predictions can serve as the basis for diagnostic checks. Let h(x 1:t+1 ) be a function of interest in x 1:t+1 . Its conditional mean given p 1:t can be expressed as where π θ (x t+1 |x t ) is the state-transition density as given by Equation (14) and π θ (x 1:t |p 1:t ) is the filtering density for x 1:t . Since the particles and IS-weights {x k 1:t , W k t } N k=1 produced by the BPF provide an MC approximation to this filtering density, the conditional mean in Equation (17) for a given value of θ can be easily estimated byÊ , where x k t+1 is obtained by propagating the BPF particle x k 1:t via the state-transition density, i.e. x k t+1 ∼ π θ (x t+1 |x k t ). In practice, the parameters θ are set equal to their estimates.
This MC approximation of a predicted mean like that in Equation (17) If the model is correctly specified, then η t+1 and η 2 t+1 are serially uncorrelated so that they can be used for diagnostic checking of the assumed dynamic structure. The conditional moments of p t+1 for the storage SSM are given by E(p t+1 |p 1:t ) = p t +E(log[f (x t+1 )/f (x t )]|p 1:t ) and Var(p t+1 |p 1:t ) = Var(log[f (x t+1 )/f (x t )]|p 1:t )+ v 2 , which can be evaluated by Equation (18), using the functions h( In order to check the capability of the storage SSM to approximate the distributional properties of the observed prices we use the probability integral transformed (PIT) residuals defined as where Pr(p t+1 ≤ p o t+1 |p 1:t ) is the predicted probability that p t+1 is less or equal to the actually ex-post observed price p o t+1 , and Φ denotes the cdf of a N (0, 1)-distribution (Kim et al., 1998). The PIT residuals ξ t+1 follow a N (0, 1)-distribution if the model is valid. For the storage SSM, the probability u t+1 can be calculated by setting the function h(

Marginal Likelihood for Model Comparison
Marginal likelihood is used to compare the storage SSM with alternative models and assess the empirical relevance of the structural storage model component in the SSM. In order to evaluate the marginal likelihood for the storage SSM we rely upon the procedure proposed by Chib and Jeliazkov (2001), which is specifically customized for Bayesian analyses implemented using MH algorithms targeting the posterior of the parameters. This procedure takes advantage of the fact that the marginal likelihood can be expressed as where πθ(p 1:T ) is the likelihood function for the observed prices evaluated at some value of the parametersθ, and π(θ) and π(θ|p 1:T ) are the corresponding ordinates of the prior and posterior of the parameters. Then it exploits that the posterior ordinate π(θ|p 1:T ) can be expressed in terms of the MH acceptance probability α(·, ·) and the proposal density Q(·|·). Namely, as the ratio of the expectation of α(θ, θ)Q(θ|θ) under the posterior π(θ|p 1:T ) relative to the expectation of α(θ,θ) under the proposal density Q(θ|θ). This implies that a consistent MC estimate for π(θ|p 1:T ) based on the MH acceptance probability defined in Equation (15) is given byπ where {θ i } M i=1 are the M simulated draws from the posterior distribution π(θ|p 1:T ) and {θ l } L l=1 are draws from the proposal distribution Q(θ|θ). For evaluating the likelihood πθ(p 1:T ) in Equation (21) we use the same BPF algorithm as applied for the computation of the MH acceptance probabilities in Equation (15) (and outlined in Section 3.3). The value of the pointθ is set equal to the posterior mean of θ.

Ability to Isolate the Trend and Storage Model Component
In order to illustrate the capability of our Bayesian storage SSM approach to empirically separate the variation in the observed prices into the variation generated by the structural storage model component and that of the stochastic trend, we conduct a simulation experiment. Prices are simulated from the storage SSM for parameters that are set equal to their posterior mean values found for the empirical application to natural gas prices, discussed further below (see Table 1). Prices are simulated for 800 periods, with the first 500 discarded as burn-in, so that the size of the simulated sample is T = 300. The storage SSM is then fitted to the time series of simulated prices by using the PMMH procedure, and the BPF is applied to produce estimates of the filtered mean for the storage model price component E(f (x t )|p 1:t ) and the stochastic trend E(k t |p 1:t ), evaluated at the posterior mean of the parameters.  correspond to the stock-out pricing regime, and values below the lower boundary correspond to the full storage capacity pricing regime. The plotted time series show that the estimated filtered means of the two price components track the time evolution of the true components quite well. We also observe that the estimated filtered means predict most of the stock-out events, which is critically important for identifying the structural parameters of the storage model. These simulation results illustrate that our approach appears to be well capable to empirically identify -from observed commodity prices -the fraction of the price variation which is due competitive storage decisions and to separate it from stochastic trend variation.

Empirical Application
In this section, we apply our Bayesian storage SSM approach to historical monthly price data for the following

Estimation Results for the Storage SSM with Stochastic Trend
For the Bayesian posterior analysis of the storage SSM, we run the PMMH algorithm for 12,000 iterations and discard the first 2,000 as burn-in. In order to evaluate the sampling efficiency of the PMMH for estimating the parameters, we compute the effective sample size (ESS) of their posterior PMMH samples (Geyer, 1992 monthly price changes for natural gas, 66% for coffee, 71% for cotton, and 81% for aluminum. As for the estimates of the depreciation rate δ, we observe that they are fully in line with the actual storage costs to be expected for the different types of commodities: For natural gas we find the largest estimated depreciation rate (1.1%), which implies that the monthly cost of storage amounts to 1.5% of the price. This relatively large estimated storage cost is in accordance with the fairly expensive storage technology for US natural gas, which is typically stored in underground salt caves and similar facilities. The second largest storage cost is found for coffee, with a monthly depreciation rate of 0.2% leading to estimated monthly costs of 0.6% of the price. The lowest storage costs are predicted for the non-food and non-energy products cotton and aluminum, for which the estimated depreciation rate is 0.1% resulting in storage costs of 0.5%. We also observe that the larger the estimated storage cost for a commodity, the larger the fraction of observed price variation which is captured by the storage decision behavior. This is in agreement with the rationality of the competitive storage model, where higher storage costs are associated with more frequent stock-out events, which in turn implies greater price volatility. The posterior mean values for the slope parameter b of the inverse demand function imply that a reduction in supply on the market by one standard deviation of production leads to a price increase of 42% for natural gas, 38% for coffee, 32% for cotton and 20% for aluminum. The size of these estimated price elasticities roughly corresponds to the size of the price peaks observed in these markets.  Table 1. peaks and drops. Beyond the periods with elevated price volatility, the contribution of this component to the price variation appears small. This reflects that when equilibrium storage is an inner solution (so that 0 < σ(x t ) < C), the resulting price is subject to an intertemporal price restriction leading to prices which behave as a stationary Markov process. Accordingly, in this no-arbitrage pricing regime, the economic storage model provides little additional information about the price evolution that goes beyond the stochastic trend. However, storage becomes empirically relevant with a significant impact on the price behavior when the normal no-arbitrage pricing mechanism collapses in the stock-out and full-capacity regime, which occurs in the storage model in periods of severe and prolonged commodity shortages or oversupply.
The limits-to-arbitrage regimes (stock-out or full-capacity) detected by the fitted storage model tend to coincide with known historical market events. For example, the time periods with peaks in the filtered storage price component for natural gas usually correspond to periods when the historical level of natural gas storage in the market was very low . The sharp drop in the storage price component for coffee in 1989 coincides with the collapse of the International Coffee Agreement (a cartel of coffee-producing countries) and oversupply in the market due to World Bank subsidies, while the 1994 peak is consistent with a negative supply shock triggered by significant frost damage in much of the coffee-growing areas of Brazil. The cotton price peak detected by the storage model in 2011 was arguably due to the severe global shortages, which were caused, inter alia, by the tightening of Indian export restrictions on cotton.
The early nineties spike in aluminum prices coincides with the collapse of the Soviet Union, and the 2008-2009 price drop is consistent with the sharp decline in global aluminum demand that created a large stock overhang during this period after the subprime crisis.

Model Comparisons
In this section, we assess the empirical relevance of the price component related to the competitive storage model for explaining the observed price variation, and compare the storage SSM model with stochastic trend to that with deterministic trend specifications. For this assessment, we rely on the marginal likelihood as well as diagnostic checks on Pearson and PIT residuals.

Alternative Models
For assessing the relevance of the storage model price component, we compare our SSM model to the restricted SSM that results in the absence of storage. The latter is obtained by letting δ → 1, making storage prohibitively costly, so that the stock process x t collapses to that of the supply shocks z t . In this case the SSM in Equations (1)-(3) with the assumed demand function P (x) = exp(−bx) reduces to This represents a standard linear Gaussian local level (LGLL) SSM (Durbin and Koopman, 2012)   same as those we assume for the unrestricted storage SSM.
As deterministic trend specifications to be compared with the stochastic trend in the storage SSM, we consider those used in the study of Gouel and Legrand (2017). They use a linear time trend, for which k t in Equation (2) is replaced by k t = α + βt. In addition, they consider restricted cubic spline trend are the basis functions of B-splines, G is the degree of freedom, and γ g are the corresponding trend parameters to be estimated. For our comparison we consider restricted cubic splines with 3 knots (RSC3) and 5 trend parameters as well as 7 knots (RSC7) and 9 trend parameters 2 . For these deterministic trends the SSM in Equations (1)-(3) reduces to a univariate, non-linear autoregression for the log-price: Analogously to the LGLL SSM, we can simulate from the posterior for the parameters of the deterministic trend models by using a standard MH algorithm. For the structural parameters (δ, b) we assume the same priors as used in the storage SSM, and to the deterministic trend parameters (α, β, γ g ) we assign independent N (0, 20 2 ) priors. For details on the computation and derivation of the Pearson and PIT residuals of the deterministic trend models, see Appendix A.2.

Marginal Likelihood Model Comparisons and Diagnostics Checks
Table 2 provides the log marginal likelihood values log π(p 1:T |model s ) for the storage SSM together with those of the LGLL SSM and the storage model combined with the deterministic trend specifications. Also reported are the resulting values for the log Bayes factor of the storage SSM relative to the four alternative models t Natural gas  log[π(p 1:T |storage SSM)/π(p 1:T |model )]. The results reveal that the storage SSM is strongly preferred over the LGLL SSM for all commodities, which suggests that the structural storage component in the SSM substantially contributes to the model fit. Hence, the non-linear price dynamics with periodically recurring increases in price volatility and price spiking, as predicted by the competitive storage model, adds significantly to explaining the price behavior. For all commodities, we also observe that the storage SSM is clearly favored over all deterministic trend specifications. Thus, the storage SSM has a trend component that is not only consistent with the rationality of the economic model, but is also much more supported by the data than the deterministic trends, such as those used by Gouel and Legrand (2017) for the estimation of the structural parameters of the competitive storage model. Our estimates of the structural parameters for the deterministic trend models are found in Appendix A.3. Figure 3 shows the time series plots of the fitted deterministic trendsk t and the smoothed mean of the stochastic trend E(k t |p 1:T ), all computed by setting the parameters to their posterior mean values 3 . Unsurprisingly, we find that the stochastic trend captures a substantially larger fraction of the observed price variations than the deterministic trends. Table 3 provides the results of diagnostic checks on the PIT residuals ξ t and the Pearson residuals η t for the storage SSM and the four alternative models considered. The PIT residuals of the storage SSM suggest that this model accounts well for the observed distributional properties of the prices for all commodities.
The skewness and kurtosis of its PIT residuals are close to their benchmark values for a normal distribution and they all pass the Jarque-Bera normality test at the 5% significance level. In contrast, the LGLL SSM  as well as the storage models with deterministic trends have difficulties approximating the distributional properties of the prices. Only the PIT residuals of the storage model with an RSC7 trend for aluminum pass the Jarque-Bera normality test at a conventional significance level.
The first-order serial correlation of the Pearson residuals η t and the p-values of the Ljung-Box test for η t and η 2 t including 12 lags reported in Table 3 show that the storage SSM successfully accounts for the observed autocorrelation in the level and volatility of the gas price, while they point towards significant residual correlation in price level and volatility for coffee, cotton and aluminum. However, all competing models cannot fully capture the serial correlation in the price levels of those three commodities either.
Only the volatility dynamics for coffee is better approximated by the linear and RSC7 trend model than by the storage SSM. Clearly, based on these results, we can not identify whether the failure of the storage SSM and the deterministic trend models to explain all of the observed dynamics in the coffee, cotton and aluminum prices is due to a potential misspecification of the trend or the competitive storage model itself, since the diagnostic tests are, as any specification test in this context, joint tests for the validity of both price components.
In sum, the results show that the storage SSM outperforms the deterministic trend models in explaining the observed distributional properties of commodity prices, and that its ability to account for the dynamics in the price levels is not worse. Only in the approximation of the volatility dynamics, the deterministic trend specifications appear to have a slight advantage.

Structural Parameter Estimates Under Stochastic and Deterministic Trends
As it is evident from Figure 6, the dynamic and distributional characteristics of the de-trended prices substantially differ depending on whether a stochastic or deterministic trend is assumed. Therefore, it can be expected that the nature of the trend has a critical impact on the estimates of the parameters that determine the storage costs (δ) and the price elasticity of demand (b), since these parameters are identified by the strength of the serial correlation and the size of the spikes in the trend-adjusted prices. The lower the storage costs in the competitive storage model are, the stronger the predicted serial correlation, while the more inelastic the demand is, the larger the resulting price spikes. As larger price spikes also imply more speculative storage activity, an inelastic demand also contributes to the strength of the predicted serial correlation in the prices.  Table 4: Estimates for the annual storage costs (net of interest costs) in percent of the average price and price elasticities of demand.
deterministic trend models used for annual data on various commodities, including coffee and cotton. This allows for some comparisons with our results for those two commodities. The annual storage costs estimates they report for their preferred trend model for coffee and cotton are, respectively, 1.4% and 0.3% of the average price. These estimates based on annual data are much lower than those we found for the storage SSM as well as the deterministic trend models fitted to monthly data. However, they argue that their estimated annual costs are possibly too small -an assessment that is consistent with our estimates for the storage costs. For the annual price elasticity of demand, the estimates of Gouel and Legrand (2017) are -0.04% for coffee and -0.03% for cotton. These estimates imply a demand for those commodities which is substantially more inelastic than that implied from our estimates. One can argue which elasticities better reflect the markets. Mehta and Chavas (2008) assume a range of plausible values for the annual elasticity of demand for coffee between -0.2% and -0.4%, while Duffy et al. (1990) argue that the annual export demand for cotton is likely fairly elastic. Hence, our elasticity estimates are more in line with these assessments than those found by Gouel and Legrand (2017).

Conclusion
In this paper, we have proposed a stochastic trend competitive storage model for commodity prices, which defines a non-linear state-space model (SSM). For the Bayesian posterior analysis of the proposed stochastic trend SSM, we use an efficient MCMC procedure. This adds to existing empirical commodity storage models based on deterministic trend specifications. Our stochastic trend approach fits into the economic rationality of the competitive storage model and is also sufficiently flexible to account for the variation in the observed prices that the competitive storage model is not intended to explain. The obvious benefit is that it makes the storage model applicable to markets with highly persistent unit root-like prices, which appears relevant for many commodity markets. Our approach aims at increasing the empirical relevance and applicability of the competitive storage model.
The MCMC procedure we propose for jointly estimating the structural and trend parameters in the SSM is a particle marginal Metropolis-Hastings algorithm based on the bootstrap particle filter. A Monte Carlo simulation experiment shows that this approach is able to disentangle the stochastic trend from the price variation due to speculative storage. The SSM is applied to monthly price data for natural gas, cotton, coffee and aluminum. Not surprisingly, the stochastic trend explains a large part of the observed variation in the commodity prices. More importantly, the competitive storage component adds short-run price volatility and price spiking, and becomes periodically relevant to explain non-linear pricing behavior related to states of market turmoil. A formal empirical comparison of the SSM to the corresponding model that results in the absence of storage suggests that the speculative storage price component significantly contributes to commodity price variation.
Which trend to apply will depend on the specific market under consideration. If a stochastic trend is not appropriate, fitting a highly flexible stochastic trend model risks overfitting the price variation and downplaying the contribution of the storage model. Consequently, the price elasticity of demand will tend to be overestimated and the estimates of the storage costs can be expected to be correspondingly biased. On the other hand, failing to account for a stochastic trend when it is appropriate will tend to underestimate the elasticity of demand. Our empirical results show that the stochastic trend model consistently estimates a higher elasticity of demand and a different amount of storage costs than existing deterministic trend models for the commodity markets investigated in this paper. Pre-testing of price characteristics can guide trend choice. For instance, unit root tests can be applied to evaluate whether a stochastic trend specification is suitable.
The empirical comparison of the stochastic trend SSM to existing deterministic trend models using the Bayes factor and model residual analysis shows that the stochastic trend fits the price data for the investigated commodity markets much better than the deterministic trends. In particular, in contrast to the deterministic trend specifications, the stochastic trend SSM captures the observed distributional properties of the prices, such as their skewness and kurtosis, quite well. While the stochastic trend SSM also accounts for the price dynamics in the observed prices on the natural gas market it is not able to fully capture all the serial dependence of the coffee, cotton and aluminum prices. This is similar to the results of financial approaches on modeling commodity term structures, showing the relevance of additional pricing factors beyond the traditional ones for the spot price and the convenience yield (Miltersen and Schwartz 1998;Schwartz 1997;Tang 2012). The stochastic trend SSM is essentially a two-factor model with one reduced-form random walk component orthogonally appended to a factor restricted by economic constraints. Increasing the flexibility in the economic model will arguably improve the explanatory power of the model, although with additional statistical challenges in separately identifying the trend behavior from the price component related to the competitive storage model.
6. Until convergence, set n ← n + 1 and go back to step 2.
Allowing the grid space to adjust to the updated functional solutions ensures that the grid can dynamically concentrate in the region of the state-space where a high precision is needed, namely the region defining the storage regime. This provides both efficient and precise numerical solutions to the pricing function.

A.2 Residuals for the Deterministic Trend Models
For the deterministic trend models as given by Equation (23)