2 Wilmott magazine
Alireza Javaheri, RBC Capital Markets
Delphine Lautier, Université Paris IX, Ecole Nationale Supérieure des Mines
de Paris
Alain Galli, Ecole Nationale Supérieure des Mines de Paris
Filtering in Finance
Further, we shall provide a mean to estimate the model parameters
via the maximization of the likelihood function.
1.1 The Simple and Extended Kalman Filters
1.1.1 Background and Notations
In this section we describe both the traditional Kalman Filter used for linear
systems and its extension to nonlinear systems known as the Extended
Kalman Filter (EKF). The latter is based upon a first order linearization of
the transition and measurement equations and therefore would coincide
with the traditional filter when the equations are linear. For a detailed
introduction, see Harvey (1989) or Welch and Bishop (2002).
Given a dynamic process xk following a transition equation
xk = f(xk-1, wk) (1)
we suppose we have a measurement zk such that
zk = h(xk, uk) (2)
where wk and uk are two mutually-uncorrelated sequences of temporallyuncorrelated
Normal random-variables with zero means and covariance
matrices Qk, Rk respectively4. Moreover, wk is uncorrelated with xk-1 and
uk uncorrelated with xk.
1 Filtering
The concept of filtering has long been used in Control Engineering and
Signal Processing. Filtering is an iterative process that enables us to estimate
a model's parameters when the latter relies upon a large quantity
of observable and unobservable data. The Kalman Filter is fast and easy
to implement, despite the length and noisiness of the input data.
We suppose we have a temporal time-series of observable data zk (e.g.
stock prices as in Javaheri (2002), Wells (1996), interest rates as in Babbs
and Nowman (1999), Pennacchi (1991), futures prices as in Lautier (2000),
Lautier and Galli (2000)) and a model using some unobservable timeseries
xk (e.g. volatility, correlation, convenience yield) where the index k
corresponds to the time-step. This will allow us to construct an algorithm
containing a transition equation linking two consecutive unobservable
states, and a measurement equation relating the observed data
to this hidden state.
The idea is to proceed in two steps: first we estimate the hidden state
a priori by using all the information prior to that time-step. Then using
this predicted value together with the new observation, we obtain a conditional
a posteriori estimation of the state.
In what follows we shall first tackle linear and nonlinear equations with
Gaussian noises. We then will extend this idea to the Non-Gaussian case.
Abstract
In this article we present an introduction to various Filtering algorithms and some of
their applications to the world of Quantitative Finance. We shall first mention the
fundamental case of Gaussian noises where we obtain the well-known Kalman Filter.
Because of common nonlinearities, we will be discussing the Extended Kalman Filter
(EKF) as well as the Unscented Kalman Filter (UKF) similar to Kushner's Nonlinear
Filter. We also tackle the subject of Non-Gaussian filters and describe the Particle
Filtering (PF) algorithm. Lastly, we will apply the filters to the term structure model of
commodity prices and the stochastic volatility model.
We denote the dimension of xk as nx, the dimension of wk as nw and
so on.
We define the a priori process estimate as
^xk
= E [xk] (3)
which is the estimation at time step k - 1 prior to the step k measurement.
Similarly, we define the a posteriori estimate
^xk = E [xk|zk] (4)
which is the estimation at time step k after the measurement.
We also have the corresponding estimation errors ek
= xk - ^xk
and
ek = xk - ^xk and the estimate error covariances
Pk
= E ek
e-t
k
Pk = E eket
k (5)
where the superscript t corresponds to the transpose operator.
In order to evaluate the above means and covariances we will need
the conditional densities p(xk|zk-1) and p(xk|zk), which are determined
iteratively via the Time Update and Measurement Update equations. The
basic idea is to find the probability density function corresponding to a
hidden state xk at time step k given all the observations z1:k up to that
time.
The Time Update step is based upon the Chapman-Kolmogorov equation
p(xk|z1:k-1) = p(xk|xk-1, z1:k-1)p(xk-1|z1:k-1)dxk-1
= p(xk|xk-1)p(xk-1|z1:k-1)dxk-1
(6)
via the Markov property.
The Measurement Update step is based upon the Bayes rule
p(xk|z1:k) =
p(zk|xk)p(xk|z1:k-1)
p(zk|z1:k-1)
(7)
with
p(zk|z1:k-1) = p(zk|xk)p(xk|z1:k-1)dxk
A proof of the above is given in the appendix.
EKF is based on the linearization of the transition and measurement
equations, and uses the conservation of Normal property within the class
of linear functions. We therefore define the Jacobian matrices of f with
respect to the state variable and the system noise as Ak and Wk respectively;
and similarly for h, as Hk and Uk respectively.
More accurately, for every row i and column j we have
Aij = fi/xj(^xk-1, 0), Wij = fi/wj(^xk-1, 0),
Hij = hi/xj(^xk
, 0), Uij = hi/uj(^xk
, 0)
^
TECHNICAL ARTICLE 1
Needless to say for a linear system, the function matrices are equal to
these Jacobians. This is the case for the simple Kalman Filter.
1.1.2 The Algorithm
The actual algorithm could be implemented as follows:
(1) Initialization of x0 and P0
For k in 1...N
(2) Time Update (Prediction) equations
^xk
= f(^xk-1, 0) (8)
and
Pk
= AkPk-1At
k + WkQ k-1Wt
k (9)
(3-a) Innovation: We define
^zk
= h(^xk
, 0)
and
k = zk - ^z-
k
as the innovation process.
(3-b) Measurement Update (Filtering) equations
^xk = ^xk
+ Kkk (10)
and
Pk = (I - KkHk)Pk
(11)
with
Kk = Pk
Ht
k(HkPk
Ht
k + UkRkUt
k)-1
(12)
and I the Identity matrix.
The above Kalman gain Kk corresponds to the mean of the conditional
distribution of xk upon the observation zk or equivalently, the matrix
that would minimize the mean square error Pk within the class of linear
estimators.
This interpretation is based upon the following observation. Having x
a Normally distributed random-variable with a mean mx and variance Pxx ,
and z another Normally distributed random-variable with a mean mz and
variance Pzz, and having Pzx = Pxz the covariance between x and z, the conditional
distribution of x|z is also Normal with
mx|z = mx + K(z - mz)
with
K = PxzP-1
zz
which corresponds to our Kalman gain.
1.1.3 Parameter Estimation
For a parameter-set in the model, the calibration could be carried out
via a Maximum Likelihood Estimator (MLE) or in case of conditionally
Gaussian state variables, a Quasi-Maximum Likelihood (QML) algorithm.
Wilmott magazine 3
4 Wilmott magazine
To do this we need to maximize N
k=1 p(zk|z1:k-1), and given the
Normal form of this probability density function, taking the logarithm,
changing the signs and ignoring constant terms, this would be equivalent
to minimizing over the parameter-set
L1:N =
N
k=1
ln(Pz k z k
) +
zk - mz k
Pz k z k
(13)
For the KF or EKF we have
mz k
= ^z-
k
and
Pz k z k
= HkPk
Ht
k + UkRkUt
k
1.2 The Unscented Kalman Filter and Kushner's
Nonlinear Filter
1.2.1 Background and Notations
Recently Julier and Uhlmann (1997) proposed a new extension of the
Kalman Filter to Nonlinear systems, different from the EKF. The new
method called the Unscented Kalman Filter (UKF) will calculate the
mean to a higher order of accuracy than the EKF, and the covariance to
the same order of accuracy.
Unlike the EKF, this method does not require any Jacobian calculation
since it does not approximate the nonlinear functions of the process and
the observation. Indeed it uses the true nonlinear models but approximates
the distribution of the state random variable xk (as well as the
observation zk ) with a Normal distribution by applying an Unscented
Transformation to it.
In order to be able to apply this Gaussian approximation (unless we
have xk = f(xk-1) + wk and zk = h(xk) + uk , i.e. unless the equations are
linear in noise) we will need to augment the state space by concatenating
the noises to it. This augmented state will have a dimension
na = nx + nw + nu .
1.2.2 The Algorithm
The UKF algorithm could be written in the following way:
(1-a) Initialization: Similarly to the EKF, we start with an initial choice
for the state vector ^x0 = E [x0] and its covariance matrix
P0 = E [(x0 - ^x0)(x0 - ^x0)t
].
We also define the weights W
(m)
i and W
(c)
i as
W
(m)
0 =

