CPM.MLE: N-Factor Model Parameter Estimation through the Kalman Filter...
In TomAspinall/TSE: N-Factor Commodity Pricing Models

Description Usage Arguments Details Value References Examples

\loadmathjax

the CPM.MLE function performs parameter estimation of an n-factor model given observable term structure data through maximum likelihood estimation and genetic algorithms. CPM.MLE allows for missing observations as well as constant or variable time to maturity of observed futures contracts.

CPM.MLE(
  log.futures,
  dt,
  TTM,
  N.factors,
  GBM = TRUE,
  S.Constant = TRUE,
  Estimate.Initial.State = FALSE,
  Richardsons.Extrapolation = TRUE,
  cluster = FALSE,
  Domains = NULL,
  ...
)

`log.futures`	Object of class `matrix` corresponding to the natural logarithm of observable futures prices. NA's are allowed within the `matrix`. Every column of the matrix must correspond to a particular futures contract, with each row corresponding to a quoted price on a given date.
`dt`	Constant, discrete time step of observations
`TTM`	Object of class 'vector' or 'matrix' that specifies the time to maturity of observed futures contracts. time to maturity can be either constant (ie. class 'vector') or time homogeneous (ie. class 'matrix'). When the time to maturity of observed futures contracts is time homogeneous, the dimensions of `TTM` must be identical to that of `log.futures`. Every element of `TTM` corresponds to the time to maturity, in years, of a futures contract at a given observation date.
`N.factors`	`numeric`. Number of state variables in the spot price process.
`GBM`	`logical`. When `TRUE`, factor 1 of the model is assumed to follow a Brownian Motion, inducing a unit-root in the spot price process.
`S.Constant`	`logical`. When `TRUE`, the white noise of observable contracts are assumed and identical (and independent).
`Estimate.Initial.State`	`logical`. When `TRUE`, the initial state vector is specified as unknown parameters of the commodity pricing model. When `FALSE`, a "diffuse" assumption is taken instead (see details)
`Richardsons.Extrapolation`	`logical`. When `TRUE`, the `grad` function from the `numDeriv` package is called to approximate the gradient within the `genoud` optimization algorithm.
`cluster`	an optional object of the 'cluster' class returned by one of the makeCluster commands in the `parallel` package to allow for parameter estimation to be performed across multiple cluster nodes.
`Domains`	an optional `matrix` of the lower and upper bounds for the parameter estimation process. The `CPM.bounds` function is highly recommended. When `Domains` is not specified, the standard bounds specified within the `CPM.bounds` function are used.
`...`	additional arguments to be passed into the `genoud` function. See `help(genoud)`

CPM.MLE is a wrapper function that uses the genetic algorithm optimization function genoud from the rgenoud package to optimize the log-likelihood score returned from the CPM.Kalman.filter function. When Richardsons.Extrapolation = T, gradients are approximated numerically within the optimization algorithm through the grad function from the numDeriv package. CPM.MLE is designed to perform parameter estimation as efficiently as possible, ensuring a global optimum is reached even with a large number of unknown parameters and state variables. Arguments passed to the genoud function can greatly influence estimated parameters and must be considered when performing parameter estimation. Recommended arguments to pass into the genoud function (see Examples for recommended arguments). All arguments of the genoud function may be passed through the CPM.Parameter.Estimation function, except for gradient.check, which is set to false.

CPM.MLE performs boundary constrained optimization of log-likelihood scores and does not allow does not allow for out-of-bounds evaluations within the genoud optimization process, preventing candidates from straying beyond the bounds provided by Domains. When Domains is not specified, the default bounds specified by the CPM.bounds function are used.

The N-factor model The N-factor model was first presented in the work of Cortazar and Naranjo (2006, equations 1-3). The N-factor framework describes the spot price process of a commodity as the correlated sum of \mjeqnNN state variables \mjeqnx_tx[t].

