\documentclass[a4paper,11pt]{article}
\usepackage[authoryear,round]{natbib}
\bibliographystyle{plainnat}
\usepackage{url}
%\usepackage{a4wide}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage[
colorlinks=true,
linkcolor=black,
anchorcolor=black,
citecolor=black,
filecolor=black,
pdfstartview={Fit},
pdfstartpage={1}
]
{hyperref}
\newcommand{\RR}{\textsf{R}}
\newcommand{\vnorm}[1]{\Vert#1\Vert} % Notation: 2-Norm of a Vector
\newcommand{\condexp}[2]{
  \mathbb{\tilde{E}}\left[#1\vert\mathcal{F}_{#2}\right]}
\newcommand{\condexprf}[2]{
  \mathbb{\tilde{E}}\left[#1\vert\mathcal{F}_{#2}\right]}
% Conditional Expectation
\newcommand{\Rpackage}{\textit{schwartz97}}
\newcommand{\sigmaS}{\sigma_S}
\newcommand{\sigmaE}{\sigma_\epsilon}
% <---------------------------------------------------------------------->
\SweaveOpts{engine=R,eps=TRUE}
%%\VignetteIndexEntry{User Guide}
\title{Using the \Rpackage{} package}
\author{David L\"{u}thi, Philipp Erb \& Simon Otziger}%\\\email{otis@zhaw.ch}}
\date{\today}
\begin{document}
\SweaveOpts{concordance=TRUE}
\maketitle
\begin{abstract}
  The purpose of this document is to show how the \RR{} package
  \Rpackage{} can be used. This is done by numerous examples and
  intuitive explanations. 
\end{abstract}
<<echo=FALSE>>=
library(schwartz97)
options(width = 60)  
@ 
%\tableofcontents
% <============================================================>
\section{Introduction}\label{sec:intro}
% <============================================================>
The package \Rpackage{} provides a set of functions to work with the
two-factor model of \cite{Gibson1990}\footnote{Because the model was
  extended in \cite{Schwartz1997} and \cite{Miltersen1998} we call it
  the \emph{Schwartz two-factor model} hereafter.}. The two-factor
model describes the joint dynamics of the state variables \emph{spot
  price} and \emph{spot convenience yield} (later simply called
\emph{convenience yield}).

We believe that the value of this package primarily lies in the
parameter fitting routine \verb=fit.schwartz2f=. Once the parameters
of the two-factor model are estimated, a number of functions can be
used to, e.g., draw samples from the model, price derivatives as
European options, or filter a futures price series to get an estimate
of the underlying state variables. The document is organized as
follows: Section \ref{sec:model} briefly describes the model and
section \ref{sec:overview} gives an overview of the classes and
functions. Then, section \ref{sec:init} to
\ref{sec:parameter-estimation} give examples and a case study.
% <============================================================>
\section{The Schwartz Two-Factor Model}\label{sec:model}
% <============================================================>
The spot price of the commodity and the instantaneous convenience
yield are assumed to follow the joint stochastic process:
\begin{eqnarray}
  dS_t &=& (\mu - \delta_t) S_t dt + \sigmaS S_t
  dW_S   \label{eq:Pdynamics-S} \\
  d\delta_t &=& \kappa (\alpha - \delta_t) dt + \sigmaE
  dW_\epsilon, \label{eq:Pdynamics-delta}
\end{eqnarray}
with Brownian motions $W_S$ and $W_\epsilon$ under the objective
measure $\mathbb{P}$ and correlation $dW_SdW_\epsilon = \rho dt$.
Under the pricing measure $\mathbb{Q}$ the dynamics are
\begin{eqnarray}
  dS_t &=& (r - \delta_t) S_t dt + \sigmaS S_t
  d\widetilde{W}_S   \label{eq:Qdynamics-S} \\
  d\delta_t &=& [\kappa (\alpha - \delta_t) - \lambda] dt + \sigmaE
  d\widetilde{W}_\epsilon, \label{eq:Qdynamics-delta}
\end{eqnarray}
where the constant $\lambda$ denotes the market price of convenience
yield risk and $\widetilde{W}_S$ and $\widetilde{W}_\epsilon$ are
$\mathbb{Q}$-Brownian motions. It may be handy to introduce a new
mean-level for the convenience yield process under $\mathbb{Q}$
\begin{equation}
  \label{eq:Q-drift}
 \tilde{\alpha} = \alpha - \lambda / \kappa.
\end{equation}
The dynamics is then
\begin{equation}
d\delta_t = \kappa (\tilde{\alpha} - \delta_t) dt + \sigmaE
  d\widetilde{W}_\epsilon.
\end{equation}
For more information on the model and pricing formulas we refer to the
other package vignette \emph{Technical Document}, or to
\cite{Schwartz1997} and \cite{Hilliard1998}.
% <============================================================>
\section{Package Overview}\label{sec:overview}
% <============================================================>
This section gives an overview of the functions and classes contained
in the package \Rpackage. In addition, the object-oriented programming
approach followed in this package is explained.
% <------------------------------------------------------------>
\subsection{Functions}\label{sec:functions}
% <------------------------------------------------------------>
The core of the package \Rpackage{} is built by the following
functions\footnote{There are also a number of utility functions as
  \emph{coef}, \emph{mean}, \emph{vcov}, \emph{plot}, \emph{resid},
  and \emph{fitted}.}:
\begin{center}
  \begin{tabular}{l|l}
    R-function & Description \\
    \hline
    \verb=dstate= & Density of the bivariate state vector. \\
    \verb=pstate= & Distribution of the bivariate state vector. \\
    \verb=qstate= & Quantile of the bivariate state vector. \\
    \verb=rstate= & Sample from the state distribution at some future time. \\
    \verb=simstate= & Generate trajectories from the bivariate state vector. \\
    \verb=dfutures= & Density of the futures price. \\
    \verb=pfutures= & Distribution of the futures price. \\
    \verb=qfutures= & Quantile of the futures price. \\
    \verb=rfutures= & Sample from the futures price distribution. \\
    \verb=pricefutures= & Calculate the futures price. \\
    \verb=priceoption= & Calculate the price of European call or put options. \\
    \verb=filter.schwartz2f= & Filter a futures price series. \\
    \verb=fit.schwartz2f= & Fit the two-factor model to data. \\
  \end{tabular}
\end{center}
Except the function \verb=fit.schwartz2f= all the above
functions are set to generic and accept three different signatures
(see section \ref{sec:oop}). 
% <------------------------------------------------------------>
\subsection{Classes}\label{sec:classes}
% <------------------------------------------------------------>
The package \Rpackage{} provides the class \verb=schwartz2f=. This
class contains all parameters which are needed to define the dynamics
of the state variables \emph{spot price} and \emph{convenience yield}
under the objective measure $\mathbb{P}$. The class
\verb=schwartz2f= has the following slots:\\[1ex]
\begin{center}
  \begin{tabular}{llll}
    Slot name & Class & Symbol & Description\\
    \hline
    s0 & numeric & $s_0$ & Initial spot price.\\
    delta0 & numeric & $\delta_0$ & Initial convenience yield.\\
    mu & numeric & $\mu$ & Drift parameter of the spot price.\\
    sigmaS & numeric & $\sigma_S$ & Diffusion parameter of the spot price.\\
    kappaE & numeric & $\kappa$ & Speed of mean-reversion of the convenience yield.\\
    alpha & numeric & $\alpha$ & Mean-level of the convenience yield.\\
    sigmaE & numeric & $\sigma_\epsilon$ & Diffusion parameter of the convenience yield.\\
    rhoSE & numeric & $\rho$ & Correlation between the two Brownian motions.\\
    call & call & & The function call.\\
  \end{tabular}
\end{center}
The above set of parameters contains the symbols appearing in
(\ref{eq:Pdynamics-S}) and (\ref{eq:Pdynamics-delta}) as well as the
initial values $s_0$ and $\delta_0$. To create an object of class
\verb=schwartz2f= the constructor with the same name can be used (see
section \ref{sec:init}).

The function \verb=fit.schwartz2f=, which estimates parameters of the
two-factor model, returns an object of class
\verb=schwartz2f.fit=. This class inherits from the class
\verb=schwartz2f= and adds the following slots.\\
\begin{center}
  \begin{tabular}{llll}
    Slot name & Class & Symbol & Description\\
    \hline
    r & numeric & $r$ & Risk-free interest rate.\\
    alphaT & numeric & $\tilde{\alpha}$ & Mean-value of the convenience yield under $\mathbb{Q}$.\\
    lambda & numeric & $\lambda$ & Market price of convenience yield risk.\\
    deltat & numeric & & Time-increment of the transition equation.\\
    n.iter & numeric & & Number of iterations.\\
    llh & numeric & & Log-likelihood value.\\
    converged & logical & & States whether the fit converged or not.\\
    error.code & numeric & & An error code or $0$.\\
    error.message & character & & Contains the error message if any.\\
    fitted.params & logical & & States which  parameters were fitted.\\
    trace.pars & matrix & & Parameter evolution during the estimation.\\
    meas.sd & numeric & & Standard deviation of the measurement equation.\\
    % lambda & numeric & $\lambda$ & Market price of convenience yield
    % risk.\\
    % alphaT & numeric & $\tilde{\alpha}$ & Drift parameter
    % $\tilde{\alpha}$ of the convenience yield under $\mathbb{Q}$.\\
    % r & numeric & $r$ & Risk-free interest-rate.\\
  \end{tabular}
\end{center}
These slots together with the ones contained in the class
\verb=schwartz2f= fully determine the dynamics of the model under
both, the objective measure and the pricing measure. Notice that one
of the parameters \verb=lambda= and \verb=alphaT= is redundant
according to equation \ref{eq:Q-drift}.
% <------------------------------------------------------------>
\subsection{Object Orientation}\label{sec:oop}
% <------------------------------------------------------------>
As mentioned earlier most of the functions dealing with the state
variables and futures prices are set to generic. The idea is to leave
some freedom to the user, who can decide whether he wants to use an
object-oriented approach or provide a fairly large set of arguments
for each function-call.

Consider the function \verb=dfutures= for example. The function
headers for different signatures are:
\begin{verbatim}
## S4 method for signature 'ANY,ANY,ANY,numeric':
dfutures(x, time = 0.1, ttm = 1, s0 = 50, delta0 = 0,
         mu = 0.1, sigmaS = 0.3, kappa = 1, alpha = 0,
         sigmaE = 0.5, rho = 0.75, r = 0.05, lambda = 0,
         alphaT = NULL, measure = c("P", "Q"), ...)

## S4 method for signature 'ANY,ANY,ANY,schwartz2f':
dfutures(x, time = 0.1, ttm = 1, s0, r = 0.05,
         lambda = 0, alphaT = NULL, measure = c("P", "Q"), ...)

## S4 method for signature 'ANY,ANY,ANY,schwartz2f.fit':
dfutures(x, time = 0.1, ttm = 1, s0, measure = c("P", "Q"), ...)
\end{verbatim}
Without object-orientation (first header) the function has 15
arguments. Ten parameters are needed to describe the dynamics under
both measures.

If a \verb=schwartz2f.fit= object is provided for \verb=s0= the only
additional arguments required are \verb=x= (quantiles), \verb=time=
(time where the futures process is evaluated), and \verb=ttm= (time to
maturity of the futures contract).
% <============================================================>
% <============================================================>
\section{Object Initialization}\label{sec:init}
% <============================================================>
A \verb=schwartz2f= object with reasonable parameters is
constructed in the following code chunk.
<<>>=
s0 <- 100
delta0 <- 0
mu <- 0.1
sigmaS <- 0.2
kappa <- 1
alpha <- 0.1
sigmaE <- 0.3
rho <- 0.4

obj <- schwartz2f(s0 = s0, delta0 = delta0, alpha = alpha,
                       mu = mu, sigmaS = sigmaS, sigmaE = sigmaE,
                       rho = rho, kappa = kappa)
obj
@

Objects of class \verb=schwartz2f.fit= are constructed via the
function \verb=fit.schwartz2f= (see section
\ref{sec:parameter-estimation}).
% <============================================================>
\section{Working with the state variables}\label{sec:state}
% <============================================================>
As soon as a \verb=schwartz2f= object is initialized, it can be passed
to the functions \verb=dstate=, \verb=pstate=, \verb=qstate=,
\verb=rstate=, and \verb=simstate=. The distribution of the state
variables depend on the horizon. Once this point in time is defined
the above functions can be used like the standard \RR{} distribution
functions for, e.g., the normal distribution (\verb=dnorm=,
\verb=pnorm=, \verb=qnorm=, \verb=rnorm=).

In this example a sample of the spot price and the convenience yield
in five years is generated by the function \verb=rstate=. Then, the
probability that the spot price is below 150 and the convenience yield
is lower than 0 in five years is computed. The mean of the state
variables in one and ten years is calculated next. Finally,
trajectories of the state variables are plotted (see
fig. \ref{fig:state-trajectory}).

<<fig=FALSE>>=
time <- 5
sample.t <- rstate(n = 2000, time, obj)
pstate(c(0, -Inf), c(150, 0), time, obj)
mean(obj, time = c(1, 10))
plot(obj, n = 30, time = 5, dt = 1 / 52)
@

\begin{figure}[h!]\centering
<<fig=TRUE,echo=FALSE,width=9,height=7>>=
plot(obj, n = 50, time = 5, dt = 1 / 52)
@
\caption{Thirty trajectories of the state variables are plotted on a
  weekly interval and a five years horizon. The initial values of the
  state variables are 100 for the spot price ($s_0$) and 0 for the
  convenience yield ($\delta_0$). The spot price has a drift $\mu$ of
  10\% and a volatility $\sigma_S$ of 20\%.  The \emph{speed of
    mean-reversion} parameter $\kappa$ of the convenience yield
  process is 1, and the long-term mean ($\alpha$) is 10\%.  The
  volatility of the convenience yield $\sigma_\epsilon$ is 30\% and
  the correlation $\rho$ between the Brownian motions driving the
  state variables is 40\%.}\label{fig:state-trajectory}
\end{figure}
% <============================================================>
\section{Working with derivatives}\label{sec:derivatives}
% <============================================================>
In this example we calculate some futures prices and plot the dynamics
of the term structure (``forward curve''). In addition prices of
European options are computed.

We work through this section by looking at corn and assuming all the
parameters are known. The current price ($s_0$) of 1000 bushels of
corn is assumed to be 80 USD. The convenience yield ($\delta_0$) is
zero at the moment but it's long-term mean ($\alpha$) is 5\%. The
drift ($\mu$) of corn is 10\% and the volatility is 30\%. The speed of
mean-reversion of the convenience yield ($\kappa$) is 1.5 and its
volatility is 40\%. Correlation is assumed to be 60\%. The risk-free
rate is 3\% and the market price of convenience yield risk ($\lambda$)
is zero.

First the object is initialized. Next a trajectory is generated based
on weekly sampling over five years. Then futures prices are calculated
with time to maturities ranging from zero (which is the spot) to two
years. Finally, a call option which matures in one year written on a
futures contract with time to maturity of two years is priced. Forward
curves are shown in fig. \ref{fig:forward-curves}.
<<>>=
s0 <- 80
delta0 <- 0.05
mu <- 0.1
sigmaS <- 0.3
kappa <- 1.5
alpha <- 0.05
sigmaE <- 0.4
rho <- 0.6
lambda <- 0.04
r <- 0.03
set.seed(1)
obj <- schwartz2f(s0, delta0, mu, sigmaS, kappa, alpha, sigmaE, rho)
state.traj <- simstate(n = 52 * time, time, obj)
pricefutures(seq(0, 2, by = 0.4), obj, lambda = lambda, r = r)
priceoption(type = "call", time = 1, Time = 2, K = 85,
            obj, r = r, lambda = lambda)
@
\setkeys{Gin}{width=\textwidth}
\begin{figure}[h!]\centering
<<fig=TRUE,echo=FALSE,height=7,width=9>>=
ttm <- seq(0, 2, by = 0.4)
forward.curves <- sapply(ttm, pricefutures, s0 = state.traj[,1], delta0 = state.traj[,2],
                         sigmaS = sigmaS, kappa = kappa, alpha = alpha, sigmaE = sigmaE,
                         lambda = lambda, r = r)
ttm.mat <- matrix(ttm, ncol = 6, nrow = 52 * time, byrow = TRUE)
dates <- Sys.Date() + seq(0, length = time * 52, by = 7)
rownames(ttm.mat) <- as.character(dates)
par(oma = c(6, 0, 5, 0), mfrow = c(2, 1), mar = c(0, 5.1, 0, 2.1))
## Prepare a list which is accepted by 'futuresplot'
futuresplot(list(price = forward.curves, ttm = ttm.mat * 260), xaxt = "n")
plot(dates, state.traj[,2], main = "", type = "l", ylab = "Convenience yield",
     xlim = c(dates[1], rev(dates)[1] + max(ttm) * 260))
abline(h = 0)

@
\caption{Forward curves with time to maturity up to two years are
  plotted for the trajectory \emph{state.traj}. The closest to
  maturity contract is in fact the spot price because the time to
  maturity is zero.}\label{fig:forward-curves}
\end{figure}
% <============================================================>
\section{Contango, Backwardation, and Hump Shapes}\label{sec:shape}
% <============================================================>
Fig. \ref{fig:forward-curves} shows the ability of the Schwartz
two-factor model to generate contango and backwardation
situations. Mixed shapes (humps and ``inverse humps'') are possible
too. E.g. an upwards sloping forward curve at the short end which
points downwards at the long end.

Looking at the $\mathbb{Q}$-dynamics in equations
(\ref{eq:Qdynamics-S}) and (\ref{eq:Qdynamics-delta}) it is obvious
that, locally, the drift of the spot price is positive when $\delta_t
< r$. This corresponds to a (local) contango situation. However, the
long-term mean of $\delta_t$, $\tilde{\alpha}$, defines the shape at
the far end of the term structure. Four different shapes are
generated in the following example and plotted in
fig. \ref{fig:shape}:
\begin{description}
\item[Pure contango:] If $\delta_t < r$ and $\tilde{\alpha} < r$.
\item[Short end backwardation, long end contango:] If $\delta_t > r$
  and $\tilde{\alpha} < r$.
\item[Pure backwardation:] If $\delta_t > r$ and $\tilde{\alpha} > r$.
\item[Short end contango, long end backwardation:] If $\delta_t < r$
  and $\tilde{\alpha} > r$.
\end{description}
<<>>=
s0 <- 1
delta0 <- 0.0
sigmaS <- 0.3
kappa <- 1
sigmaE <- 0.4
rho <- 0.5
r <- 0.03
ttm <- 0:4
## Pure contango
pricefutures(ttm, s0 = s0, delta0 = 0, sigmaS = sigmaS,
             kappa = kappa, sigmaE = sigmaE, rho = rho,
             r = r, alphaT = 0)
## Backwardation and then contango
pricefutures(ttm, s0 = s0, delta0 = 2 * r, sigmaS = sigmaS,
             kappa = kappa, sigmaE = sigmaE, rho = rho,
             r = r, alphaT = 0)
## Pure backwardation
pricefutures(ttm, s0 = s0, delta0 = r, sigmaS = sigmaS,
             kappa = kappa, sigmaE = sigmaE, rho = rho,
             r = r, alphaT = 2 * r)
## Contango and then backwardation
pricefutures(ttm, s0 = s0, delta0 = -r, sigmaS = sigmaS,
             kappa = kappa, sigmaE = sigmaE, rho = rho,
             r = r, alphaT = 2 * r)
@
             
\begin{figure}\centering             
<<fig=TRUE,echo=FALSE,width=9,height=7>>=
s0 <- 1
delta0 <- 0.0
sigmaS <- 0.3
kappa <- 1
sigmaE <- 0.4
rho <- 0.5
r <- 0.03
ttm <- seq(0, 4, by = 0.05)
## Pure contango
plot(ttm, ttm, type = "n", ylim = c(0.9, 1.18), xlab = "Time to maturity [y]",
     ylab = "Futures price", main = "Term structure shapes")
lines(ttm, pricefutures(ttm, s0 = s0, delta0 = 0, sigmaS = sigmaS,
                        kappa = kappa, sigmaE = sigmaE, rho = rho,
                        lambda = 0, r = r, alpha = 0), col = "blue", lty = 2)
## Backwardation and then contango
lines(ttm,pricefutures(ttm, s0 = s0, delta0 = 2 * r, sigmaS = sigmaS,
                       kappa = kappa, sigmaE = sigmaE, rho = rho,
                       lambda = 0, r = r, alpha = 0), col = "blue")
## Pure backwardation
lines(ttm,pricefutures(ttm, s0 = s0, delta0 = r, sigmaS = sigmaS,
                       kappa = kappa, sigmaE = sigmaE, rho = rho,
                       lambda = 0, r = r, alpha = 2 * r), lty = 2)
## Contango and then backwardation
lines(ttm,pricefutures(ttm, s0 = s0, delta0 = -r, sigmaS = sigmaS,
                       kappa = kappa, sigmaE = sigmaE, rho = rho,
                       lambda = 0, r = r, alpha = 2 * r))
legend("topleft", c("delta0 = 0%, alpha = 0%",
                    "delta0 = 6%, alpha = 0%",
                    "delta0 = 3%, alpha = 6%",
                    "delta0 = -3%, alpha = 6%"),
       title = "Common parameters: s0 = 1, sigmaS = 30%, kappa = 1, sigmaE = 40%, rho = 50%, lambda = 0, r = 3%",
       lty = c(2, 1, 2, 1), col = c(rep("blue", 2), rep("black", 2)), lwd = 2, bty = "n")
@
\caption{Four term structure shapes which can be generated by the
  Schwartz two-factor model. For a given risk-free interest rate $r$,
  the shape is determined by th value of the convenience yield
  $\delta_t$ and by the mean level of the convenience yield
  $\tilde{\alpha}$ under the pricing measure $\mathbb{Q}$ (see
  (\ref{eq:Qdynamics-S}) and
  (\ref{eq:Qdynamics-delta})).}\label{fig:shape}
\end{figure}       
\clearpage
% <============================================================>
\section{Parameter Estimation}\label{sec:parameter-estimation}
% <============================================================>
As mentioned in section \ref{sec:intro}, we believe that the package's
most valuable piece of code is the function
\verb=fit.schwartz2f=. This function estimates the parameters involved
in equations (\ref{eq:Pdynamics-S}) - (\ref{eq:Qdynamics-delta})
including the initial values of the state variables $s_0$ and
$\delta_0$. Because log-futures prices linearly depend on a bivariate
Gaussian random vector (the log-spot price and the convenience yield),
it is most straightforward to use a linear state-space model.
Therefore, the estimation procedure is based on the Kalman filter as
proposed in \cite{Schwartz1997}.

The header of the function \verb=fit.schwartz2f= looks like

<<>>=
  args(fit.schwartz2f)
@

The data inputs are \verb=data= and \verb=ttm=. \verb=data= must be a
regularly spaced time-series matrix of futures prices and \verb=ttm= a
matrix giving the time-to-maturity.

The time-to-maturity matrix admits the following interpretation:
\verb=data[i,j]= denotes the futures price whose time to maturity was
\verb=ttm[i,j]= when it was observed. The unit is defined by
\verb=deltat= which is the time between observations \verb=data[i,j]=
and \verb=data[i+1,j]=.

The arguments from \verb=s0= to \verb=lambda= are initial values of
the parameters. 

\verb=meas.sd= gives (initial) values of the measurement error
standard deviations. 

\verb=opt.pars= states which parameters shall be estimated. Observe
that some parameters are held constant by default.

\verb=opt.meas.sd= specifies how measurement uncertainty is treated in
the fit: According to the model there should be a one-to-one
correspondance 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 \verb=meas.sd=. All
components of \verb=meas.sd= can be fitted. However, it might be
sufficient to fit only a scalar where the measurement uncertainty is
parametrized by \verb=scalar * meas.sd= (this is the default). In this
case define the vector \verb=meas.sd= and set \verb=opt.meas.sd= to
``scalar''.  \verb=meas.sd= can be set to a vector with each component
set to the same value, thereby giving each point on 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 \verb=meas.sd=. If
all components of \verb=meas.sd= shall be fitted choose ``all''. If
the measurement uncertainty is known set \verb=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). This is reflected in zero
off-diagonals of the measurement error covariance matrix.

Finally, the risk-free rate \verb=r= must be given.

% <------------------------------------------------------------>
\subsection{Statistical and Computational
  Considerations}\label{sec:computational-considerations}
% <------------------------------------------------------------>
Estimation of the Schwartz two-factor model parameters is
statistically fragile and computationally demanding. Multiple local
maxima of the likelihood may exist which can result in absurd
parameter estimates as, e.g., a yearly drift of 300\% and or a market
price of convenience yield risk of -200\%. Therefore, a reasonable
parameter estimation is most likely an iteration where several initial
values are used and different combinations of parameters are held
constant during estimation. Also, simulation studies showed that a
fairly large sample is required to get adequate estimates (e.g. 20000
daily observations, depending on the number of parameters which shall
be estimated). For this reason the default is to hold \verb=s0=,
\verb=delta0=, and \verb=lambda= constant.

Several utility functions as \verb=fitted=, \verb=resid=, \verb=plot=,
and \verb=coef= may help to investigate the quality of the fit (see
example below).

The fitting procedure generally requires a large number of iterations
to achieve a reasonable tolerance level. Each optimization iteration
involves the filtering of the data set by the Kalman
filter. Therefore, an efficient implementation of the Kalman filter is
key. 

% <------------------------------------------------------------>
\subsection{Example: Estimating Wheat Parameters}\label{sec:estim-example}
% <------------------------------------------------------------>
This section takes you through a ``real-world'' example of a Schwartz
two-factor parameter estimation.  There are daily observations of the
five closest to maturity wheat futures prices from Jan. 1995 to April
2010 (approx. 4000 observations). 

The default parameters of \verb=fit.schwartz= are used, i.e., all
parameters except the initial values of the state variables and the
market price of convenience yield risk \verb=lambda= are estimated.
The maximum number of iterations is limited to 300 to save (build and
check) time. Then the object is printed and the parameter evolution is
plotted.
<<<fig=FALSE>>=
data(futures)
wheat.fit <- fit.schwartz2f(futures$wheat$price, futures$wheat$ttm / 260,
                          deltat = 1 / 260, control = list(maxit = 300), silent = TRUE)
wheat.fit
plot(wheat.fit, type = "trace.pars")
@ 

\begin{figure}[h!]\centering
<<fig=TRUE,echo=FALSE,width=9,height=7>>=
plot(wheat.fit, type = "trace.pars")
@ 
\caption{This figure shows the parameter evolution of the first 300
  iterations of the unconstrained parameter estimation of wheat. While
  the relative tolerance gets below $10^{-4}$ after around 200
  iterations the parameter values still fluctuate
  substantially.}\label{fig:param-evolution-raw}
\end{figure}

Observe that the fit did not converge as the maximum of 300 iterations
is not sufficient to achieve to required tolerance (here optim's
default). Let's discuss the parameters: A $\mu$ of 25\% is probably
not far off and the spot price volatility $\sigma_S$ of 37\% seems to
be fine too. The speed of mean reversion of the convenience yield
process $\kappa$ is alarmingly close to zero which means that the
convenience yield can drift far away from it's mean. The mean level of
the convenience yield $\alpha$ of -12\% is not intuitive to say the
least. The correlation $\rho$ of 90\% is very high, making the
parameters even harder to interpret and the process dynamic less
intuitive. In addition, the parameter evolution shown in
fig. \ref{fig:param-evolution-raw} is concerning.

Two measures are introduced as an attempt to make parameters more
appealing: $\kappa$ is constrained to 1 and the measurement error
standard deviations are set proportional to the average traded
volumes. The average initial value of the measurement error standard
deviations is set to 1\%.

<<fig=FALSE>>=
vol.std <- colSums(futures$wheat$vol, na.rm = TRUE) / sum(futures$wheat$vol, na.rm = TRUE)
wheat.fit.constr <- fit.schwartz2f(futures$wheat$price, futures$wheat$ttm / 260,
                                   kappa = 1,
                                   opt.pars = c(s0 = FALSE, delta0 = FALSE, mu = TRUE, 
                                     sigmaS = TRUE, kappa = FALSE, alpha = TRUE, 
                                     sigmaE = TRUE, rho = TRUE, lambda = FALSE),
                                   meas.sd = 1 / vol.std  / sum(1 / vol.std) * length(vol.std) * 0.01,
                                   deltat = 1 / 260, control = list(maxit = 300), silent = TRUE)

wheat.fit.constr
plot(wheat.fit.constr, type = "trace.pars")
@
\begin{figure}[h!]\centering
<<fig=TRUE,echo=FALSE,width=9,height=7>>=
plot(wheat.fit.constr, type = "trace.pars")
@ 
\caption{This figure shows the parameter evolution of the first 300
  iterations of the constrained parameter estimation of wheat. The
  relative tolerance gets below $10^{-6}$ after 150 iterations and the
  parameter values get more and more stationary. ``mu'' is an
  exception: Between iteration number 250 and 300 it goes from
  negative territory to 20\% and back to
  10\%.}\label{fig:param-evolution}
\end{figure}
Parameters are more reasonable now: $\mu$, $\alpha$, and $\rho$ seem
to be fine at 10\%, 4.3\%, and 54\%, respectively, as well as the
parameter evolution presented in fig. \ref{fig:param-evolution}. Also,
as a quick check, simulated trajectories in
fig. \ref{fig:simulated-trajectories} look plausible. Real term
structures are compared to model generated term structures in
fig. \ref{fig:term-structures}.
\begin{figure}[h!]\centering
<<fig=TRUE,width=9,height=7,echo=FALSE>>=
plot(wheat.fit.constr, type = "sim", time = 5, n = 30)
@ 
\caption{Thirty trajectories based on the constrained parameter
  estimates of wheat on a five years horizon and weekly
  sampling.}\label{fig:simulated-trajectories}
\end{figure}


<<fig=FALSE>>=
wheat.2007 <- lapply(futures$wheat,
                     function(x)x[as.Date(rownames(x)) > "2007-01-01" & as.Date(rownames(x)) < "2008-07-01",])
par(mfrow = c(1, 2))
futuresplot(wheat.2007, type = "forward.curve")
plot(wheat.fit.constr, type = "forward.curve", data = wheat.2007$price,
     ttm = wheat.2007$ttm / 260)
##xx <- filter.schwartz2f(data = wheat.2007$price,
##                        ttm = wheat.2007$ttm / 260, wheat.fit.constr)
##plot(ts(xx$state, frequency = 260, start = 2007))
@ 

\begin{figure}[h!]\centering
<<fig=TRUE,echo=FALSE,width=12,height=9>>=
par(mfrow = c(1, 2))
futuresplot(wheat.2007, type = "forward.curve")
plot(wheat.fit.constr, type = "forward.curve", data = wheat.2007$price,
     ttm = wheat.2007$ttm / 260)
@ 
\caption{Real term structures (left panel) and model generated term
  structures (right panel) from Jan. 2007 to June 2008. The
  backwardation starting at mid of 2007 is not captured by the model
  initially, and underestimated later in 2007. The shape of the term
  structure at the peak starting in 2008 can not be produced by the
  Schwartz two-factor model. The model's prediction looks reasonable
  in Q1 2007 and Q2 2008.}\label{fig:term-structures}
\end{figure}
\clearpage
% <------------------------------------------------------------>
\subsubsection{Residual Analysis}\label{sec:resid}
% <------------------------------------------------------------>
Model validation is done here via graphical residual
analysis. ``Residuals'' refer to prediction errors of the Kalman
filter's measurement equation. According to the model residuals should
be serially independent Gaussian random variables.

Different types of residulas can be obtained by the generic
\verb=resid= function and the specific argument \verb=type=, which can
be ``filter'' (raw filter residuals), ``filter.std'' (standardized
filter residuals), and ``real'' (observed minus fitted prices).
Standardized residuals are of interest here, hence the argument is
``filter.std''. First, serial independence is checked and then
normality of residuals. Both assumptions are violated as shown in
fig. \ref{fig:resid-acf} and fig. \ref{fig:resid-normality}.

Before rejecting the two-factor model one should first increase the
maximum number of iterations from 300 to, e.g., 3000 and investigate
the parameters estimates again. Then, different settings for the
measurement error standard deviations (argument \verb=meas.sd=) could
be tried, e.g. ``all''. Beside that different parameters could be hold
constant.

<<fig=FALSE>>=
model.resid <- resid(wheat.fit.constr, data = futures$wheat$price, ttm = futures$wheat$ttm / 260,
                     type = "filter.std")
acf(model.resid, na.action = na.pass)

par(mfrow = c(3, 2))
invisible(apply(model.resid, 2, function(x)plot(density(na.omit(x)))))
@ 
%<---------------------------------->
\begin{figure}[h!]\centering
<<fig=TRUE,echo=FALSE,width=12,height=9>>=
acf(model.resid, na.action = na.pass)
@ 
\caption{Residual's auto- and crosscorrelation estimates. Recall that
  off-diagonals of the measurement error covariance matrix are not
  estimated, hence the crosscorrelation is not relevant. Residuals of
  the closest to maturity futures show insignificant
  autocorrelation. However, residuals of all other futures are heavily
  autocorrelated. As the measurement error standard deviations are set
  proportional to the average traded volumes of the wheat futures, the
  two closest to maturity futures clearly get highest
  weights.}\label{fig:resid-acf}
\end{figure}
%<---------------------------------->

%<---------------------------------->
\begin{figure}[h!]\centering
<<fig=TRUE,echo=FALSE,width=12,height=9>>=
par(mfrow = c(3, 2))
dummy <- apply(model.resid, 2, function(x)plot(density(na.omit(x))))
@
\caption{Non of the residual distributions look normal. The first
  distribution looks normally shaped but exhibits outliers on both
  tails. The second to fifth distributions look asymmetric and fat
  tailed.}\label{fig:resid-normality}
\end{figure}
\clearpage
% <------------------------------------------------------------>
\subsubsection{Confidence Intervals and Value-at-Rsik}\label{sec:CI}
% <------------------------------------------------------------>
Confidence intervals and values-at-risk are computed here regardless
of the misspecification of the model estimated above. In order to
estimate the most recent values of the spot price and the convenience
yield the function \verb=filter.schwartz2f= is called first. Then, the
5\% and 95\% quantiles are computed and plotted (see
fig. \ref{fig:CI}) for a one week horizon . Note that the 5\% quantile
is the 95\% value-at-risk.

<<fig=FALSE>>=
state <- filter.schwartz2f(data = futures$wheat$price, ttm = futures$wheat$ttm / 260, wheat.fit.constr)$state
coefs <- coef(wheat.fit.constr)
n <- nrow(futures$wheat$price)
q.fut <- sapply(futures$wheat$ttm[n,] / 260, function(ttm, ...)
                qfutures(ttm = ttm, ...),
                p = c(0.05, 0.95), time = 5 / 260, s0 = state[n,1], delta0 = state[n,2],
                mu = coefs$mu, sigmaS = coefs$sigmaS, kappa = coefs$kappa, alpha = coefs$alpha,
                sigmaE = coefs$sigmaE, rho = coefs$rho, r = coefs$r, lambda = coefs$lambda)

plot(futures$wheat$ttm[n,], futures$wheat$price[n,], ylim = c(650, 850), type = "b",
     xlab = "Time to maturity [d]", ylab = "Price")
points(futures$wheat$ttm[n,], q.fut[1,], col = "blue", type = "b")
points(futures$wheat$ttm[n,], q.fut[2,], col = "blue", type = "b")
legend("topleft", c("Current observed futures price", "One week ahead 90% confidence interval"),
       fill = c("black", "blue"))
## q.fut <- sapply(1:nrow(state), function(i, s0, delta0, ttm, ...)
##                 qfutures(s0 = s0[i], delta0 = delta0[i], ttm = ttm[i], ...),
##                 p = 0.05, time = 1 / 12, ttm = futures$wheat$ttm / 260, s0 = state[,1],
##                 delta0 = state[,2], mu = coefs$mu, sigmaS = coefs$sigmaS, kappa = coefs$kappa, alpha = coefs$alpha,
##                 sigmaE = coefs$sigmaE, rho = coefs$rho, r = coefs$r, lambda = coefs$lambda)

## par(mfrow = c(1, 2))
## plot(ts(cbind(state[,1], futures$wheat$price[,1], q.fut), freq = 260, start = 1995),
##      plot.type = "single", col = c("black", "blue", "red"))
## plot(ts(state[,2], freq = 260, start = 1995))
## abline(h = coefs$alphaT)
@

\begin{figure}[t]\centering
<<fig=TRUE,echo=FALSE,width=9,height=6>>=
plot(futures$wheat$ttm[n,], futures$wheat$price[n,], ylim = c(650, 850), type = "b",
     xlab = "Time to maturity [d]", ylab = "Price")
points(futures$wheat$ttm[n,], q.fut[1,], col = "blue", type = "b")
points(futures$wheat$ttm[n,], q.fut[2,], col = "blue", type = "b")
legend("topleft", c("Current observed futures price", "One week ahead 90% confidence interval"),
       col = c("black", "blue"), lty = 1)
@
\caption{Current observed futures prices and one week ahead confidence
  intervals for the Sept. 2010, Dec. 2010, Mar. 2011, May 2011, and
  July 2011 wheat constracts as of Sept. 7, 2010. The July 2011
  futures is already at the lower bound of the confidence
  interval. This is due to the fact that the current shape of the term
  structure cannot be captured fully by the Schwartz two-factor
  model}\label{fig:CI}
\end{figure}
% <---------------------------------->
\clearpage
% <------------------------------------------------------------>
\subsection{Example: Estimating Soybean Meal Parameters}\label{sec:estim-example-soybean}
% <------------------------------------------------------------>
Soybean meal parameters are estimated in this example based on weekly
observations. All measurement error standard deviations are estimated
(\verb=opt.meas.sd == "all"=), but kappa and lambda are held constant.
Time-to-maturity (\verb=ttm=) is divided by 260 as it is in unit of
days and \verb=deltat= is set to 1/52 because weekly price
observations are used here.

First, data has to be made weekly. Wednesday observations are
taken. Second, the estimation is carried out. Finally, real and fitted
term structures are plotted in fig. \ref{fig:term-structure-soybean}.
Residuals could be analysed as outlined in sec. \ref{sec:resid}.
<<fig=FALSE>>=
futures.w <- rapply(futures, function(x)x[format(as.Date(rownames(x)), "%w") == 3,],
                    classes = "matrix", how = "list")
soybean.meal.fit <- fit.schwartz2f(data = futures.w$soybean.meal$price,
                                   ttm = futures.w$soybean.meal$ttm / 260,
                                   kappa = 1, mu = 0,
                                   opt.pars = c(s0 = FALSE, delta0 = FALSE, mu = TRUE, 
                                     sigmaS = TRUE, kappa = FALSE, alpha = TRUE, 
                                     sigmaE = TRUE, rho = TRUE, lambda = FALSE),
                                   opt.meas.sd = "all", deltat = 1 / 52,
                                   control = list(maxit = 1000), silent = TRUE)
soybean.meal.fit

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)
@ 

\begin{figure}[h!]\centering
<<fig=TRUE,echo=FALSE,width=9,height=7>>=
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)
@
\caption{Real and fitted term structures of soybean meal
  data. Directionally fitted term structures mostly match observed
  term structures.}\label{fig:term-structure-soybean}
\end{figure}

\clearpage
% <---------------------------------------------------------------------->
\bibliography{SpotCommodity}
\end{document}