na + 
and
W
(c)
0 =

na + 
+ (1 - 2
+ )
and for i = 1...2na
W
(m)
i = W
(c)
i =
1
2(na + )
(14)
where the scaling parameters , ,  and  = 2
(na + ) - na will be chosen
for tuning.
For k in 1...N
(1-b) State Space Augmentation: As mentioned earlier, we concatenate
the state vector with the system noise and the observation noise, and create
an augmented state vector for each time-step
xa
k-1 =


xk-1
wk-1
uk-1


and therefore
^xa
k-1 = E [xa
k-1|zk] =


^xk-1
0
0


and
Pa
k-1 =


Pk-1 Pxw (k - 1|k - 1) 0
Pxw (k - 1|k - 1) Pww (k - 1|k - 1) 0
0 0 Puu (k - 1|k - 1)


(1-c) The Unscented Transformation: Following this, in order to use the
Normal approximation, we need to construct the corresponding Sigma
Points through the Unscented Transformation:
a
k-1(0) = ^xa
k-1
For i = 1...na
a
k-1(i) = ^xa
k-1 + ( (na + )Pa
k-1)i
and for i = na + 1...2na
a
k-1(i) = ^xa
k-1 - ( (na + )Pa
k-1)i-na
(15)
where the above subscripts i and i - na correspond to the ith
and i - nth
a
columns of the square-root matrix5.
(2) Time Update equations are
k|k-1(i) = f(x
k-1(i), w
k-1(i)) (16)
for i = 0...2na + 1 and
^xk
=
2na
i=0
W
(m)
i k|k-1(i) (17)
and
Pk
=
2na
i=0
W
(c)
i (k|k-1(i) - ^xk
)(k|k-1(i) - ^xk
)t
(18)
The superscripts x and w correspond to the state and system-noise portions
of the augmented state respectively.
(3-a) Innovation: We define
zk|k-1(i) = h(k|k-1(i), u
k-1(i)) (19)
and
^zk
=
2na
i=0
W
(m)
i Zk|k-1(i) (20)
^
Wilmott magazine 5
TECHNICAL ARTICLE 1
and as before
k = zk - ^z-
k
(3-b) Measurement Update
Pz k z k
=
2na
i=0
W
(c)
i (Zk|k-1(i) - ^zk
)(Zk|k-1(i) - ^zk
)t
and
Px k z k
=
2na
i=0
W
(c)
i (k|k-1(i) - ^xk
)(Zk|k-1(i) - ^zk
)t
(21)
which gives us the Kalman Gain
Kk = Px k z k
P-1
z k z k
and we have as previously
^xk = ^xk
+ Kkk (23)
as well as
Pk = Pk
- KkPz k z k
Kt
k (23)
which completes the Measurement Update Equations.
1.2.3 Parameter Estimation
The maximization of the likelihood could be done exactly as in (13) taking
^zk
and Pz k z k
defined as above.
1.2.4 Analogy with Kushner's Nonlinear Filter
It would be interesting to compare this algorithm to Kushner's Nonlinear
Filter6 (NLF) based on an approximation of the conditional distribution
as explained in Kushner (1967), Kushner and Budhiraja (2000). In this
approach, the authors suggest using a Normal approximation to the densities
p(xk|zk-1) and p(xk|zk). They then use the fact that a Normal distribution
is entirely determined via its first two moments, which reduces
the calculations considerably.
They finally rewrite the moment calculation equations (3), (4) and (5)
using the above p(xk|zk-1) and p(xk|zk), after calculating these conditional
densities via the time and measurement update equations (6) and
(7). All integrals could be evaluated via Gaussian Quadratures7.
Note that when h(x, u) is strongly nonlinear, the Gauss Hermite integration
is not efficient for evaluating the moments of the measurement
update equation, since the term p(zk|xk) contains the exponent zk - h(xk).
The iterative methods based on the idea of importance sampling proposed
in Kushner (1967), Kushner and Budhiraja (2000) correct this problem at
the price of a strong increase in computation time. As suggested in Ito and
Xiong (2000), one way to avoid this integration would be to make the additional
hypothesis that xk, h(xk)|z1:k-1 is Gaussian.
When nx = 1 and  = 2, the numeric integration in the UKF will correspond
to a Gauss-Hermite Quadrature of order 3. However in the UKF
we can tune the filter and reduce the higher term errors via the previously
mentioned parameters  and .
As the Kushner paper indicates, having an N-dimensional Normal
random-variable X = N (m, P) with m and P the corresponding mean
and covariance, for a polynomial G of degree 2M - 1 we can write
E[G(X)] =
1
(2)
N
2 |P|
1
2 RN
G(y)exp[(y
- m)t
P-1
(y - m)
2
]dy
which is equal to
E[G(X)] =
M
i1 =1
...
M
iN =1
wi1
...wiN
G(m +

P)
where t
= ( i1
... iN
) is the vector of the Gauss-Hermite roots of
order M and wi1
...wiN
are the corresponding weights.
Note that even if both Kushner's NLF and UKF use Gaussian
Quadratures, UKF only uses 2N + 1 sigma points, while NLF needs MN
points for the computation of the integrals.
More accurately, for a Quadrature-order M and an N-dimensional (possibly
augmented) variable, the sigma-points are defined for j = 1...N and
ij=1...M as
a
k-1(i1, ..., iN ) = ^xa
k-1 + Pa
k-1(i1, ..., iN )
where this square-root corresponds to the Cholesky factorization.
Similarly to the UKF, we have the Time Update equations
k|k-1(i1, ..., iN ) = f x
k-1(i1, ..., iN ), w
k-1(i1, ..., iN )
but now
^xk
=
M
i1 =1
...
M
iN =1
wi1
...wiN
k|k-1(i1, ..., iN )
and
Pk
=
M
i1 =1
...
M
iN =1
wi1
...wiN
(k|k-1(i1, ..., iN ) - ^xk
)(k|k-1(i1, ..., iN ) - ^xk
)t
and similarly for the measurement update equations.
1.3 The Non-Gaussian Case: The Particle Filter
It is possible to generalize the algorithm for the fundamental Gaussian
case to one applicable to any distribution.
1.3.1 Background and Notations
In this approach, we use Markov-Chain Monte-Carlo simulations instead
of using a Gaussian approximation for (xk|zk) as the Kalman or Kushner
Filters do. A detailed description is given in Doucet, De Freitas and
Gordon (2001).
The idea is based on the Importance Sampling technique:
We can calculate an expected value
E[ f(xk)] = f(xk)p(xk|z1:k)dxk (24)
by using a known and simple proposal distribution q().
6 Wilmott magazine
More precisely, it is possible to write
E[f(xk)] = f(xk)
p(xk|z1:k)
q(xk|z1:k)
q(xk|z1:k)dxk
which could be also written as
E[f(xk)] = f(xk)
wk(xk)
p(z1:k)
q(xk|z1:k)dxk (25)
with
wk(xk) =
p(z1:k|xk)p(xk)
q(xk|z1:k)
defined as the filtering non-normalized weight as step k.
We therefore have
E[f(xk)] =
Eq[wk(xk)f(xk)]
Eq[wk(xk)]
= Eq[ ~wk(xk)f(xk)] (26)
with
~wk(xk) =
wk(xk)
Eq[wk(xk)]
defined as the filtering normalized weight as step k.
Using Monte-Carlo sampling from the distribution q(xk|z1:k) we can
write in the discrete framework:
E[f(xk)] 
Nsims
i=1
~wk(x
(i)
k )f(x
(i)
k ) (27)
with again
~wk(x
(i)
k ) =
wk(x
(i)
k )
Nsims
j=1 wk(x
(j)
k )
Now supposing that our proposal distribution q() satisfies the Markov
property, it can be shown that wk verifies the recursive identity
w
(i)
k = w
(i)
k-1
p(zk|x
(i)
k )p(x
(i)
k |x
(i)
k-1)
q(x
(i)
k |x
(i)
k-1, z1:k)
(28)
which completes the Sequential Importance Sampling algorithm. It is important
to note that this means that the state xk cannot depend on future
observations, i.e. we are dealing with Filtering and not Smoothing8.
One major issue with this algorithm is that the variance of the
weights increases randomly over time. In order to solve this problem, we
could use a Resampling algorithm which would map our unequally
weighted xk's to a new set of equally weighted sample points. Different
methods have been suggested for this9.
Needless to say, the choice of the proposal distribution is crucial.
Many suggest using
q(xk|xk-1, z1:k) = p(xk|xk-1)
since it will give us a simple weight identity
w
(i)
k = w
(i)
k-1p(zk|x
(i)
k )
However this choice of the proposal distribution does not take into
account our most recent observation zk at all and therefore could
become inefficient.
Hence the idea of using a Gaussian Approximation for the proposal,
and in particular an approximation based on the Kalman Filter, in order
to incorporate the observations.
We therefore will have
q(xk|xk-1, z1:k) = N (^xk, Pk) (29)
using the same notations as in the section on the Kalman Filter. Such filters
are sometimes referred to as the Extended Particle Filter (EPF) and
the Unscented Particle Filter (UPF). See Haykin (2001) for a detailed
description of these algorithms.
1.3.2 The Algorithm
Given the above framework, the algorithm for an Extended or Unscented
Particle Filter could be implemented in the following way:
(1) For time step k = 0 choose x0 and P0 > 0.
For i such that 1  i  Nsims take
x
(i)
0 = x0 + P0Z(i)
where Z(i)
is a standard Gaussian simulated number.
Also take P
(i)
0 = P0 and w
(i)
0 = 1/Nsims
While 1  k  N
(2) For each simulation-index i
^x
(i)
k = KF(x
(i)
k-1)
with P
(i)
k the associated a posteriori error covariance matrix.
(KF could be either the EKF or the UKF)
(3) For each i between 1 and Nsims
~x
(i)
k = ^x
(i)
k + P
(i)
k Z(i)
where again Z(i)
is a standard Gaussian simulated number.
(4) Calculate the associated weights for each i
w
(i)
k = w
(i)
k-1
p(zk|~x
(i)
k )p(~x
(i)
k |x
(i)
k-1)
q(~x
(i)
k |x
(i)
k-1, z1:k)
with q() the Normal density with mean ^x
(i)
k and variance P
(i)
k .
(5) Normalize the weights
~w
(i)
k =
w
(i)
k
Nsims
i=1 w
(i)
k
(6) Resample the points ~x
(i)
k and get x
(i)
k and reset w
(i)
k = ~w
(i)
k = 1/Nsims .
(7) Increment k, Go back to step (2) and Stop at the end of the While loop.
^
Wilmott magazine 7
1.3.3 Parameter Estimation
As in the previous section, in order to estimate the parameter-set we
can use an ML Estimator. However since the particle filter does not necessarily
assume Gaussian noise, the likelihood function to be maximized
has a more general form than the one used in previous sections.
Given the likelihood at step k
lk = p(zk|z1:k-1) = p(zk|xk)p(xk|z1:k-1)dxk
the total likelihood is the product of the lk's above and therefore the loglikelihood
to be maximized is
ln(L1:N ) =
N
k=1
ln(lk) (30)
Now lk could be written as
lk = p(zk|xk)
p(xk|z1:k-1)
q(xk|xk-1, z1:k)
q(xk|xk-1, z1:k)dxk
and given that by construction the ~x
(i)
k 's are distributed according to q(),
considering the resetting of w
(i)
k to a constant 1/Nsims during the resampling
step, we can approximate lk with
~lk =
Nsims
i=1
w
(i)
k
which will provide us with an interpretation of the likelihood as the total
weight.
2 Term Structure Models of Commodity
Prices
In this section, the Kalman filter is applied to a well known term structure
model of commodity prices developed by Schwartz (1997). We first
present the Schwartz model, and show how it can be transformed into a
state-spaced model for a simple filter as well as an extended filter. We
then explain how the iteration process can be initiated, and how it can
be stabilized. Lastly, we compare the performance obtained with the simple
and extended filters, and make a sensitivity analysis.
2.1 The Schwartz model
The Schwartz model supposes that two state variables, namely the spot
price S and the convenience yield C, can explain the behavior of the
futures prices F. The convenience yield can briefly be defined as the comfort
associated with the possession of physical stocks. There are usually
no empirical data for these two variables, because most of the time there
are no reliable time series for the spot price, and the convenience yield is
not a traded asset.
2.1.1 Presentation of the model
The dynamics of these state variables are the following:
dS = ( - C)Sdt + SSdzS
dC = [( - C)]dt + CdzC
with:
E[dzS × dzC] = dt
, S, C > 0
where:
--  is the immediate expected return for the spot price S,
-- S is the spot price's volatility,
-- dzS is the increment of the Brownian motion associated with S,
--  is the long run mean of the convenience yield C,
--  represents the convergence of the convenience yield towards ,
-- C is the convenience yield's volatility,
-- dzC is the increment of the Brownian motion associated with C.
--  is the correlation between the two Brownian motions associated
with S and C,
The model's solution expresses the relationship between an observable
futures price F for a delivery in T, and the state variables. This solution
is:
F(S, C, t, T) = S(t) × exp -C(t)
1 - e-

