N-factor model parameter estimation through the Kalman filter and maximum likelihood estimation
the NFCP_MLE
function performs parameter estimation of an n-factor model given observable term structure futures data through maximum likelihood estimation.
NFCP_MLE
allows for missing observations as well as constant or variable time to maturity of observed futures contracts.
NFCP_MLE( log_futures, dt, futures_TTM, N_factors, N_ME = 1, ME_TTM = NULL, GBM = TRUE, estimate_initial_state = FALSE, Richardsons_extrapolation = TRUE, cluster = FALSE, Domains = NULL, ... )
log_futures |
Object of class |
dt |
Constant, discrete time step of observations |
futures_TTM |
Object of class |
N_factors |
|
N_ME |
|
ME_TTM |
vector of maturity groupings to consider for observed futures prices. The length of |
GBM |
|
estimate_initial_state |
|
Richardsons_extrapolation |
|
cluster |
an optional object of the 'cluster' class returned by one of the makeCluster commands in the |
Domains |
an optional |
... |
additional arguments to be passed into the |
NFCP_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 NFCP_Kalman_filter
function. When Richardsons_extrapolation = TRUE
, gradients are approximated
numerically within the optimization algorithm through the grad
function from the numDeriv
package. NFCP_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 are included within the vignette of NFCP
. All arguments of the genoud
function may be passed through the NFCP_MLE
function (except for gradient.check
, which is hard set to false).
NFCP_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 NFCP_domains
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 \(N\) state variables \(x_t\).
When GBM = TRUE
:
\[log(S_{t}) = \sum_{i=1}^N x_{i,t}\]
When GBM = FALSE
:
\[log(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
:
\[dx_{1,t} = \mu^*dt + \sigma_{1} dw_{1}t\]
When GBM = FALSE
:
\[dx_{1,t} = - (\lambda_{1} + \kappa_{1}x_{1,t})dt + \sigma_{1} dw_{1}t\]
And: \[dx_{i,t} =_{i\neq 1} - (\lambda_{i} + \kappa_{i}x_{i,t})dt + \sigma_{i} dw_{i}t\]
where: \[E(w_{i})E(w_{j}) = \rho_{i,j}\]
The following constant parameters are defined as:
param
\(\mu\): long-term growth rate of the Brownian Motion process.
param
\(E\): Constant equilibrium level.
param
\(\mu^*=\mu-\lambda_1\): Long-term risk-neutral growth rate
param
\(\lambda_{i}\): Risk premium of state variable \(i\).
param
\(\kappa_{i}\): Reversion rate of state variable \(i\).
param
\(\sigma_{i}\): Instantaneous volatility of state variable \(i\).
param
\(\rho_{i,j} \in [-1,1]\): Instantaneous correlation between state variables \(i\) and \(j\).
Disturbances - Measurement Error:
The Kalman filtering algorithm assumes a given measure of measurement error or disturbance in the measurement equation (ie. matrix \(H\)). 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). These disturbances are typically assumed independent.
var
\(ME_i\) measurement error of contract \(i\).
where the measurement error of futures contracts \(ME_i\) is equal to 'ME_'
[i] (i.e. 'ME_1'
, 'ME_2'
, ...) specified in arguments parameter_values
and parameter_names
.
There are three particular cases on how the measurement error of observations can be treated in the NFCP_Kalman_filter
function:
Case 1: Only one ME is specified. The Kalman filter assumes that the measurement error of observations are independent and identical.
Case 2: One ME is specified for every observed futures contract. The Kalman filter assumes that the measurement error of observations are independent and unique.
Case 3: A series of ME's are specified for a given grouping of maturities of futures contracts. The Kalman filter assumes that the measurement error of observations are independent and unique to their respective time-to-maturity.
Grouping of maturities for case 3 is specified through the ME_TTM
argument. This is a vector that specifies the maximum maturity to consider for each respective ME parameter argument.
in other words, ME_1 is considered for observations with TTM less than ME_TTM[1], ME_2 is considered for observations with TTM less than ME_TTM[2], ..., etc.
The first case is clearly the simplest to estimate, but can be a restrictive assumption. The second case is clearly the most difficult to estimate, but can be an infeasible assumption when considering all available futures contracts that make up the term structure of a commodity.
Case 3 thus serves to ease the restriction of case 1, and allow the user to make the modeling of measurement error as simple or complex as desired for a given set of maturities.
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 \(Q\)
Initial states of factors that follow an Ornstein-Uhlenbeck process 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).
NFCP_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 |
Filtered 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 |
Term Structure Volatility Fit |
matrix . The theoretical and empirical volatility of futures returns for each observed contract as returned from the TSFit.Volatility function |
proc_time |
list . The real and CPU time (in seconds) the NFCP_MLE function has taken. |
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/.
##Perform One Generation of Maximum Likelihood Estimation on the ##first 20 weekly observations of the Schwartz and Smith (2000) Crude Oil Data: SS_2F_estimated_model <- NFCP_MLE( ####Arguments log_futures = log(SS_oil$contracts)[1:20,1:5], dt = SS_oil$dt, futures_TTM= SS_oil$contract_maturities[1:20,1:5], N_ME = 1, N_factors = 1, GBM = TRUE, ####Genoud arguments: hessian = TRUE, Richardsons_extrapolation = FALSE, pop.size = 4, optim.method = "L-BFGS-B", print.level = 0, max.generations = 0, solution.tolerance = 10)
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.