schwartz97 (version 0.0.6)

parameter-estimation: Schwartz 1997 two factor parameter estimation

Description

Fit the Schwartz 1997 two factor commodity model to futures data.

Usage

fit.schwartz2f(data, ttm, deltat = 1 / 260, s0 = data[1,1], delta0 = 0, mu = 0.1, sigmaS = 0.3, kappa = 1, alpha = 0, sigmaE = 0.3, rho = 0.7, lambda = 0, meas.sd = rep(0.1, ncol(data)), opt.pars = c(s0 = FALSE, delta0 = FALSE, mu = TRUE, sigmaS = TRUE, kappa = TRUE, alpha = TRUE, sigmaE = TRUE, rho = TRUE, lambda = FALSE), opt.meas.sd = c("scalar", "all", "none"), r = 0.03, silent = FALSE, ...)

Arguments

data

A matrix with futures prices. NA-values are allowed.

ttm

A matrix with the corresponding time to maturity (see Details).

deltat

Time increment (see Details).

s0

Initial value of the commodity spot price.

delta0

Initial value of the convenience yield.

mu

Initial value of the drift parameter of the commodity spot price.

kappa

Initial parameter of the speed of mean-reversion of the convenience yield process.

alpha

Initial parameter of the mean-level of the convenience yield process.

sigmaS

Initial parameter of the diffusion parameter of the spot price process.

sigmaE

Initial parameter of the diffusion parameter of the convenience yield process.

rho

Initial parameter of the correlation coefficient between the Brownian motion driving the spot price and the convenience yield process.

meas.sd

Initial parameter of the standard deviation of the measurement equation (see Details).

lambda

Initial value of the market price of convenience yield risk.

opt.pars

A logical vector giving the parameters which shall be fitted. The order is as given in the function header, names are discarded.

opt.meas.sd

States how the standard deviation of the measurement equation should be treated (see Details).

r

Instantaneous risk-free interest rate.

silent

If FALSE the log-likelihood and parameters will be printed in each iteration.

...

Arguments passed to optim.

Value

An object of class schwartz2f.fit.

Details

The elements of data and ttm have the following interpretation: data[i,j] denotes the futures price whose time to maturity was ttm[i,j] when it was observed (in units of deltat).

The time increment between observation data[i,j] and data[i + 1,j] is denoted with deltat. Note that this specification requires a regularly spaced data series.

opt.meas.sd specifies how measurement uncertainty is treated in the fit: According to the model there should be a one-to-one correspondence between the spot and the futures price. In reality, the term structure does not fully match for any set of parameters. This is reflected in the measurement uncertainty-vector meas.sd. All components of meas.sd can be fitted. However, it might be sufficient to fit only a scalar where the measurement uncertainty is parametrized by scalar * meas.sd. In this case define the vector meas.sd and set opt.meas.sd to “scalar”. meas.sd can be set to a vector with each component set to, e.g., 2%, giving each point in the term structure equal weight. Another reasonable specification takes open interest or volumes into account: The higher the volume, the higher the weight and therefore the lower the corresponding component of meas.sd. If all components of meas.sd shall be fitted choose “all”. If the measurement uncertainty is known set meas.sd to “none”. Note that the measurement errors are assumed to be independent in this implementation (even though the model and the filter do not require independence).

The model and its parameters are described in the Details section of the schwartz2f-class documentation and in the package vignette Technical Document.

References

The Stochastic Behavior of Commodity Prices: Implications for Valuation and Hedging by Eduardo S. Schwartz Journal of Finance 52, 1997, 923-973

See Also

schwartz2f.fit class, fitted, resid, plot, coef), pricefutures, futures-data.

Examples