+ B()
with:
B() = r -  +
2
C
22
-
SC

×  +
2
C
4
×
1 - e-2
3
+   + SC -
2
C

×
1 - e-
2
,
 =  - (/)
where:
-- r is the risk-free interest rate10,
--  is the risk premium associated with the convenience yield,
--  = T - t is the maturity of the futures contract.
The model risk-neutral parameter-set is therefore:  = (S, ,
, C, , )
2.1.2 Applying the simple filter to the Schwartz model
To apply the simple filter, the Schwartz model must be expressed under a
linear form:
ln(F(S, C, t, T)) = ln(S(t)) - C(t) ×
1 - e-

+ B()
Considering the relationship G = ln(S), we also have:
dG =  - C - 1
2
2
S dt + SdzS
dC = [k( - C)]dt + CdzC
Then, to apply the Kalman filter, the model must be expressed under its
state-spaced form. The transition equation is:
Gt
Ct
= c + F ×
Gt-1
Ct-1
+ Twt, t = 1, . . . NT
TECHNICAL ARTICLE 1
8 Wilmott magazine
where:
-- c =
 - 1
2
2
S t
 t
is a (2 × 1) vector, and t is the period separating
2 observation dates
-- F =
1 - t
0 1 -  t
is a (2 × 2) matrix,
-- T is an identity matrix, (2 × 2),
-- wt is a sequence of uncorrelated errors, with:
E[wt] = 0, and Q = Var[wt] =
2
S t SC t
SC t 2
C t
The measurement equation is directly written from the solution of
the model:
zt = d + H ×
Gt
Ct
+ ut, t = 1, . . . NT
where:
-- The ith
component of the vector zt of size N is ln(F(i)). N is the
number of maturities that were used for the estimation.
-- The ith
component of the vector d of size N is B(i)
-- H is the (N × 2) matrix whose ith
line is 1, -
1 - e-i