When GBM = TRUE: \mjdeqnlog(S_t) = \sum_i=1^N x_i,tlog(S[t]) = sum_i=1^n x[i,t] When GBM = FALSE: \mjdeqnlog(S_t) = E + \sum_i=1^N x_i,tlog(S[t]) = E + sum_i=1^n x[i,t]

Additional factors within the spot-price process are designed to result in additional flexibility, and possibly fit to the observable term structure, in the spot price process of a commodity. The fit of different N-factor models, represented by the log-likelihood can be directly compared with statistical testing possible through a chi-squared test.

Flexibility in the spot price under the N-factor framework allows the first factor to follow a Brownian Motion or Ornstein-Uhlenbeck process to induce a unit root. In general, an N-factor model where GBM = T allows for non-reversible behaviour within the price of a commodity, whilst GBM = F assumes that there is a long-run equilibrium that the commodity price will revert to in the long-term.

State variables are thus assumed to follow the following processes:

When GBM = TRUE: \mjdeqndx_1,t = \mu^*dt + \sigma_1 dw_1tdx[1,t] = mu^* dt + sigma[1] dw[1]t

When GBM = FALSE: \mjdeqndx_1,t = - (\lambda_1 + \kappa_1x_1,t)dt + \sigma_1 dw_1tdx[1,t] = - (lambda[1] + kappa[1] x[1,t]) dt + sigma[1] dw[t]t

And: \mjdeqndx_i,t =_i\neq 1 - (\lambda_i + \kappa_ix_i,t)dt + \sigma_i dw_itdx[i,t] =_(i != 1) - (lambda[i] + kappa[i] x[i,t]dt + sigma[i] dw[i]t)

where: \mjdeqnE(w_i)E(w_j) = \rho_i,jE(w[i])E(w[j])

The following constant parameters are defined as:

var \mjeqn\mumu: long-term growth rate of the Brownian Motion process.

var \mjeqnEE: Constant equilibrium level.

var \mjeqn\mu^*=\mu-\lambda_1mu^* = mu-lambda[1]: Long-term risk-neutral growth rate

var \mjeqn\lambda_ilambda[i]: Risk premium of state variable \mjeqnii.

var \mjeqn\kappa_ikappa[i]: Reversion rate of state variable \mjeqnii.

var \mjeqn\sigma_isigma[i]: Instantaneous volatility of state variable \mjeqnii.

var \mjeqn\rho_i,j \in [-1,1]rho[i,j] in [-1,1]: Instantaneous correlation between state variables \mjeqnii and \mjeqnjj.

Measurement Error:

The Kalman Filtering algorithm assumes a given measure of measurement error (ie. matrix \mjeqnHH). Measurement errors can be interpreted as error in the model's fit to observed prices, or as errors in the reporting of prices (Schwartz and Smith, 2000).

var \mjeqns_is[i] Observation error of contract \mjeqnii.

When S.Constant = T, the values of parameter \mjeqns_is[i] are assumed constant and equal to object 'sigma.contracts'. When S.Constant = F, the observation error of futures contracts \mjeqns_is[i] is equal to parameter 'sigma.contract_' [i].

Diffuse Kalman Filtering

If Estimate.Initial.State = F, a 'diffuse' assumption is used within the Kalman filtering algorithm. Factors that follow an Ornstein-Uhlenbeck are assumed to equal zero. When Estimate.Initial.State = F and GBM = T, the initial value of the first state variable is assumed to equal the first element of log.futures. This is an assumption that the initial estimate of the spot price is equal to the closest to maturity observed futures price.

The initial covariance of the state vector for the Kalman Filtering algorithm assumed to be equal to matrix \mjeqnQQ

Initial states of factors that follow an Ornstein-Uhlenbeck are generally not estimated with a high level of precision, due to the transient effect of the initial state vector on future observations, however the initial value of a random walk variable persists across observations (see Schwartz and Smith (2000) for more details).

CPM.MLE returns a list with 10 objects. 9 objects are returned when the user has specified not to calculate the hessian matrix at solution.

