Given a set of parameters of the N-factor model, filter term structure data using the Kalman filter.
NFCP.Kalman.filter(
parameter.values,
parameters,
log.futures,
dt,
TTM,
verbose = FALSE,
debugging = FALSE
)Vector of parameter values of an N-factor model. The NFCP.Kalman.filter function is designed for
application to optim type functions, and thus parameter values and
corresponding parameters different inputs within the function.
Vector of parameter names. Each element of parameters must correspond to its respective value
element in object parameter.values.
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.
Constant, discrete time step of observations
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.
logical. Should additional information be output? see values. When verbose = F, the NFCP.Kalman.filter function is significantly faster, see details
logical. Should additional filtering information be output? see values
NFCP.Kalman.filter returns a numeric object when verbose = F, which corresponds to the log-likelihood of observations.
When verbose = T, the NFCP.Kalman.filter function returns a list object of length seven with the following objects:
|
Log-Likelihood of observations |
|
vector. The final observation of the state vector |
|
matrix. All observations of the state vector, after the updating equation has been applied |
|
matrix. Estimated futures prices at each observation |
|
matrix. Estimation error of each futures contracts at each observation |
|
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 |
|
matrix. The theoretical and empirical volatility of futures returns for each observed contract as returned from the TSFit.Volatility function |
When debugging = T, 9 objects are returned in addition to those returned when verbose = T:
|
array. The covariance matrix at each observation point, with the third dimension indexing across time |
|
vector. The function of the Kalman filter covariance matrix at each observation point, with the third dimension indexing across time |
|
matrix. The Kalman Gain at each observation point, with the third dimension indexing across time |
|
matrix. d_td[t] (see details) |
|
matrix. Z_tz[t] (see details) |
|
matrix. G_tG[t] (see details) |
|
vector. C_tc[t] (see details) |
|
matrix. Q_tQ[t] (see details) |
|
matrix. HH (see details) |
NFCP.Kalman.filter applies the Kalman filter algorithm for observable log.futures prices against the input parameters of an N-factor model
provided through the parameter.values and parameters input vectors.
The NFCP.Kalman.filter function is
designed for subsequent input into optimization functions and is called within the N-factor parameter estimation function NFCP.MLE. The first input to the
NFCP.Kalman.filter function is a vector of parameters of an
N-factor model, with elements of this vector corresponding to the parameter names within the elements of input vector parameters.
When logical input verbose = F, the NFCP.Kalman.filter function calls the fkf_SP function of the FKF_SP package, which itself is a wrapper
of a routine of the Kalman filter written in C utilizing Sequential Processing for maximum computational efficiency (see fkf_SP for more details). When verbose = T,
the NFCP.Kalman.filter instead applies a Kalman filter algorithm written in base R and outputs several other list objects, including filtered values and
measures for model fit and robustness (see Returns)
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 NN state variables x_tx[t].
When GBM = TRUE:
log(S_t) = _i=1^N x_i,tlog(S[t]) = sum_i=1^n x[i,t]
When GBM = FALSE:
log(S_t) = E + _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:
dx_1,t = ^*dt + _1 dw_1tdx[1,t] = mu^* dt + sigma[1] dw[1]t
When GBM = FALSE:
dx_1,t = - (_1 + _1x_1,t)dt + _1 dw_1tdx[1,t] = - (lambda[1] + kappa[1] x[1,t]) dt + sigma[1] dw[t]t
And: dx_i,t =_i 1 - (_i + _ix_i,t)dt + _i dw_itdx[i,t] =_(i != 1) - (lambda[i] + kappa[i] x[i,t]dt + sigma[i] dw[i]t)
where: E(w_i)E(w_j) = _i,jE(w[i])E(w[j])
The following constant parameters are defined as:
param mu: long-term growth rate of the Brownian Motion process.
param EE: Constant equilibrium level.
param ^*=-_1mu^* = mu-lambda[1]: Long-term risk-neutral growth rate
param _ilambda[i]: Risk premium of state variable ii.
param _ikappa[i]: Reversion rate of state variable ii.
param _isigma[i]: Instantaneous volatility of state variable ii.
param _i,j [-1,1]rho[i,j] in [-1,1]: Instantaneous correlation between state variables ii and jj.
Measurement Error:
The Kalman filtering algorithm assumes a given measure of measurement error (ie. matrix HH). 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) and are assumed independent.
var s_is[i] Observation error of contract ii.
When S.Constant = T, the values of parameter s_is[i] are assumed constant and equal to parameter object 'sigma.contracts'. When S.Constant = F, observation errors are assumed unique, where the error of
futures contracts s_is[i] is equal to parameter object 'sigma.contract_' [i] (i.e. 'sigma.contract_1', 'sigma.contract_2', ...).
Kalman Filtering
The following section describes the Kalman filter equations used to filter the N-factor model.
The Kalman filter iteration is characterised by a transition and measurement equation. The transition equation develops the vector of state variables between discretised time steps (whilst considering a given level of covariance between state variables over time). The measurement equation relates the unobservable state vector to a vector of observable measurements (whilst also considering a given level of measurement error). The typical Kalman filter algorithm is a Gaussian process state space model.
Transition Equation: x_t|t-1 = c_t + G_t x_t-1 + Q_t _t hat(x)[t|t-1] = c[t] + G[t] * hat(x)[t-1] Measurement Equation: y_t = d_t + Z_t x_t|t-1 + H_t _that(y)[t] = d[t] + Z[t] * hat(x)[t|t-1]
t = 1, , n t = 1, ..., n
Where _teta[t] and _tepsilon[t] are IID N(0,I(m))N(0,I(m)) and iid N(0,I(d))N(0,I(d)) respectively.
The state vector follows a normal distribution, x_1 N(a_1, P_1)x[1] ~ N(a[1], P[1]), with a_1a[1] and P_1P[1] as the mean vector and variance matrix of the initial state vector x_1x[1], respectively.
When verbose = F, the NFCP.Kalman.filter function performs Kalman filtering through the fkf.SP function of the FKF.SP package for maximum filtering efficiency, which is
crucial when filtering and estimating parameters of term structure data under the N-factor model, which generally has many observations, contracts and unknown parameters. When verbose = T,
the NFCP.Kalman.filter function instead performs Kalman filtering in R, returning greater details of filtered objects (see Value)
The Kalman filter can be used for parameter estimation through the maximization of the Log-Likelihood value. See NFCP.MLE.
Filtering the N-factor model
let mm represent the number of observations at time tt
let nn represent the number of factors in the N-factor model
observable futures prices: y_t = [ln(F(t,T_1)), ln(F(t,T_2)), , ln(F(t,T_m))]'y[t] = [ln(F(t,T[1])), ln(F(t,T[2])), ..., ln(F(t,T[m]))]'
State vector: x_t=[x_1t,x_2t,,x_nt ]'x[t] = [x[1t], x[2t], ..., x[nt]]'
Measurement error: diag(H) = [s_1^2, s_2^2, , s_n^2]diag(H) = [s[1]^2, s[2]^2, ..., s[n]^2]
Where s_is[i] == "sigma.contract_" [i] when the logical of function NFCP.Parameters S.constant = F
When S.constant = T, s_1 = s_2 = = s_n = s[1] = s[2] = ... = s[n] = "sigma.contracts"
var ZZ is an m nm X n matrix, where each element [i,j][i,j] is equal to:
Z_i,j = e^-_i T_jZ[i,j] = e^(-kappa[i] * T[j])
var d_td[t] is an m 1m X 1 vector:
d_t=[A(T_1), A(T_2), , A(T_m)]'d[t]=[A(T[1]), A(T[2]), ..., A(T[m])]'
Under the assumption that Factor 1 follows a Brownian Motion, \(A(T)\) is given by: A(T) = ^*T-_i=1^N - 1-e^-_i T_i_i+12(_1^2T + _i.j 1 _i _j _i,j 1-e^-(_i+_j)T_i+_j)A(T) = mu^* * T - sum_i=1^N (1-e^(-kappa[i] T) lambda[i])/(kappa[i]) + 1/2 (sigma[1]^2 * T) + sum_i.j != 1 sigma[i] sigma[j] rho[i,j] (1 - e^(-(kappa[i] + kappa[j]) * T)) / (kappa[i] + kappa[j])
var v_tv[t] is a n 1n X 1 vector of serially uncorrelated Guassian disturbances with
E(V_t) = 0E(V[t]) = 0 and cov(v_t)=R^2cov(v[t])=R^2
Where:
diag(G_t) = [e^-_1 , e^-_2 , , e^-_n ]diag(G[t]) = [e^-kappa[1] tau, e^-kappa[2] tau, ..., e^-kappa[n] tau
Where =T-ttau = T - t
var w_tw[t] is an n 1n X 1 vector of serially uncorrelated Guassian disturbances where:
E(w_t) = 0E(w[t]) = 0 and cov(w_t) = Q_tcov(w[t]) = Q[t]
var c_t=[ t,0,,0]'c[t] = [mu * Delta t, 0, ..., 0]' is an N 1N X 1 vector of the intercept of the transition equation.
var Q_tQ[t] is equal to the covariance function, given by:
Cov_1,1(x_1,t,x_1,t) = _1^2tCov[1,1](x[1,t],x[1,t]) = sigma[1]^2 * t
Cov_i,j(x_i,t,x_j,t) = _i_j_i,j1-e^-(_i+_j)t_i+_jCov[i,j](x[i,t],x[j,t]) = sigma[i] sigma[j] rho[i,j] (1-e^-(kappa[i]+kappa[j]) * t) / (kappa[i] + kappa[j])
(see also cov_func)
Penalising poorly specified models
The Kalman filter returns non-real log-likelihood scores when the function of the covariance matrix becomes singular or its determinant becomes negative.
This occurs when a poorly specified parameter set is input. Non-real log-likelihood scores can break optimization algorithms. To circumvent this, the NFCP.Kalman.filter
returns a heavily penalized log-likelihood score whilst also returning a warning. Penalized log-likelihood scores are calculated by:
stats::runif(1, -1.5e6, -1e6)
Diffuse Kalman filtering
If the initial values of the state vector are not supplied within the parameters and parameter.values vectors (ie. Initial.State = F within the
NFCP.Parameters function), 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 QQ
Initial states of factors that follow an Ornstein-Uhlenbeck have a transient effect on future observations, however the initial value of a random walk variable persists across observations and therefore influencing model fit more (see Schwartz and Smith (2000) for more details).
Anderson, B. D. O. and J. B. Moore, (1979). Optimal filtering Englewood Cliffs: Prentice-Hall.
Fahrmeir, L. and G. tutz,(1994) Multivariate Statistical Modelling Based on Generalized Linear Models. Berlin: Springer.
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.
Durbin, J., and S. J. Koopman, (2012). Time series analysis by state space methods. Oxford university press.
# NOT RUN {
##Example 1 - complete, stitched data.
##Replicating the Schwartz and Smith (2000)
##Two-Factor commodity pricing model applied to crude oil:
Schwartz.Smith.Oil.Stitched <- NFCP.Kalman.filter(
parameter.values = SS.Oil$Two.Factor,
parameters = names(SS.Oil$Two.Factor),
log.futures = log(SS.Oil$Stitched.Futures),
TTM = SS.Oil$Stitched.TTM,
dt = SS.Oil$dt,
verbose = FALSE)
##Example 2 - incomplete, contract data.
##Replicating the Schwartz and Smith (2000)
##Two-Factor commodity pricing model applied to all available
##crude oil contracts:
SS.Oil.2F <- SS.Oil$Two.Factor
##omit stitched contract white noise
SS.Oil.2F <- SS.Oil.2F[!grepl("sigma.contract",
names(SS.Oil.2F))]
##Assume constant white noise in observable contracts:
SS.Oil.2F["sigma.contracts"] <- 0.01
#Kalman filter
Schwartz.Smith.Oil.Contracts <- NFCP.Kalman.filter(
parameter.values = SS.Oil.2F,
parameters = names(SS.Oil.2F),
## All available contracts are considered
log.futures = log(SS.Oil$Contracts),
## Respective 'TTM' of these contracts are input
TTM = SS.Oil$Contract.Maturities,
dt = SS.Oil$dt,
verbose = FALSE)
# }
Run the code above in your browser using DataLab