-- ut is a white noise vector, of size N, with no serial correlation:
E[ut] = 0 and R = Var[ut]. R is (N × N)
2.1.3 Applying the extended filter to the Schwartz model
The transition equation is directly obtained from the dynamics of the
state variables:
St
Ct
= f(St-1, Ct-1) + T(St-1, Ct-1)wt
where:
--
St
Ct
is the state vector, (2 × 1),
-- f(St-1, Ct-1) is a (2 × 1) vector:
f(St-1, Ct-1) =
St-1(1 +  t - Ct-1 t)
 t + Ct-1(1 -  t)
--T(St-1, Ct-1) is a (2 × 2) matrix: T(St-1, Ct-1) =
St-1 0
0 1
-- Q is a (2 × 2) matrix: Var(wt) =
2
S SC
SC 2
C
The measurement equation becomes:
zt = h(St, Ct) + ut
where:
-- The ith
component of the vector zt is F(i)
-- h(St, Ct) is a vector whose ith
component is:
[St × exp(-ziCt/t-1 + B(i))]
with: Zi =
1 - e-i

B(i) = r -  +
2
C
22
-
SC

× i +
2
C
4
×
1 - e-2i
3
+  + SC -
2
C

×
1 - e-i
2
 =  - /
-- ut is a white noise vector, (N × 1), with no serial correlation:
E[ut] = 0 and R = Var[ut]. R is (N × N)
Lastly A and H, the derivatives of the functions f and h with respect to the
state variables, are the following:
-- A is a (2 × 2) matrix:
A(St-1, Ct-1) =
1 +  t - Ct-1 t -St-1 t
0 (1 -  t)
-- H is a (N × 2) matrix whose ith
line is:
e(-Zi Ct-1 +B(i ))
- Zi × e(-zi Ct-1 +B(i ))
B(i) = r -  +
2
C
22
-
SC

× i +
2
C
4
×
1 - e-2i
3
+  + SC -
2
C

×
1 - e-i
2
 =  - /
2.2 Practical difficulties associated with the empirical
study
To perform the empirical study, some difficulties must be overcome.
First, there are choices to make when the iterative process is started.
Second, if the model has been expressed under the logarithmic form for
the simple Kalman filter, some precautions must be taken when the performance
is being considered. Third, the stability of the iteration process
and the model's performance are extremely sensitive to the covariance
matrix R.
2.2.1 Starting the iterative process
To start the iterative process, there is a need for the initial values of the
non-observable variables and for their covariance matrix. In the case of
^
Wilmott magazine 9
term structure models of commodity prices, the nearest futures price is
generally used as the spot price S, and the convenience yield C is calculated
as the solution of the Brennan and Schwartz (1985) model. This
solution requires the use of two observed futures prices, for delivery in
T1 and in T2:
c = r ln(F(S,
t, T1)) - ln(F(S, t, T2))
T1 - T2
where T1 is the nearest delivery, and T2 is the one immediately afterwards.
We choose the first value of the estimation period for the nonobservable
variables.
The covariance matrix associated with the state variables must also be
initialized. We choose a diagonal matrix, and calculate the variances
with the first 30 points of the estimation period.
2.2.2 Measuring the performance
As we have transformed the equations taking the logarithms, in order to
use the simple Kalman filter, care has to be taken not to introduce a bias
when transforming back the results. Indeed in that case, the innovations
are calculated with the logarithms of the futures prices.
zt = ^zt
+ V
where  is the standard deviation of the innovations and V is a standardized
Gaussian residual independent to ^zt
. The exponential of ^zt
is a
biased estimator of the future prices, because ezt
= e^zt
× e V
and taking
the expectation we find:
E[ez
t ] = E[e^zt
] × e
 2