# data(futures)
# 
# ## Estimate parameters for wheat data.
# ## (little precision required with reltol = 1e-3)
# fit.obj <- fit.schwartz2f(futures$wheat$price, futures$wheat$ttm / 260,
#                           deltat = 1 / 260, control = list(reltol = 1e-3))
# 
# ## See how parameter values evolved during the fit
# plot(fit.obj, type = "trace.pars")
# 
# ## Do the same but with lower tolerance level:
# high.precision.fit <- fit.schwartz2f(futures$wheat$price, futures$wheat$ttm / 260,
#                                      control = list(maxit = 3000, reltol = 5e-8))
# 
# high.precision.fit
# 
# plot(high.precision.fit, type = "trace.pars") # ...concerning parameter evolution.
# 
# ## Alpha becomes nonsensically high, kappa (speed of mean-reversion
# ## of the convenience yield) goes to zero. Solution: Choose different
# ## initial values and/or hold kappa constant at 1.
# 
# constrained.fit <- fit.schwartz2f(futures$wheat$price, futures$wheat$ttm / 260,
#                                   opt.pars = c(s0 = FALSE, delta0 = FALSE, mu = TRUE, 
#                                                sigmaS = TRUE, kappa = FALSE, alpha = TRUE, 
#                                                sigmaE = TRUE, rho = TRUE, lambda = TRUE),
#                                   alpha = 0, kappa = 1, 
#                                   control = list(maxit = 3000, reltol = 5e-8))
# 
# constrained.fit
# 
# plot(constrained.fit, type = "trace.pars")
# 
# ## The above parameters based on a fit - where kappa was held constant at 1 -
# ## look more reasonable.
# 
# ## These residuals should be iid standard normal
# model.resid <- resid(fit.obj, data = futures$wheat$price, ttm = futures$wheat$ttm / 260,
#                      type = "filter.std")
# acf(model.resid, na.action = na.pass)
# par(mfrow = c(3, 2))
# apply(model.resid, 2, function(x)plot(density(na.omit(x))))
# ## ...but are anything but iid standard normal.
# 
# ## ...though the fitted values look fine:
# fitted.futures <- fitted(fit.obj, futures$wheat$price, futures$wheat$ttm / 260)
# par(mfrow = c(1, 3))
# ### Plot futures prices
# plot(as.ts(futures$wheat$price), plot.type = "single", ylab = "Futures prices",
#      col = gray(seq(0.1, 0.9, length = 4)))
# ## Plot fitted values
# plot(as.ts(fitted.futures), plot.type = "single", ylab = "Fitted values",
#      col = gray(seq(0.1, 0.9, length = 4)))
# ## Plot relative difference
# plot(as.ts((fitted.futures - futures$wheat$price) / fitted.futures), plot.type = "single",
#      ylab = "Relative difference", col = gray(seq(0.1, 0.9, length = 4)))
# 
# 
# ## Try with kappa = 1, alpha = 0, and flexible standard deviationss of
# ## the measurement errors. Stop at 200 iterations.
# fit.obj.2 <- fit.schwartz2f(futures$wheat$price, futures$wheat$ttm / 260,
#                             opt.pars = c(s0 = FALSE, delta0 = FALSE, mu = TRUE, 
#                                          sigmaS = TRUE, kappa = FALSE, alpha = FALSE, 
#                                          sigmaE = TRUE, rho = TRUE, lambda = TRUE),
#                             alpha = 0, kappa = 1, opt.meas.sd = "all",
#                             deltat = 1 / 260, control = list(maxit = 200))
# 
# 
# model.resid.2 <- resid(fit.obj.2, data = futures$wheat$price, ttm = futures$wheat$ttm / 260,
#                        type = "filter.std")
# ## Residuals seem to be better:
# acf(model.resid.2, na.action = na.pass)
# par(mfrow = c(3, 2))
# apply(model.resid.2, 2, function(x)plot(density(na.omit(x))))
# 
# 
# ## The schwartz2f.fit-object 'fit.obj' can be used to do further calculations as
# ## pricing a put option on the wheat futures which matures in 1.1
# ## years. The option expires in 1 year. The strike price is the
# ## expected futures price in 1.1 years.
# priceoption("put", time = 1, Time = 1.1, K = pricefutures(1.1, fit.obj),
#             fit.obj)
# 
# 
# 
# ## Parameter estimation for weekly soybean meal data:
# ## Get Wednesday observations:
# futures.w <- rapply(futures, function(x)x[format(as.Date(rownames(x)), "%w") == 3,],
#                     classes = "matrix", how = "list")
# 
# ## Estimate soybean meal parameters (stop after 500 iterations).
# ## ttm (time-to-maturity) is divided by 260 as it is in unit of days and
# ## deltat = 1/52 because weekly price observations are used.
# ## Estimate all measurement error standard deviations (opt.meas.sd == "all"),
# ## but hold kappa, alpha, and lambda constant.
# soybean.meal.fit <- fit.schwartz2f(data = futures.w$soybean.meal$price,
#                                    ttm = futures.w$soybean.meal$ttm / 260,
#                                    opt.meas.sd = "all",
#                                    mu = 0, kappa = 1, alpha = 0.03, r = 0.04,
#                                    opt.pars = c(s0 = FALSE, delta0 = FALSE, mu = TRUE, 
#                                                 sigmaS = TRUE, kappa = FALSE, alpha = FALSE, 
#                                                 sigmaE = TRUE, rho = TRUE, lambda = FALSE),
#                                    deltat = 1 / 52, control = list(maxit = 500))
# 
# soybean.meal.fit 
# 
# plot(soybean.meal.fit, type = "trace.pars") # plot the parameter evolution
# 
# ## Plot real and predicted forward curves:
# par(mfrow = c(1, 2))
# futuresplot(futures.w$soybean.meal, type = "forward.curve")
# plot(soybean.meal.fit, type = "forward.curve", data = futures.w$soybean.meal$price,
#      ttm = futures.w$soybean.meal$ttm / 260)