`MLE`	`numeric` The Maximum-Likelihood-Estimate of the solution
`Estimated.Parameters`	`vector`. The estimated parameters
`Standard.Errors`	`vector`. Standard error of the estimated parameters. Returned only when `hessian = T` is specified
`X.t`	`vector`. The final observation of the state vector
`X`	`matrix`. All observations of the state vector, after the updating equation has been applied
`Y`	`matrix`. Estimated futures prices at each observation
`V`	`matrix`. Estimation error of each futures contracts at each observation
`TSFit.Error`	`matrix`. The Mean Error (Bias), Mean Absolute Error, Standard Deviation of Error and Root Mean Squared Error (RMSE) of each observed contract, matching the column names of `log.futures`
`TSFit.Volatility`	`matrix`. The theoretical and empirical volatility of futures returns for each observed contract as returned from the `TSFit.Volatility` function
`genoud.value`	`list`. The output of the called `genoud` function.

Schwartz, E. S., and J. E. Smith, (2000). Short-Term Variations and Long-Term Dynamics in Commodity Prices. Manage. Sci., 46, 893-911.

Cortazar, G., and L. Naranjo, (2006). An N-factor Gaussian model of oil futures prices. Journal of Futures Markets: Futures, Options, and Other Derivative Products, 26(3), 243-268.

Mebane, W. R., and J. S. Sekhon, (2011). Genetic Optimization Using Derivatives: The rgenoud Package for R. Journal of Statistical Software, 42(11), 1-26. URL http://www.jstatsoft.org/v42/i11/.

##Example 1 - Perform One Generation of Maximum Likelihood Estimation on the
##first 20 weekly observations of the Schwartz and Smith (2000) Crude Oil Data:
##Example 1 - Stitched Futures Data:
##Setting a seed before parameter estimation is recommended:
set.seed(1)
Schwartz.Smith.Two.Factor.Stitched <- CPM.MLE(
 ####Arguments
 log.futures = log(SS.Oil$Stitched.Futures)[1:5,1],
 dt = SS.Oil$dt,
 TTM= SS.Oil$Stitched.TTM[1],
 S.Constant = FALSE, N.factors = 2, GBM = TRUE,
 cluster = NULL,
 hessian = TRUE,
 ####Genoud arguments:
 Richardsons.Extrapolation = FALSE,
 pop.size = 4, optim.method = "L-BFGS-B", print.level = 0,
 max.generations = 0, solution.tolerance = 0.1)

##Example 2 - All Available Contracts:
Schwartz.Smith.Two.Factor.Contract <- CPM.MLE(
 ####Arguments
 log.futures = log(SS.Oil$Contracts)[1:20,1:5],
 dt = SS.Oil$dt,
 TTM= SS.Oil$Contract.Maturities[1:20,1:5],
 S.Constant = TRUE,
 N.factors = 2, GBM = TRUE,
 cluster = NULL,
 hessian = TRUE,
 ####Genoud arguments:
 Richardsons.Extrapolation = FALSE,
 pop.size = 4, optim.method = "L-BFGS-B", print.level = 0,
 max.generations = 0, solution.tolerance = 0.1)

##Not Run - Perform Full Maximum Likelihood Estimation of a two-factor
##crude oil model, using suggested Genoud arguments:
##This function call can take several minutes

# Schwartz.Smith.Two.Factor.Stitched <- CPM.MLE(
#       ####Arguments
#       log.futures = log(SS.Oil$Stitched.Futures),
#       dt = SS.Oil$dt,
#       TTM= SS.Oil$Stitched.TTM,
#       S.Constant = FALSE, N.factors = 2, GBM = TRUE,
#       cluster = NULL,
#       Richardsons.Extrapolation = TRUE,
#       ####Genoud arguments:
#       pop.size = 1500, optim.method = "L-BFGS-B", print.level = 1,
#       max.generations = 100, solution.tolerance = 0.1)

##Compare results to that presented in Schwartz and Smith (2000):
# print(Schwartz.Smith.Two.Factor.Stitched$Estimated.Parameters)
# print(SS.Oil$Two.Factor)