2
We therefore have to correct it by the exponential factor.
2.2.3 Stabilizing the iteration process
Another important choice must be made before initiating the Kalman
iteration process, concerning the estimation of the covariance matrix
associated with the errors introduced in the measurement equation. This
matrix R is crucial for the iteration's stability, because it is added to the
innovations covariance matrix, during the innovation phase. Therefore,
the updating of the non-observable variable is strongly affected by the
matrix R, and if the terms of this matrix are too high, the iteration can
become unstable.
Most of the time, the easiest way to estimate this matrix is to calculate
the variances and the covariances of the estimation database. This
method was chosen to measure the model's performance and is presented
in paragraph 2.3. But it is important to know how much the empirical
results are affected by this choice. In order to answer this question, a few
simulations will be run and analyzed.
2.3 Comparison between the two filters
The comparison between the performance of the Schwartz model measured
with the two filters allows us to appreciate the influence of the
linearization on the results. We first present the empirical data. Then
the performance criteria are exposed. Finally, the results are delivered
and commented upon.
2.3.1 The empirical data
The data used for the empirical study corresponds to daily crude oil
prices for the settlement of the Nymex's WTI futures contracts,
between the 18th
of May 1998, and the 15th
of October 2001. They have
been arranged such that the first futures price's maturity 1 is actually
the one month maturity, and the second futures price's corresponds
to the two months maturity 2 ,. . . Keeping the first observation
of each group of five, these daily data points were transformed
into weekly data. For the parameters estimation, and for the measure
of the model's performance, four series of futures prices11 were
retained, corresponding to the one, the three, the six and the nine
months maturities. The interest rates are set to the three-month T-bill
rates. Because interest rates are supposed to be constant in the model,
we use the mean of all the observations between 1995 and 2002.
2.3.2 The performance criteria
To measure the model's performance, two criteria were retained: the
mean pricing errors and the root mean squared errors.
The mean pricing errors (MPE) are defined in the following way:
MPE =
1
N
N
n=1
(^F(n, ) - F(n, ))
where N is the number of observations, ^F(n, ) is the estimated futures price
for a maturity  at time-step n, and F(n, ) is the observed futures price.
Retaining the same notations, the root mean squared error (RMSE), is
defined in the following way, for one given maturity :
RMSE =
1
N
N
n=1
(^F(n, ) - F(n, ))2
2.3.3 The empirical results
First, the optimal parameters obtained with the two filters are compared.
Then the model's ability to represent the price curve and its dynamics is
appreciated. Finally, the sensitivity of the results to the errors covariance
matrix is presented.
2.3.3.1 Optimal parameters12
The differences between the optimal parameters obtained with the two
filters show that the linearization has a significant influence.
Nevertheless, the parameters have the same order size as those Schwartz
obtained in 1997 on the crude oil market, on different periods.
2.3.3.2 The model's performance
The simple filter is always more precise than the extended one. This is true
for all the maturities14. These measures also always decrease with maturity,
TECHNICAL ARTICLE 1
10 Wilmott magazine
which is consistent with Schwartz's results on other periods. Nevertheless,
Schwartz has worked with longer maturities, and shown that the root
mean squared error increases again for deliveries after 15 months.
To be rigorous, the model's performance associated with the simple
Kalman filter should be corrected when, as is the case here, the logarithms
of the estimations are used to obtain the estimations themselves.
The correction slightly improves the performance, as shown in table 3:
the root mean squared errors and the mean pricing errors diminish a
bit for almost all the maturities.
Finally, the innovations range diminishes with the futures contracts
maturity. Figure 1 illustrates the innovations behavior for the
one-month maturity case. It shows that they tend to return to zero for
the two filters. Nevertheless as the figure illustrates it, even if the
mean pricing errors are low for the two filters, the pricing errors can
be rather important at certain specific dates. They reach USD 6 in
absolute value for the extended filter, which represents 25% of the
mean futures price for the one-month maturity case. It is a bit lower
for the simple filter: USD 5.
Consequently we can say that there is clearly an impact of the linearization
introduced in the extended filter: it can be observed from the
optimal parameters, from the performance, and from the innovations.
Nevertheless with an extended filter, the model's ability to represent the
price curve is still quite good for this case.
Another way to assess a model's performance is to see if it is able to
reproduce the price dynamics, which can be shown graphically.
Considering this, the first important conclusion is that the model is
able to reproduce the price dynamics quite precisely, even if as in 1998--
2001, there are very large fluctuations in the futures prices. Figure 2
shows the results obtained for the one-month maturity case. During
that period, the crude oil price goes from USD 11 per barrel to USD 37
per barrel! Even if the Kalman filters are often suspected to be unstable,
these results show that they can be used even with extremely volatile
data. The graph also shows that the two filters attenuate the range of
TABLE 1. OPTIMAL PARAMETERS, 1998­200113
Simple filter Extended filter
Parameters Gradients Parameters Gradients
Pull back force:  1.59171 -0.003631 1.258133 0.000628
Trend:  0.379926 0.000497 0.352014 -0.001178
Spot price's volatility: s 0.263525 -0.000448 0.320235 -0.000338
Long run mean:  0.252260 -0.012867 0.232547 0.004723
Convenience yield's volatility: c 0.237071 -0.000602 0.288427 -0.001070
Correlation coefficient:  0.938487 -0.001533 0.969985 0.000008
Risk premium:  0.177159 0.009272 0.181955 -0.002426
TABLE 2. THE MODEL'S PERFORMANCE
WITH THE SIMPLE AND THE EXTENDED
FILTERS, 1998­2001
Simple filter Extended filter
Maturity MPE RMSE MPE RMSE
1 month -0.060423 2.319730 0.09793 2.294503
3 months -0.107783 1.989428 0.057327 2.120727
6 months -0.054536 1.715223 0.109584 1.877654
9 months -0.007316 1.567467 0.141204 1.695222
Average -0.057514 1.897962 0.101511 1.997027
Unit: USD/b.
FIGURE 1: Innovations, 1998­2001
-7
-6
-5
-4
-3
-2
-1
0
1
2
3
4
5
6
7
18/05/9818/07/9818/09/9818/11/9818/01/9918/03/9918/05/9918/07/9918/09/9918/11/9918/01/0018/03/0018/05/0018/07/0018/09/0018/11/0018/01/0118/03/0118/05/0118/07/0118/09/01
($/b)
Simple filter Extended filter
TABLE 3. THE COMPARISON BETWEEN
THE MODEL'S PERFORMANCE ASSOCIATED
WITH THE SIMPLE FILTER, WITH
AND WITHOUT CORRECTIONS FOR THE
LOGARITHM, 1998­2001
Simple filter Simple filter corrected
Maturity MPE RMSE MPE RMSE
1 month -0.060423 2.319730 0.065644 2.314178
3 months -0.107783 1.989428 0.006419 1.981453
6 months -0.054536 1.715223 0.026010 1.709931
9 months -0.007316 1.567467 0.061301 1.564854
Average -0.057514 1.897962 0.036637 1.892604
Unit: USD/b.
^
Wilmott magazine 11
Most of the time, the terms of this matrix correspond to the variances
and the covariances of the estimation database, namely, in the case studied
here, the variances and covariances between futures prices for different
maturities. But one must know that the results obtained with the
Kalman filter can be more precise if these terms are (artificially) lowered,
as shown in table 4. This table exposes the different results obtained during
1998--2001 with the extended Kalman filter. This period is especially
interesting because the data strongly fluctuates. The performance is
achieved, first with the matrix based on the observations, then with the
artificially lowered one.
Simulations 1 to 4 correspond respectively to the model's performance
obtained by multiplying the system's matrix R by (1/2), (1/16),
(1/160), and (1/1600). As the matrix elements are lowered, the model's
performance strongly improves: from the initial simulation to the fourth
one, the root-mean-squared-error is almost divided by two. The comparison
between the third and the fourth simulations also illustrates the fact
that there is a limit to the performance amelioration. Figure 3 portrays
the main results of these simulations.
2.4 Conclusion
The main conclusions of this section are the following: Firstly, the
approximation introduced in the extended Kalman filter has clearly an
influence on the model's performance; the extended filter generally
leads to less precise estimations than the simple one. Nevertheless, the
difference between the two filters is quite low and the extended filter is
still acceptable. Secondly, the estimation results are sensitive to the system's
matrix containing the errors of the measurement equation, and
this matrix can be used to obtain more precise results on the estimation
database. Thirdly, the approximation made in the extended Kalman filter
is not a real problem until the model starts to become highly nonlinear.
The following synthetic example illustrates the situation.
TECHNICAL ARTICLE 1
FIGURE 2: Estimated and observed futures prices for the one month maturity,
1998­2001
10
15
20
25
30
35
18/05/9818/07/9818/09/9818/11/9818/01/9918/03/9918/05/9918/07/9918/09/9918/11/9918/01/0018/03/0018/05/0018/07/0018/09/0018/11/0018/01/0118/03/0118/05/0118/07/0118/09/01
($/b)
Futures one month Simple filter Extended filter
TABLE 4. SIMULATIONS WITH DIFFERENT
SYSTEM'S MATRIX R
1 month 3 months 6 months 9 months Average
Observations
MPE 0.0979 0.0573 0.1096 0.1412 0.1015
RMSE 2.2945 2.1207 1.8777 1.6952 1.9970
Simulation 1
MPE 0.0013 0.0935 0.1501 1.6506 0.4739
RMSE 1.8356 1.5405 1.2478 2.6602 1.8210
Simulation 2
MPE 0.0073 0.0152 0.0612 0.0137 0.0244
RMSE 1.4759 1.1686 0.9386 0.8317 1.1037
Simulation 3
MPE 0.0035 -0.0003 0.0383 0.0005 0.0105
RMSE 1.3812 1.0950 0.8647 0.7499 1.0227
Simulation 4
MPE 0.0131 0.0067 0.0415 0.0075 0.0172
RMSE 1.3602 1.0919 0.8697 0.7591 1.0202
price fluctuations. This phenomenon can actually be observed for
every maturity.
2.3.3.3 Simulations
The last results presented in this sesction are simulations. They show
how the model's performance is affected by the choice of the system's
matrix R. This matrix represents the errors in the measurement equation
and the way it is estimated has a strong influence on the empirical
results.
FIGURE 3: One month futures prices observed/estimated, 1998­2001
10,00
15,00
20,00
25,00
30,00
35,00
18/05/9818/07/9818/09/9818/11/9818/01/9918/03/9918/05/9918/07/9918/09/9918/11/9918/01/0018/03/0018/05/0018/07/0018/09/0018/11/0018/01/0118/03/0118/05/0118/07/0118/09/01
($/b)
Observed futures prices for a one month's maturity
Estimated futures prices for a one month's maturity, with a matrix R corresponding to the observations
Estimated futures prices for a one month's maturity and an artificially lowered matrix (simulation 4)
12 Wilmott magazine
initial series and the estimations we obtain with the two versions of
the Kalman filter when a = 0.8.
3 Stochastic Volatility Models
In this section, we apply the different filters to a few stochastic volatility
models including the Heston, the GARCH and the 3/2 models. To
test the performance of each filter, we use five years of S&P500 time-
series.
The idea of applying the Kalman Filter to Stochastic Volatility models
goes back to Harvey, Ruiz & Shephard (1994), where the authors attempt
to determine the system parameters via a QML Estimator. This approach
has the obvious advantage of simplicity, however it does not account for
the nonlinearities and non-Gaussianities of the system.
More recently, Pitt & Shephard (Doucet, De Freitas and Gordon (2001))
suggested the use of Auxiliary Particle Filters to overcome some of these
difficulties. An alternative method based upon the Fourier transform has
been presented in Bates (2002).
3.1 The State Space Model
Let us first present the state-space form of the stochastic volatility models:
3.1.1 The Heston Model
Let us study the Euler-discretized Heston (1993) Stochastic Volatility
model in a risk-neutral framework15
lnSk = lnSk-1 + (rk-1 - 1
2
vk-1) t +

vk-1

t Bk-1 (31)
vk = vk-1 + ( - vk-1) t + 

vk-1

t Zk-1 (32)
where Sk is the stock price at time-step k, t the time interval, rk the
risk-free rate of interest (possibly netted by a dividend-yield), vk the
stock variance and Bk , Zk two sequences of temporally-uncorrelated
Gaussian random-variables with a mutual correlation . The model riskneutral
parameter-set is therefore = (, , , ).
Considering vk as the hidden state and lnSk+1 as the observation, we
can subtract from both sides of the transition equation
xk = f(xk-1, wk-1), a multiple of the quantity h(xk-1, uk-1) - zk-1 which is
equal to zero. This would allow us to eliminate the correlation between
the system and the measurement noises.
Indeed, if the system equation is
vk = vk-1 + ( - vk-1) t + 

vk-1

t Zk-1 - [lnSk-1
+ (rk - 1
2
vk-1) t +

vk-1

t Bk-1 - lnSk]
posing for every k
~Zk =
1
1 - 2
(Zk - Bk)
FIGURE 4: Real process versus its estimates: Top EKF, bottom Kushner's NLF
EKF Xk
Kushner's NLF
TABLE 5. PERFORMANCE OF THE
TWO FILTERS
EKF Kushner's NLF
a E(X-X*) Var(X-X*) E(X-X*) Var(X-X*)
0.2 0.01 0.94 0.02 0.92
0.4 -0.1 0.65 -0.06 0.60
0.6 -0.23 0.47 -0.11 0.41
0.8 0.4 0.62 0.14 0.31
1 0.68 1.67 -0.12 0.26
Let x follow an Ornstein-Uhlenbeck or mean-reverting diffusion, like
the convenience yield in the Schwartz model, and z be an exponential:
z = exp(ax) + 0.4 N(0, 1)
The advantage of the example chosen is that we can easily simulate
x and it is therefore easy to assess the performance of each version of
the Kalman filter. It also allows us to make the nonlinearity vary with a.
We first illustrate the effect of an increasing nonlinearity by making
the parameter a vary from 0.2 to 1. The statistics for the expectation and
the variance of the errors have been obtained on a series of 2000 points.
Table 5 shows clearly that when a is small, the bias is small as well.
The latter increases with a. The estimation of the variance is quite high
when the values of a are smaller than or equal to 0.4, because the relative
weight of the noise is large. These estimations of the variance also
increase with the nonlinearity.
To illustrate the performance of the filters and for a better understanding
of the bias, we show on Figure 4 the 200 first points of the
^
Wilmott magazine 13
we will have as expected ~Zk uncorrelated with Bk and
xk = vk = vk-1 + [ - rk -  -
1
2
 vk-1] t
+ ln
Sk
Sk-1
+  1 - 2

vk-1

t ~Zk-1
(33)
and the measurement equation would be
zk = lnSk+1 = lnSk + (rk - 1
2
vk) t +

vk

t Bk (34)
3.1.2 Other Stochastic Volatility Models
It is easy to generalize the above state space model to other stochastic
volatility approaches. Indeed we could replace (32) with
vk = vk-1 + ( - vk-1) t + v
p
k-1

t zk-1 (35)
where p = 1/2 would naturally correspond to the Heston (Square-Root)
model, p = 1 to the GARCH diffusion-limit model, and p = 3/2 to the 3/2
model. These models have all been described and analyzed in Lewis
(2000).
The new state transition equation would therefore become
vk = vk-1 +  - rkv
p- 1
2
k-1 -  -
1
2
v
p- 1
2
k-1 vk-1 t
+ v
p- 1
2
k-1 ln
Sk
Sk-1
+  1 - 2v
p
k-1

t ~Zk-1
(36)
where the same choice of state space xk = vk is made.
3.1.3 Robustness and Stability
In this state-space formulation, we only need to choose a value for v0
which could be set to the historic-volatility over a period preceding our
time-series. Ideally, the choice of v0 should not affect the results enormously,
i.e. we should have a robust system.
As we saw in the previous section, the system stability greatly
depends on the measurement noise. However in this case the system
noise is precisely

tvk, and therefore we do not have a direct control on
this issue. Nevertheless we could add an independent Normal measurement
noise with a given variance R
zk = lnSk+1 = lnSk + (rk - 1
2
vk) t +

vk

tBk + RW k
which would allow us to tune the filter.16
3.2 The Filters
We can now apply the Gaussian and the Particle Filters to our problem:
3.2.1 Gaussian Filters
For the EKF we will have
Ak = 1 - rk p -
1
2
v
p- 3
2
k-1 +  -
1
2
 p +
1
2
v
p- 1
2
k-1 t
+ p -
1
2
v
p- 3
2
k-1 ln
Sk
Sk-1
and
Wk =  1 - 2v
p
k-1

t
as well as
Hk = - 1
2
t
and
Uk =

vk

t
The same time update and measurement update equations could be used
with the UKF or Kushner's NLF.
3.2.2 Particle Filters
We could also apply the Particle Filtering algorithm to our problem.
Using the same notations as in section 1.3.2 and calling
n(x, m, s) =
1

2s
exp((x
- m)2
2s2
)
the Normal density with mean m and standard deviation s, we will have
q(~x
(i)
k |x
(i)
k-1, z1:k) = n ~x
(i)
k , m = ^x
(i)
k , s = P
(i)
k
as well as
p(zk|~x
(i)
k ) = n zk, m = zk-1 + (rk - 1
2
~x
(i)
k ) t, s = ~x
(i)
k

t
and
p(~x
(i)
k |x
(i)
k-1) = n ~x
(i)
k , mx, s =  1 - 2(x
(i)
k-1)p

t
with
mx = x
(i)
k-1 +  - rk(x
(i)
k-1)p- 1
2 -  - 1
2
(x
(i)
k-1)p- 1
2 x
(i)
k-1 t
+ (x
(i)
k-1)p- 1
2 (zk-1 - zk-2)
and as before we have
w
(i)
k = w
(i)
k-1
p(zk|~x
(i)
k )p(~x
(i)
k |x
(i)
k-1)
q(~x
(i)
k |x
(i)
k-1, z1:k)
which provides us with what we need for the filter implementation.
TECHNICAL ARTICLE 1
14 Wilmott magazine
3.3 Parameter Estimation and Back-Testing
For the Gaussian MLE we will need to minimize
(, , , ) =
K
k=1
ln(Fk) +
2
k
Fk
with k = zk - ^zk
and Fk = HkPk
Ht
k + UkUt
k for the EKF.
For the Particle MLE, as previously mentioned, we need to maximize
N
k=1
ln
Nsims
i=1
w
(i)
k
By maximizing the above likelihood functions, we will find the optimal
parameter-set ^ = ( ^, ^, ^, ^). This calibration procedure could then be
used for pricing of derivatives instruments or forecasting volatility.
We could perform a back-testing process in the following way.
Choosing an original parameter-set

= (0.02, 0.5, 0.05, -0.5)
and using a Monte-Carlo simulation, we generate an artificial time-series.
We take S0 = 1000 USD, v0 = 0.04, rk = 0.027 and t = 1. Note that we
are taking a large t in order to have meaningful errors. The time-series
is generated via the above transition equation (32) and the usual log-normal
measurement equation.
We find the following optimal17 parameter-sets:
^ EKF = (0.036335, 0.928632, 0.036008, -1.000000)
^ UKF = (0.033104, 0.848984, 0.033263, -0.983985)
^ EPF = (0.019357, 0.500021, 0.033354, -0.581794)
^ UPF = (0.019480, 0.489375, 0.047030, -0.229242)
which show the better performance of the Particle Filters. However, it
should be reminded that the Particle Filters are also more computationintensive
than the Gaussian ones.
3.4 Application to the S&P500 Index
The above filters were applied to five years of S&P500 time-series (1996 to
2001) in Javaheri (2002) and the filtering errors were considered for the
Heston model, the GARCH model and the 3/2 model. Daily index closeprices
were used for this purpose, and the time interval was set to
t = 1/252. The appropriate risk-free rate was applied and was adjusted
by the index dividend yield at the time of the measurement.
As in the previous section, the performance could be measured via
the MPE and the RMSE. We could also refer to figures 5 to 11 for a visual
interpretation of the performance measurements.
MPEEKF-Heston = 3.58207e - 05 RMSEEKF-Heston = 1.83223e - 05
MPEEKF-GARCH = 2.78438e - 05 RMSEEKF-GARCH = 1.42428e - 05
MPEEKF- 3
2
= 2.63227e - 05 RMSEEKF- 3
2
= 1.74760e - 05
FIGURE 5: Comparison of EKF Filtering errors for Heston, GARCH and 3/2
Models
0
1e-05
2e-05
3e-05
4e-05
5e-05
6e-05
7e-05
0 200 400 600 800 1000 1200
errors
Days
EKF errors for different SV Models
Heston
GARCH
3/2
FIGURE 6: Comparison of UKF Filtering errors for Heston, GARCH and 3/2
Models
0
2e-05
4e-05
6e-05
8e-05
0.0001
0.00012
0 200 400 600 800 1000 1200
errors
Days
UKF errors for different SV Models
Heston
GARCH
3/2
^
Wilmott magazine 15
TECHNICAL ARTICLE 1
FIGURE 7: Comparison of EPF Filtering errors for Heston, GARCH and 3/2
Models
5e-06
1e-05
1.5e-05
2e-05
2.5e-05
3e-05
3.5e-05
4e-05
4.5e-05
5e-05
5.5e-05
0 200 400 600 800 1000 1200
error
Days
EPF errors for different SV Models
Heston
GARCH
3/2
FIGURE 8: Comparison of UPF Filtering errors for Heston, GARCH and 3/2
Models
1.2e-05
1.4e-05
1.6e-05
1.8e-05
2e-05
2.2e-05
2.4e-05
2.6e-05
2.8e-05
0 200 400 600 800 1000 1200
errors
Days
UPF errors for different SV Models
Heston
GARCH
3/2
FIGURE 9: Comparison of Filtering errors for the Heston Model
0
1e-05
2e-05
3e-05
4e-05
5e-05
6e-05
7e-05
200 400 600 800 1000 1200
errors
Days
Heston model errors for different Filters
EKF
UKF
EPF
UPF
FIGURE 10: Comparison of Filtering errors for the GARCH Model
0
1e-05
2e-05
3e-05
4e-05
5e-05
6e-05
7e-05
8e-05
200 400 600 800 1000 1200
errors
Days
GARCH model errors for different Filters
EKF
UKF
EPF
UPF
16 Wilmott magazine
As expected, we observe an improvement when Particle Filters are
used. What is more, given the tests carried out on the S&P500 data it
seems that, despite its vast popularity, the Heston model does not perform
as well as the 3/2 representation.
This suggests further research on other existing models such as Jump
Diffusion (Merton (1976)), Variance Gamma (Madan, Carr and Chang
(1998)) or CGMY (Carr, Geman, Madan and Yor (2002)). Clearly, because of
the non-Gaussianity of these models, the Particle Filtering technique
will need to be applied to them18.
Finally it would be instructive to compare the risk-neutral parameterset
obtained from the above time-series based approaches, to the parameter-set
resulting from a cross-sectional approach using options market
prices at a given point in time19. Inconsistent parameters between the
two approaches would signal either arbitrage opportunities in the market,
or a misspecification in the tested model20.
4 Summary
In this article, we present an introduction to various filtering algorithms:
the Kalman filter, the Extended Filter (EKF), as well as the
Unscented Kalman Filter (UKF) similar to Kushner's Nonlinear filter. We
also tackle the subject of Non-Gaussian filters and describe the Particle
Filtering (PF) algorithm.
We then apply the filters to a term structure model of commodity
prices. Our main results are the following: Firstly, the approximation
introduced in the Extended filter has an influence on the model performances.
Secondly, the estimation results are sensitive to the system
matrix containing the errors of the measurement equation. Thirdly, the
approximation made in the extended filter is not a real issue until the
model becomes highly nonlinear. In that case, other nonlinear filters
such as those described in section 1.2 may be used.
Lastly, the application of the filters to stochastic volatility models
shows that the Particle Filters perform better than the Gaussian ones,
however they are also more expensive. What is more, given the tests carried
out on the S&P500 data it seems that, despite its vast popularity, the
Heston model does not perform as well as the 3/2 representation. This
suggests further research on other existing models such as Jump
Diffusion, Variance Gamma or CGMY. Clearly, because of the nonGaussianity
of these models, the Particle Filtering technique will need to
be applied to them.
Appendix
The Measurement Update equation is
p(xk|z1:k) =
p(zk|xk)p(xk|z1:k-1)
p(zk|z1:k-1)
FIGURE 11: Comparison of Filtering errors for the 3/2 Model
0
1e-05
2e-05
3e-05
4e-05
5e-05
6e-05
7e-05
200 400 600 800 1000 1200
errors
Days
3/2 model errors for different Filters
EKF
UKF
EPF
UPF
MPEUKF -Heston = 3.00000e - 05 RMSEUKF -Heston = 1.91280e - 05
MPEUKF -GARCH = 2.99275e - 05 RMSEUKF -GARCH = 2.58131e - 05
MPEUKF - 3
2
= 2.82279e - 05 RMSEUKF - 3
2
= 1.55777e - 05
MPEEPF-Heston = 2.70104e - 05 RMSEEPF-Heston = 1.34534e - 05
MPEEPF-GARCH = 2.48733e - 05 RMSEEPF-GARCH = 4.99337e - 06
MPEEPF- 3
2
= 2.26462e - 05 RMSEEPF- 3
2
= 2.58645e - 06
MPEUPF-Heston = 2.04000e - 05 RMSEUPF-Heston = 2.74818e - 06
MPEUPF-GARCH = 2.63036e - 05 RMSEUPF-GARCH = 8.44030e - 07
MPEUPF- 3
2
= 1.73857e - 05 RMSEUPF- 3
2
= 4.09918e - 06
Two immediate observations can be made: On the one hand the Particle
Filters have a better performance than the Gaussian ones, which reconfirms
what one would anticipate. On the other hand for most of the
Filters, the 3/2 model seems to outperform the Heston model, which is in
line with the findings of Engle & Ishida (2002).
3.5 Conclusion
Using the Gaussian or Particle Filtering techniques, it is possible to estimate
the stochastic volatility parameters from the underlying asset timeseries,
in the risk-neutral or real-world context.
^
Wilmott magazine 17
where the denominator p(zk|z1:k-1) could be written as
p(zk|z1:k-1) = p(zk|xk)p(xk|z1:k-1)dxk
and corresponds to the Likelihood Function for the time-step k.
Indeed using the Bayes rule and the Markov property, we have
p(xk|z1:k) =
p(z1:k|xk)p(xk)
p(z1:k)
=
p(zk, z1:k-1|xk)p(xk)
p(zk, z1:k-1)
=
p(zk|z1:k-1, xk)p(z1:k-1|xk)p(xk)
p(zk|z1:k-1)p(z1:k-1)
=
p(zk|z1:k-1, xk)p(xk|z1:k-1)p(z1:k-1)p(xk)
p(zk|z1:k-1)p(z1:k-1)p(xk)
=
p(zk|xk)p(xk|z1:k-1)
p(zk|z1:k-1)
Note that p(xk|zk) is proportional to
exp - 1
2
(zk - h(xk))t
R-1
k (zk - h(xk))
under the hypothesis of additive measurement noises.
TECHNICAL ARTICLE 1
13. An alternative approach would consist in choosing a volatility proxy such as f(ln
Sk, ln Sk+1) = ln Sk+1 - ln Sk and ignoring the stock-drift term, given that t =
o(

t). We would therefore write zk = ln(|f(ln Sk, ln Sk+1)|) = E[ln(|f(lnS*k, ln S*k+1
)|)] + 1
2
ln vk + k with St* the same process as St but with a volatility of one, and k
corresponding to the measurement noise. See Alizadeh, Brandt and Diebold (2002)
for details.
14. In this section, all optimizations were made via the Direction-Set algorithm as
described in Press et al. (1997). The precision was set to 1.0e-6.
15. A study on Filtering and Lévy processes has recently been done in Barndorff--Nielsen
and Shephard N. (2002).
16. This idea is explored in At-Sahalia, Wang and Yared (2001) in a Non-parametric
fashion.
17. This comparison supposes that the Girsanov theorem is applicable to the tested
model.
I At-Sahalia Y., Wang Y., Yared F. (2001) "Do Option Markets Correctly Price the
Probabilities of Movement of the Underlying Asset?" Journal of Econometrics, 101
I Alizadeh S., Brandt M.W., Diebold F.X. (2002) "Range-Based Estimation of Stochastic
Volatility Models" Journal of Finance, Vol. 57, No. 3
I Anderson B. D. O., Moore J. B. (1979) "Optimal Filtering" Englewood Cliffs, Prentice
Hall
I Arulampalam S., Maskell S., Gordon N., Clapp T. (2002) "A Tutorial on Particle Filters
for On-line Nonlinear/Non-Gaussian Bayesian Tracking" IEEE Transactions on Signal
Processing, Vol. 50, No. 2
I Babbs S. H., Nowman K. B. (1999) "Kalman Filtering of Generalized Vasicek Term
Structure Models" Journal of Financial and Quantitative Analysis, Vol. 34, No. 1
I Barndorff-Nielsen O. E., Shephard N. (2002) "Econometric Analysis of Realized
Volatility and its use in Estimating Stochastic Volatility Models" Journal of the Royal
Statistical Society, Series B, Vol. 64
I Bates D. S. (2002) "Maximum Likelihood Estimation of Latent Affine Processes"
University of Iowa & NBER
I Brennan M.J., Schwartz E.S. (1985) "Evaluating Natural Resource Investments" The
Journal of Business, Vol. 58, No. 2
I Carr P., Geman H., Madan D., Yor M. (2002) "The Fine Structure of Asset Returns"
Journal of Business, Vol. 75, No. 2
I Doucet A., De Freitas N., Gordon N. (Editors) (2001) "Sequential Monte-Carlo
Methods in Practice" Springer-Verlag
I Engle R. F., Ishida I. (2002) "The Square-Root, the Affine, and the CEV GARCH
Models" Working Paper, New York University and University of California San Diego
I Harvey A. C. (1989) "Forecasting, Structural Time Series Models and the Kalman Filter"
Cambridge University Press
I Harvey A. C., Ruiz E., Shephard N. (1994) "Multivariate Stochastic Variance Models"
Review of Economic Studies, Vol. 61, No. 2
I Haykin S. (Editor) (2001) "Kalman Filtering and Neural Networks" Wiley Inter-Science
I Heston S. (1993) "A Closed-Form Solution for Options with Stochastic Volatility with
Applications to Bond and Currency Options" Review of Financial Studies, Vol. 6, No. 2
I Ito K., Xiong K. (2000) "Gaussian Filters for Nonlinear Filtering Problems" IEEE
Transactions on Automatic Control, Vol. 45, No. 5
I Javaheri A. (2002) "The Volatility Process" Ph.D. Thesis in progress, Ecole des Mines de
Paris
I Julier S. J., Uhlmann J.K. (1997) "A New Extension of the Kalman Filter to Nonlinear
Systems" The University of Oxford, The Robotics Research Group
FOOTNOTES & REFERENCES
1. Some prefer to write xk = f(xk-1 , wk-1). Needless to say, the two notations are
equivalent.
2. The square-root can be computed via a Cholesky factorization combined with a
Singular Value Decomposition. See Press et al. (1997) for details.
3. This analogy between Kushner's NLF and the UKF has been studied in Ito and Xiong
(2000).
4. A description of the Gaussian Quadrature could be found in Press et al. (1997).
5. See Harvey (1989) for an explanation on Smoothing.
6. See for instance Arulampalam et al. (2002) or Van der Merwe, et al. (2000) for details.
7. The same exact methodology could be used in a non risk-neutral setting. We are supposing
we have a small time-step t in order to be able to apply the Girsanov theorem.
8. In that model, interest rates are supposed to be constant.
9. This means that N = 4 in the case we study.
10. In the whole empirical study, optimizations have been made with a precision of 1e-5
on the gradients.
11. For the two filters, the parameters values retained to initiate the optimization are the
same. These values are the following:  = 0.5;  = 0.1; S = 0.3;  = 0.1; C = 0.4;  =
0.5;  = 0.1.
12.The MPE and the RMSE presented here can not be directly compared with those proposed
by Schwartz in 1997 because his computations were done using the logarithm of
the futures prices.
W
18 Wilmott magazine
TECHNICAL ARTICLE 1
I Kushner H. J. (1967) "Approximations to Optimal Nonlinear Filters" IEEE Transactions
on Automatic Control, Vol. 12
I Kushner H. J., Budhiraja A. S. (2000) "A Nonlinear Filtering Algorithm based on an
Approximation of the Conditional Distribution" IEEE Transactions on Automatic Control,
Vol. 45, No. 3
I Lautier D. (2000) "La Structure par Terme des Prix des Commodités : Analyse
Théorique et Applications au Marché Pétrolier" Thése de Doctorat, Université Paris IX
I Lautier D., Galli A. (2000) "Un Mode`le de Structure par Terme des Prix des
Matie`res Premie`res avec Comportement Asymétrique du Rendement d'Opportunité"
Finéco, Vol. 10
I Lewis A. (2000) "Option Valuation under Stochastic Volatility" Finance Press
I Madan D., Carr P., Chang E.C. (1998) "The Variance Gamma Process and Option
Pricing" European Finance Review, Vol. 2, No. 1
I Merton R. C. (1976) "Option Pricing when the Underlying Stock Returns are
Discontinuous" Journal of Financial Economics, No. 3
I Pennacchi G. G. (1991) "Identifying the Dynamics of Real Interest Rates and Inflation :
Evidence Using Survey Data" The Review of Financial Studies, Vol. 4, No. 1
I Press W. H., Teukolsky S.A., Vetterling W.T., Flannery B. P. (1997) "Numerical Recipes
in C: The Art of Scientific Computing" Cambridge University Press, 2nd Edition
I Schwartz E. S. (1997) "The Stochastic Behavior of Commodity Prices : Implications for
Valuation and Hedging" The Journal of Finance, Vol. 52, No. 3
I Van der Merwe R., Doucet A., de Freitas N., Wan E. (2000) "The Unscented Particle
Filter" Oregon Graduate Institute, Cambridge University and UC Berkeley
I Welch G., Bishop G. (2002) "An Introduction to the Kalman Filter" Department of
Computer Science, University of North Carolina at Chapel Hill
I Wells C. (1996) "The Kalman Filter in Finance" Advanced Studies in Theoretical and
Applied Econometrics, Kluwer Academic Publishers, Vol. 32
ACKNOWLEDGEMENT
The second part of the study was realized with the support of the French Institute for
Energy (IFE). The data were provided by TotalFinaElf.