Sequential Monte Carlo for fractional Stochastic Volatility

arXiv:1508.02651v2 [stat.ME] 25 Feb 2017

Models
Alexandra Chronopoulou ∗
Department of Industrial & Enterprise Systems Engineering
University of Illinois at Urbana-Champaign
Konstantinos Spiliopoulos †
Department of Mathematics and Statistics
Boston University
February 28, 2017

Abstract
In this paper we consider a fractional stochastic volatility model, that is a model
in which the volatility may exhibit a long-range dependent or a rough/antipersistent
behavior. We propose a dynamic sequential Monte Carlo methodology that is applicable
to both long memory and antipersistent processes in order to estimate the volatility as
well as the unknown parameters of the model. We establish a central limit theorem
for the state and parameter filters and we study asymptotic properties (consistency
and asymptotic normality) for the filter. We illustrate our results with a simulation
study and we apply our method to estimating the volatility and the parameters of a
long-range dependent model for S& P 500 data.

Keywords: long memory stochastic volatility, rough stochastic volatility, parameter estimation, particle filtering.

1

Introduction

Empirical studies show that the volatility may exhibit correlations of the squared log returns
that decay at a hyperbolic rate, instead of an exponential rate as the lag increases, see for
∗

Research of A.C. supported in part by a start-up fund from the University of Illinois at UrbanaChampaign and by the Simons Foundation.
†
Research of K.S. supported in part by a start-up fund from Boston University and by the National
Science Foundation (DMS 1550918).

1

example [6], [16]. This slow decay cannot be explained by a GARCH stochastic volatility
model or by a stochastic volatility model with jumps. In the literature, this behavior has
been described by a class of models that exhibit long-range dependence in the volatility.
In contrast to the models where the dependence is introduced in the stock returns, [32],
the standard assumption of absence of arbitrage opportunities is preserved. The first long
memory stochastic volatility (LMSV) model was introduced in discrete time by [6] and [26].
In particular, the authors model the stock returns by:

 r =σ
t
t
 σ = σ exp {u } , where (1 − B)d u = η
t

t

t

t

where σ > 0, t are independent and identically distributed (iid) N (0, 1), ηt are iid N (0, ση )2
and independent with t and d ∈ (0, 1/2) is the memory parameter. The feature of longmemory here stems from the fractional difference (1−B)d , where B denotes the lag operator,
i.e. BXt = Xt−1 .
The continuous analogue of the LMSV model was introduced by Comte and Renault in
[10] as a continuous-time mean reverting process in the Hull-White setting. More specifically, their model resembles a classical stochastic volatility model in which the log returns
follow a Geometric Brownian Motion , however, the volatility process is described by a fractional Ornstein-Uhlenbeck process; that is the standard Ornstein-Uhlenbeck process where
the Brownian motion is replaced by a fractional Brownian motion.
A large number of recent papers have considered modeling the volatility in terms of longrange dependent or antipersistent processes. In [11], the authors propose an affine version of
the long memory model in [10] for option pricing and they argue that long-range dependence
provides an explanation for observations of non-flat term structure in long maturities. In
[8], [9], the authors propose a pricing tree algorithm in order to compute option prices in
the discrete and continuous time frameworks. They also propose a calibration method for
determining the memory parameter. In [3], [21], the authors propose a rough fractional
stochastic volatility model, which is used to explain strong skews or smiles in the implied
volatility surface for very short time maturities. In a more recent paper, [20], the authors
discuss the case of small volatility fluctuations of both short and long memory models and
their impact on the implied volatility and derive a corrected Black-Scholes formula.
In this article, we address two problems: filtering of the unobserved volatility and parameter estimation in the case of a stochastic volatility model that is either rough or antipersistent. Specifically, we adapt a Sequential Monte Carlo algorithm in the non-Markovian
framework and then we estimate the parameters online, following an idea initially introduced by [30].
Sequential Monte Carlo (SMC) methods, also known as particle filters or recursive
Monte Carlo filters, arose from the seminal work of Gordon, Salmond and Smith (1993),
2

[22]. SMC methods are iterative algorithms that are based on the dynamic evolution and
update of discrete sets of sampled state vectors, that are referred to as particles, and are
associated with properly selected weights. For an introduction to the SMC literature, there
are several reviews and tutorials, including, the edited volume by Doucet, Freitas and
Gordon (2001), [17] and Kantas et al. (2010), [28].
In order to keep the adaptive nature of the SMC algorithm, we incorporate the parameter estimation procedure in the algorithm by treating the parameters as artificially random,
sampled from a kernel density. Generally speaking, the approach relies on augmenting the
unobserved state by considering the parameter as an unobserved state. We do not assume
that the unknown parameter is fully random, and thus we do not have an additional MCMC
step in the algorithm. We emphasize here that this is an important point, due to the heavy
computational overhead caused by the long-memory issue. Also, as it is well known, using SMC for maximum likelihood based methods for estimating parameters like the Hurst
parameter is difficult.
The output of the algorithm is a provably asymptotically consistent estimator along
with the corresponding variance. One of the appeals of the proposed method is the simplicity in its application, which becomes in particularly important due to the memory issue.
However, we do mention here that overdispersion due to the artificial random evolution of
the parameter is an issue that comes up with such methods, see for example [13]. But, as
we shall see, the overdispersion that the artificial random evolution of the parameter introduces, can be minimized by properly tuning the parameters of the mixture distribution
used for sampling the parameter. Related works using different methods can be found in
[2, 31].
The study of the asymptotic properties of SMC methods and the asymptotic properties
of the filter as the number of particle increases is a quite hard problem. A form of consistency
of the filter, when the number of particles tends to infinity is a common result in the majority
of the literature, while central limit theorem type of results are fewer. One can refer, for
example, to the work of Del Moral and Guionnet (1999), [15]. To the best of our knowledge,
the most general results in the literature can be found in Chopin (2004), [12], Douc and
Moulines (2008), [18], and Künsch (2005), [29]. In these works and under slightly different
conditions, law of large numbers and central limit theorems are derived that apply to most
sequential Monte Carlo techniques in the literature.
The rest of the paper is organized as follows: In Section 2, we introduce the mathematical formulation of the problem. In Section 3, we introduce the SISR (Sequential
Importance Sampling with Resampling) method with parameter learning and we present
the theoretical results for the proposed parameter estimators. In Section 4, we study the
performance of both methods using simulated data. In Section 5, we apply our method

3

in estimating the unobserved volatility of a discrete-time stochastic volatility model with
long-range dependence for S& P 500 data. Finally, we summarize our results in Section 6.

2

Mathematical Framework

Consider a state-space model in which the state vector is denoted by {Xt }t≥1 and the
observations {Yt }t≥1 are obtained sequentially in time. In addition, we assume that the
state vector depends on an unknown, but fixed, parameter vector that we denote by θ.
In the sequel, we use the notation X or Y for the random variables and x or y for the
corresponding realized values.
Unlike other models in the literature, we do not assume that the state vector {Xt }t≥1 is
a Markovian process. Instead, we consider the case in which the unobserved process is not
necessarily Markovian, with particular interest in the long-range dependent case. Formally,
long-range dependence or long-memory is defined as follows:
Definition 1. For a stationary process {Xt }t≥1 , there exists a parameter (Hurst index)
H ∈ ( 12 , 1), such that
Corr(Xt , X1 )
lim
= 1,
(1)
t→∞
ct2−2H
where ρ(h) := Corr(Xt , Xt−h ) is the autocorrelation function of the process.
When H = 21 , then the process is Markovian, so this is a generalization of the models
that are treated in the SMC literature. Equivalently, long-range dependence implies that
the autocorrelation function ρ(h) of a long-range dependent process is non-summable, that
P
is ∞
h=1 ρ(h) = ∞. If the auto-correlation function is summable, then the process has what
is called antipersistence, in which case H ∈ (0, 12 ).
Formally, at time t, the state-space model is specified by the observation equation that
is determined by the observation density
p(yt |xt ; θ)
and the state equation given by the conditional density
p(xt |xt−1 , . . . , x1 ; θ),
where θ ∈ Θ is an unknown vector of parameters, and Θ ⊂ Rd is open. We assume that the
observations {Yt }t≥1 are conditionally independent given {Xt }t≥1 and that the long-range
dependent process {Xt }t≥1 has known initial density X1 ∼ µ(x; θ).
In this article, our goal is to use simulation for online filtering. In other words, we want
to learn about the current state Xt given available information up to time t, which reduces
4

to estimating the probability distribution function
p(xt |y1 , . . . , yt ; θ), t = 1, . . . , n,
where {yt }t≥1 are the observations up to time t. However, since we assume that the parameter θ is unknown, at the same time we also want to estimate θ.

3

Sequential Monte Carlo Filtering

As we mentioned above, apart from filtering for the unobserved states we also want to
estimate the unknown parameter vector θ on which the state vector depends. Our approach
will be to consider θ as an additional state and thus our goal will be to estimate the posterior
distribution {p(x1:t , θ|y1:t )}t≥1 given by
p(x1:t ; θ|y1:t ) ∝ p(x1:t , y1:t , θ)
∝ p(x1 ) · p(x2 |x1 ; θ) · . . . · p(xn |xn−1 , . . . , x1 ; θ) ·

t
Y

p(yi |xi ; θ) · p(θ),

(2)

i=1

where p(θ|yt ) is a prior density for the parameter vector θ, see Subsection 3.1. If the
parameter is known, then the density is degenerate. Therefore, the additional difficulty
here is that we need to compute or approximate the theoretical density function p(θ|yt ).

3.1

On-line Parameter Estimation

One approach in the literature ([19], [22], [34], [30]), is to consider that θ is not fixed and
assume that it artificially evolves in time, for example
θt = θt−1 + et
where et is an artificial white noise with decreasing variance. Then, at each time t, p(θ|yt )
will be updated inside the SISR algorithm in order to incorporate the additional information
that is obtained.
As it was discussed in [30], this approach leads to artificial variance inflation, since the
parameter is not truly random. However, the use of a kernel density estimate with shrinkage
correction can control this artificial over-dispersion.
(j)
More specifically, standing at time t, we approximate p(θ|yt ) by a set of samples θ(t)
(j)

and weights ωt using a discrete Monte Carlo. The index t in θ is in parenthesis to indicate
that θ does not evolve in time, but that its value is drawn from the posterior density p(θ|yt )
at time t. Then, the smooth kernel density with shrinkage correction will be of the form
p(θ|yt ) ≈

N
X



(j)
(j)
wt N θ|mt , h2 Vt ,

j=1

5

where N (·|m, S) denotes a multivariate Normal density with mean m and variance S. So,
(j)
essentially, p(θ|yt ) is approximated by a mixture of normals with mean mt and variance
(j)
h2 Vt , weighted by sample weights ωt . The kernel location is specified by
(j)

mt

(j)

= α θt + (1 − α)θ̄t ,

√
where α = 1 − h2 and θ̄t denotes the average over all parameter samples (essentially a
sum over i). Regarding h, a typical choice would be a decreasing function of the sample
size, but if one wants to control the loss of information then h2 = q − ((3δ − 1)/2δ)2 , where
δ is a discount factor typically around 0.95 − 0.99 and α becomes α = (3δ − 1)/2δ.
Therefore, the SISR algorithm is adjusted in order to incorporate the update of θ.
The key idea of our approach that also allows us to establish asymptotic consistency and
normality of the estimators,
 is to re-formulate the weights so that they represent the joint
(j) (i)
posterior p xt , θ(t) |y1:t and update θ along with the state vector X.
Before presenting the algorithm, let us define the un-normalized weight functions. For
i = 1, . . . , N , and t = 1 set

 

 


(i) (i)
(i)
(i) (i)
(i) (i)


p1 θ(1) µ X1,1 |θ(1) p y1 X̃1,1 , θ(1)
p X1,1 , θ(1) y1
(i)
(i)
 

 
∝

w1 X̃1,1 ; θ(1) = 
(i)
(i)
(i)
(i)
q1 X̃1,1 p1 θ(1)
q1 X̃1,1 p1 θ(1)

 

(i) (i)
(i) (i)
µ X1,1 |θ(1) p y1 X̃1,1 , θ(1)


.
(3)
=
(i)
q1 X̃1,1
Generally for t ≥ 1


(i)
(i)
(i)
(i)
wt = wt X1:t−1,t−1 , X̃t,t ; θ(t)


(i)
(i) (i)
p X1:t−1,t−1 , X̃t,t , θ(t) y1:t
 

= 
(i) (i)
(i)
(i)
(i)
p X1:t−1,t−1 |y1:t−1 ; θ(t) qt X̃t,t , θ(t) |X1:t−1,t−1

 
 

 
(i)
(i)
(i)
(i)
(i)
(i)
(i)
(i) (i)
p yt |X1:t−1,t−1 , X̃t,t , θ(t) p X̃t,t |X1:t−1,t−1 , y1:t−1 , θ(t) p X1:t−1,t−1 |y1:t−1 , θ(t) p θ(t) |y1:t−1 ,
 
  

=
(i)
(i)
(i)
(i)
(i)
(i)
p X1:t−1,t−1 |y1:t−1 ; θ(t) qt X̃t,t X1:t−1,t−1 , θ(t) p θ(t)
 
 


(i)
(i)
(i) (i)
(i)
(i)
(i)
(i)
p yt |X1:t−1,t−1 , X̃t,t , θ(t) p X̃t,t |X1:t−1,t−1 , y1:t−1 , θ(t) p X1:t−1,t−1 |y1:t−1 , θ(t)

 

=
.
(4)
(i)
(i)
(i)
(i)
(i)
p X1:t−1,t−1 |y1:t−1 ; θ(t) qt X̃t,t X1:t−1,t−1 , θ(t)

Then, the algorithm is given by:
At time t = 1
(a) Sampling
(i)
(i)
For i = 1, . . . , N , sample X̃1,1 ∼ q1 (·) and θ(1) ∼ p1 (·).
6

(b) Re-Sampling


(i)
(i)
(i)
For i = 1, . . . , N , let w1 X̃1,1 ; θ(1) be defined by (3), normalize W1 =
P
(i)
such that N
i=1 W1 = 1 and re-sample
(i)
X1,1

∼

N
X

(i)

w1
PN
i=1

(i)

w1

,

(j)

W1 δX̃ (j) (dx1 ) .
(1,1)

j=1

At time t, t ≥ 2 (step t − 1 → t)
(a) Sampling
For i = 1, . . . , N , set
(i)

(i)

X̃t,1:t−1 = Xt−1,1:t−1
(i)

mt−1 = αθ(t−1) + (1 − α)θ̄(t−1)
where θ̄(t−1) =

(i) (i)
i=1 Wt−1 θ(t−1) ,

PN

sample

(i)

(i)

θ(t) ∼ N (mt−1 |h2 Vt−1 ),
and

where Vt−1



(i)
(i)
(i)
X̃t,t ∼ qt · X1:t−1,t−1 ; θ(t)
2
P  (i) (i)
.
= N 1−1 N
i=1 Wt−1 θ(t−1) − θ̄(t−1)

(b) Re-Sampling


(i)
(i)
(i)
(i)
For i = 1, . . . , N , let wt = wt X1:t−1,t−1 , X̃t,t ; θ(t) be defined by (4), and normal(i)
P
(i)
(i)
w
ize Wt = PN t (i) , such that N
i=1 Wt = 1.
i=1

wt

For i = 1, . . . , N , re-sample
X1:t,t ∼

πtN (dx1:t ),

where

πtN (dx1:t )

=

N
X
j=1

and set
θ̄(t) =

N
X

(j)

Wt δX (j)

(j)
1:t−1,t−1 ,X̃t,t

(dx1:t ).

(i) (i)

Wt θ(t)

i=1

Output The filtering distribution p(dx1:t |y1:t ) is approximated by
N

π (dx1:t ) =

N
X
j=1

(j)
Wt δX (j)

(j)
1:t−1,t−1 ,X̃t,t

(dx1:t ),

N
1 X
or π̃ (dx1:t ) =
δX (j) (dx1:t ).
N
1:t,t
N

j=1

7

and the estimator for θ is θ̄(t) . We also record the approximation for the combined distribution π θ (dx1:t , dθ(t) ) = p(dx1:t , dθ(t) |y1:t ) which is approximated by
π N,θ (dx1:t , dθ(t) ) =

N
X

(j)

Wt δX (j)

(j) (j)
1:t−1,t−1 ,X̃t,t ,θ(t)

j=1

3.2

(dx1:t dθ(t) ).

Convergence Results

Let us now study the convergence properties of this algorithm. Let φ : X × Θ 7→ R be
an appropriate test function. Notice now that the SISR algorithm provides us with the
estimator
φ̂N
t =

Z

φ(x1:t , θ(t) )π N,θ (dx1:t dθ(t) ) =

N
X



(i)
(i) (i)
Wti φ X1:t−1,t−1 , X̃t,t , θ(t) .

i=1

It is relatively straightforward to see that φ̂N
t is estimating
Z
φ̄t = φ(x1:t , θ(t) )p(x1:t , θ(t) |Y1:t )dx1:t dθ(t)
where θ(t) is the value that is drawn from the posterior density p(θ|yt ) at time t.
So, it is natural to quantify the performance of the algorithm by studying the convergence of φ̂N
t to φ̄t as N → ∞. We define the set of appropriate test functions under which a
central limit theorem can be established, appropriately formulated for our case of interest:
n
Φt = φ : X 7→ R measurable : there exists δ > 0 such that Eπt kWt φk2+δ < ∞,
o
and x1:(t−1) 7→ Eρ̄(x1:(t−1),· ) (Wt φ)2+δ is in Φt−1 .
(5)
Following the proof of the central limit theorem results of [12, 27] for the Markovian case,
without the parameter estimation aspect, the following result is derived.
Proposition 1. Let us assume that there exists δ > 0 such that for every t < ∞ Eπθ kWt k2+δ <
t
∞ and consider φ ∈ Φt . Then, we get

√  N

N φ̂t − φ̄t ⇒ N 0, σt2 (φ)
as N → ∞, where at time t = 1
σ12 (φ)

Z
=


p2 (x1 , θ(1) |y1 )
 φ(x1 , θ(1) ) − φ̄1 2 dx1 dθ(1)
q1 (x1 )p1 θ(1)

8

and for t > 1
Z
2
Z 2
p (x1 , θ(1) |y1:t )
2

σt (φ) =
φ(x1:t , θ(t) )p(x2:t , θ(t) |y2:t , x1 )dx2:t dθ(t) − φ̄t dx1 dθ(1)
q1 (x1 )p1 θ(1)
t−1 Z
X
p2 (x1:k , θ(k) |y1:t )
+
×
p(x1:(k−1) |y1:(k−1) ; θ(k) )qk (xk , θ(k) |x1:(k−1) )
k=2
Z
2
×
φ(x1:t , θ(t) )p(x(k+1):t , θ(t) |y(k+1):t , x1:k )dx(k+1):t dθ(t) − φ̄t dx1:k dθ(k)
Z
+

2
p2 (x1:t , θ(t) |y1:t )
φ(x1:t , θ(t) ) − φ̄t dx1:t dθ(t) .
p(x1:(t−1) |y1:(t−1) ; θ(t) )qt (xt , θ(t) |x1:(t−1) )

Proof. The proof follows from Proposition A.1.1 of [27] after making the adequate identifications. Indeed, for a general sequential importance sampling algorithm with weights Wt ,
Proposition A.1.1 of [27] implies that the formula for the variance in question is given by
Vt (φ) =

t−1
X

h

i


Varρk Wk Eπθ [φt |X1:k ] − φ̄t + Varρt Wt (φ(x1:t , θ(t) ) − φ̄t ) .
t

k=1

In our case we have
πtθ (dx1:t , dθ(t) ) = p(dx1:t , dθ(t) |y1:t )
ρt (dx1:t , dθ(t) ) = p(dx1:(t−1) |y1:(t−1) )qt (dxt , dθ(t) |x1:(t−1) )
and the weights take the form
Wt (x1:t , θ(t) ) =

p(x1:t , θ(t) |y1:t )
dπtθ
(x1:t , θ(t) ) =
dρt
p(x1:(t−1) |y1:(t−1) )qt (xt , θ(t) |x1:(t−1) )

Plugging these expressions in the formula for Vt (φt ) one immediately recovers the form of
σ 2 (φt ), completing the proof of the proposition.

3.3

Statistical properties of the parameter estimator

Essentially, θ is viewed as an augmented state variable. Proposition 1 quantifies the conN . Let us recall
vergence of the filter, but it does not discuss the statistical properties of θ̄(t)
that
N
X
(i) (N,i)
N
θ̄(t) =
Wt θ(t)
i=1

where
(N,i)

θ(t)

(N,i)

(N,i)

N
∼ N (mt−1 |h2 Vt−1
),
(N,i)

N
mt−1 = αθ(t−1) + (1 − α)θ̄(t−1)
N

N
Vt−1

2
1 X  (i) (N,i)
N
=
Wt−1 θ(t−1) − θ̄(t−1)
N −1
i=1

9

N as
By inspecting the algorithm it becomes clear that the convergence properties of θ̄(t)
N → ∞ is described by a statement very similar to that of Proposition 1 after making
(N )
the choice φ(x1:t , θ(t) ) = θ(t) . Proposition 2 shows that at time t, the estimator for θ, θ̄(t)
R
converges to θ̄(t) = θp(θ(t) |y1:t )dθ as N → ∞.

Proposition 2. Let us assume that there exists δ > 0 such that for every t < ∞ Eπθ kWt k2+δ <
t
∞ and consider the function identity θ(t) 7→ θ(t) , assuming that it belongs to the set of appropriate test functions Φt defined in (5). Let us also define the mean of the posterior
distribution p(θ(t) |y1:t )
Z
θ̄(t) =
Then, we get

√

θp(θ(t) |Y1:t )dθ.




(N )
N θ̄(t) − θ̄(t) ⇒ N 0, σt2 (θ)

as N → ∞. The asymptotic variance σ 2 (θ(t) ) is defined as follows. At time t = 1
σ12 (θ) =

Z


p2 (x1 , θ(1) |y1 )
 θ(1) − θ̄(1) 2 dx1 dθ(1)
q1 (x1 )p1 θ(1)

and for t > 1
σt2 (θ)
+

=

p2 (x1 , θ(1) |y1:t )

q1 (x1 )p1 θ(1)

2

Z
θ(t) p(x2:t , θ(t) |y2:t , x1 )dx2:t dθ(t) − θ̄(t)

dx1 dθ(1)

t−1 Z
X

p2 (x1:k , θ(k) |y1:t )
×
p(x1:(k−1) |y1:(k−1) , θ(k) )qk (xk , θ(k) |x1:(k−1) )
k=2
Z
2
×
θ(t) p(x(k+1):t , θ(t) |y(k+1):t , x1:k )dx(k+1):t dθ(t) − θ̄(t) dx1:k dθ(k)

Z
+

Z

2
p2 (x1:t , θ(t) |y1:t )
θ(t) (x1:t ) − θ̄(t) dx1:t dθ(t) .
p(x1:(t−1) |y1:(t−1) , θ(t) )qt (xt , θ(t) |x1:(t−1) )

(6)

Proof. This proposition is essentially a special case of Proposition 1 with φ(x1:t , θ(t) ) = θ(t) .
We notice that
Z
θ̄(t) = θ(t) p(x1:t , θ(t) |y1:t )dx1:t dθ(t)
Z
= θ(t) p(x1:t |θ(t) , Y1:t )p(θ(t) |y1:t )dx1:t dθ(t)
Z
= θp(θ(t) |y1:t )dθ
which is the mean of the posterior distribution p(θ(t) |y1:t ), as claimed.
The next natural question to ask is whether θ̄(t) is a consistent estimator of θ as t → ∞.
Let us recall from Subsection 3.1 that choosing δ around 0.95 − 0.99 implies for the tuning
10

parameters (α, h) ≈ (1, 0). Subsequently, this means, for every fixed t > 0, that the
(N,i)
(N,i)
mean mt
≈ θt
and the variance h2 Vt ≈ 0, which then implies that for every t > 0,
(N,i)
(N,i)
θt
≈ θt−1 . Notice that if the tuning parameters (α, h) = (1, 0) then the distribution
p(θ(t) |y1:t ) coincides with the posterior distribution p(θ|y1:t ), as the parameter does not
involve artificially in time any more. Hence, in this case, by Doob’s consistency theorem, see
for example Theorem 10.10 in [33], we would have that the sequence of posterior measures
Pθ̄|y1:t is consistent under θ if the model is identifiable, i.e., if Pθ1 6= Pθ2 for θ1 6= θ2 . In other
words, in such a case, for every prior probability measure Π on Θ the sequence of posterior
measures Pθ̄|y1:t is consistent for Π−almost every θ. However, in the algorithm and for the
purposes of dealing with the issue of degeneracy of the particles, we artificially evolve the
unknown parameter which means that we take the tuning parameters to be (α, h) ≈ (1, 0)
but not exactly equal to (1, 0). As we shall see from the simulation studies of Section 4, this
is sufficient to guarantee that the parameter is being consistently estimated, as expected.

4
4.1

Simulation Results
Fractional ARIMA process

The fractional ARIMA (AutoRegressive Integrated Moving Average) process was proposed
by Box and Jenkins, [5], in 1970 and has been very popular in applied time series. A
fractional ARIMA(p, d, q) process is formally defined as follows (due to Granger and Joyeux,
[23]):
Definition 2. Let ϕ(·) and ϑ(·) be polynomials of orders p and q respectively and Xt a
stationary process such that
ϕ(B)(1 − B)d Xt = ϑ(B)ηt .
d ∈ (−1/2, 1/2) and and (ηt )t≥0 is a sequence of iid variables with mean 0 and variance 1.
Then, the process {Xt }t≥0 is called a fractional ARIMA(p, d, q) process.
In contrast to the classical ARIMA(p, d, q) process, where the parameter d is an integer,
in the fractional case d is a real valued parameter with values between (−1/2, 1/2). It is
called the fractional integration parameter and is related to the Hurst index, H in (1), via
d = H − 12 . B denotes the lag or backshift operator, and
∞  
X
d
d
(1 − B) =
(−1)k B k ,
k
k=0

where the sum is taken over an infinite number of indices. The fractional ARIMA process
is long-range dependent when d > 0, while the upper bound on d is needed to ensure that
the process is stationary. More details regarding these models can be found in Beran [4].
11

In our framework, we consider a state-space model in which the unobserved process
is modeled by a Fractional ARIMA(p, d, q) process. Specifically, the state-space model is
defined as follows

 Y = σ Xt  
t
t
2
(7)
 ϕ(B) (1 − B)d X = ϑ(B) η ,
t

t

where t and ηt are two independent iid sequences of Gaussian random variables and σ(·)
is a known function.
4.1.1

SISR for Fractional ARIMA process with known parameter

We apply our algorithm to simulated data from an ARIMA(1, 0.3, 0) model with parameter
ϕ = 0.8. That is

 Y = |X | 
t
t t
(8)
 (1 − ϕ B) (1 − B)d X = η ,
t
t
The simulated model is shown in Figure 1(a) and the estimated filter using the SISR
algorithm is depicted in Figure 1(b). We choose as tuning parameters δ = 0.98 which then
implies (α, h) = (0.9506, 0.096).

(a) Simulated Model

(b) SISR Filter

Figure 1: SISR filter for the fractional ARIMA(1, 0.3, 0) model. Here, the parameter ϕ is
assumed to be known, ϕ = 0.8.
From the SISR plot (Figure 1b), we can see that the approximation of the unobserved
process (solid line) follows closely the true process (dashed line). In addition, the true data,
all lie within the 95% confidence interval (dotted lines) estimated using the SISR algorithm.
In our simulation study, we tried different values of d and the results are similar.

12

4.1.2

SISR for Fractional ARIMA process with unknown parameter

Consider again model (8), but now assume that the parameter ϕ is unknown. Our goal is
to estimate ϕ using the SISR algorithm. The results are summarized in Figures 4 and 3.
From Figure 4, we can see that the approximation of the unobserved process remains good,
and the 95% confidence interval still captures the process.

Figure 2: SISR filter for the ARIMA(1, 0.3, 0) model. In this case, the parameter ϕ is
assumed to be unknown and is estimated from the algorithm.
In Figure 3, we investigate the convergence of the parameter to the true value. The
estimated parameter ϕ is slightly noisy, but it converges to the true value 0.8. The difference
in the two graphs in Figure 3, is that the second one has a significantly larger number of
simulated particles, and it seems that this slightly improves the smoothness of the curve.
In Figure 4, we compare the empirical variance of the estimator for N = 500 versus
N = 2500 across all t ∈ [1, 1000]. It is clear that the variance decreases as N increases, but
at the same time the variance does increase over time. The latter is consistent with the
theoretical limiting variance as obtained in (6).

5

Application to S&P 500 Data

In this section, we apply our method to real data. As an example, we are working with a
long-range dependent state-space model in finance. The observed process are the returns of
the underlying asset (S& P 500 index to be precise) and the unobserved process is the asset’s
volatility. Based on the financial literature ([6, 26, 9]), we assume that the volatility process
is long-range dependent, and the model we focus on is the long memory stochastic volatility
model in discrete time that is described by (7). This model was introduced simultaneously
by Breidt et al. [6] and Harvey, [26] in 1993.

13

(a) Number of Particles N =500

(b) Number of Particles N =2,500

Figure 3: Convergence of the estimated parameter as a function of time, for two different
choices of number of particles.

(a) Number of Particles N =500

(b) Number of Particles N =2,500

Figure 4: Empirical Variance of the estimated parameters as a function of time for N=500
and N=2,500
To further specify this model in practice, we need to determine the order of the polynomials φ(·) and θ(·). This is a common task in time series analysis and for details we refer
to Hamilton [24]. Based on a preliminary analysis, we choose to work with a Fractional
ARIMA(1, d, 1) model, which is also in accordance to the model suggested by [1]. To be
precise, the model we will be working with is

 Y = σ Xt  
t
t
2
 (1 − ϕB) (1 − B)d X = ϑ η
+η ,
t

t−1

t

where σ(x) = log x.
The data set we consider contains daily returns of the S&P 500 for one year, that is
14

about 252 observations, starting in January 2010 until December 2010.
One assumption that we made in the SMC algorithm is that the parameter d is known.
However, when it comes to real data, this is something that we need to estimate. Here, we
used the Geweke and Porter-Hudak estimator, [7], which yields d = 0.2. Then, we apply
the SISR algorithm to estimate the remaining unknown parameters of the model ϕ and ϑ.
We choose as tuning parameters δ = 0.986 which then implies (α, h) = (0.965, 0.068).
The algorithm has two outputs. The first one is the distribution of the unobserved
volatility, which is given in Figure 5, using 500 trajectories, and the second one is the
estimated vector of parameters, which are plotted as a function of time in Figure 6.

Figure 5: Histogram of the empirical distribution of the unobserved volatility. As a reference, the implied volatility for the same period was 0.355.
Using the model parameters we estimated, we also do an out-of-sample prediction of
the values of the underlying asset, which is shown in Figure 7. By doing an one-step
prediction each time, we forecast 20 daily values of the index. The 95% confidence intervals
are computed using boostrap. In Figure 7 we also present the empirical variance for the
estimators ϕ̂ and θ̂. In Figure 8, we also present the residuals of the fitted model.

6

Conclusion

To summarize, in this article we extended the standard SISR algorithm to incorporate the
case that the observations are long-range dependent. Our findings show that the results
are very close to the case that the observations are independent or Markov. However,
the main drawback of this method is the computational time that is required to perform
the iterations. Since we need to take into account, and technically speaking to store all
past values of the trajectory, this increases the computational time and complexity of the
method. In addition, by naturally extending existing results in the literature, we proved
that the filter converges to the true distribution of the unobserved process.
Our second outcome, was the development of an SISR algorithm that along with the
estimation of the unobserved distribution of the hidden process, also estimated unknown
15

]

(a) Estimator of ϕ

(b) Estimator of ϑ

(c) Empirical Variance of ϕ̂

(d) Empirical Variance of ϑ̂

Figure 6: Estimators for ϕ and ϑ for the S&P 500 data set. Solid lines represent the
estimated values for ϕ̂ = 0.842 and θ̂ = 0.01 for φ and θ respectively.
model parameters. Our approach was dynamic, in the sense that the parameter was regarded as “time-varying” and thus the parameter estimators were updated at every step of
the algorithm. We also showed that the proposed estimators for the unknown parameter are
consistent and asymptotically normal and we corroborated these results with a simulation
study.
There are quite a few open problems that we would like to investigate in the future.
The first one is to study ways to improve the computational efficiency of the algorithm.
In our approach, we used all the history of the trajectory to run the algorithm, which
severely affected the computational efficiency, but it would be interesting to investigate if
a “window” approach would provide us with a reasonable estimator for the filter and/or
the parameter, and possibly quantify the loss that one might have by doing so in terms of
accuracy.
In addition, one question that we did not address in this paper, is what happens with the
16

Figure 7: Out-of-Sample Prediction of the S& P 500 values. The (black) solid continuous
line are the true S& P 500 values, while the (red) solid step line is our estimation. The
(red) dotted lines form the 95% prediction interval.

]

(a) Standard Residuals

(b) ACF of the Residuals

Figure 8: Residuals of the fitted fractional ARIMA (1, 0.2, 1) model with ϕ = 0.842 and
ϑ = 0.01.
long memory parameter in practice. In our approach, we assumed that d (or equivalently
H) is known (given or estimated from the data). However, it is an open question how
one would consistently estimate the memory parameter in the scenario that the long-range
dependent process is hidden.
The goal of this first paper on the topic is to lay down the algorithm, its properties
and to understand the main issues that the presence of memory brings. In this paper, we
formulated the algorithm and established some baseline theoretical properties. We also
performed numerical experiments using both simulation data and real data, in order to test
the performance of the proposed algorithm in practice. We plan to investigate the open
computational and theoretical issues mentioned above in a systematic way in future works.
17

References
[1] Baillie, R. T., Bollerslev, T., and Mikkelsen, H. O. (1996). Fractionally Integrated
Generalized Autoregressive Conditional Heteroskedasticity. J. Econom. 74(1):3–30.
[2] Beskos, A., Dureau, J, and Kalogeropoulos, K. (2013), Bayesian inference for partially
observed SDEs driven by fractional Brownian motion, arXiv: 1307.0238.
[3] Bayer, C., Friz, P. and Gatheral, J. (2015) Pricing under rough volatility. Preprint.
[4] Beran, J. (1994): Statistics for Long-Memory Processes. Chapman and Hall.
[5] Box, G.E.P. and Jenkins (1970) Time series analysis: forecasting and control. Holden
Day, San Fransisco.
[6] Breidt, F. J., Crato, N. and De Lima, P. (1998). The Detection and Estimation of
Long-Memory in Stochastic Volatility. Journal of Econometrics. 83:325–348.
[7] Casas, I. and Gao, J. (2008) Econometric estimation in long-range dependent volatility
models: Theory and practice. Journal of Econometrics 147:72–83.
[8] Chronopoulou, A. and Viens, F. (2012). Stochastic volatility and option pricing with
long-memory in discrete and continuous time. Quantitative Finance 12(4):635–649.
[9] Chronopoulou, A. and Viens, F. (2010). Estimation and pricing under long-memory
stochastic volatility. Annals of Finance, 8 (2-3):379–403.
[10] Comte, F. and Renault, E. (1998). Long Memory in Continuous-time Stochastic
Volatility Models. Mathematical Finance, 8(4):291–323.
[11] Comte, F., Coutin, L. and Renault, E. (2012) Affine fractional stochastic volatility
models. Ann Finance 8: 337.
[12] Chopin, N., (2004) Central limit theorem for sequential Monte Carlo method and its
application to Bayesian inference. The Annals of Statistics, 32(6):2385–2411.
[13] Chopin, N., Jacob, P.E., and Papaspiliopoulos, O. (2013) SMC2 : A sequential Monte
Carlo algorithm with particle Markov chain Monte Carlo updates. J. R. Statist. Soc.
B, 75(3), 397–426.
[14] Del Moral, P. (2004), Feynman-Kac Formulae: Genealogical and interacting particle
systems with applications, In Series Probability and Applications, Springer-Verlag, New
York.

18

[15] Del Moral, P. and Guionnet, A. (1999) Central limit theorem for nonlinear filtering
and interacting particle systems. Annals of Applied Probability, 9:275–297.
[16] Ding,Z., C., Granger, W., J. and Engle, R. F. (1993). A Long Memory Property of
Stock Market Returns and a New Model. J. Empirical Finance 1, 1.
[17] Doucet, A., De Freitas, J.F.G. and Gordon, N.J. (eds.) (2001). Sequential Monte Carlo
Methods in Practice. New York: Springer-Verlag.
[18] Douc, R. and Moulines, E. (2008), Limit theorems for weighted samples with applications to sequential Monte Carlo methods, Annals of Statistics, 36(5):2344–2376.
[19] Flury, T. and Shephard, N. (2011). Bayesian inference based only on simulated likelihood: particle filter analysis of dynamic economic models. Econometric Theory,
27:933–956.
[20] Garnier, J. and Solna, K. (2015). Correction to Black-Scholes formula due to fractional
stochastic volatility. Preprint.
[21] Gatheral, J., Jaisson, Th., and Rosenbaum, M. (2014) Volatility is rough. Preprint.
[22] Gordon, N., Salmond, D., and Smith, A. (1993) Novel approach to non-linear/nonGaussian Bayesian state estimation. IEE Proceedings-F, 140:107–113.
[23] Granger, C.W., Joyeux, R., (1980). An introduction to long-memory time series models
and fractional differencing. Journal of Time Series Analysis 1:15–29.
[24] Hamilton, J.D. (1994) Time Series Analysis. Princeton University Press.
[25] Handschin, J. and Mayne, D. (1969), Monte Carlo techniques to estimate the conditional expectation in multi-stage non-linear filtering. International Journal of Control,
9:547–559.
[26] Harvey, A.C. (1998). Long Memory in Stochastic Volatility, in J.Knight and S. Satchell
(eds.) Forecasting Volatility in Financial Markets, Butterworth-Haineman, Oxford, pp.
307-320.
[27] Johansen, A. and Doucet, A., Auxiliary variable sequential Monte Carlo methods.
Statistics group technical report, 07:09, University of Bristol, July 2007.
[28] Kantas, N., Doucet, A., Singh, S. S., Maciejowski, J. M., Chopin, N. (2014), An
overview of sequential Monte Carlo methods for parameter estimation in general statespace models. Submitted to Statistical Science.

19

[29] Künsch, H.R. (2005), Recursive Monte Carlo filters: algorithms and theoretical analysis, Annals of Statistics, 33(5): 1983–2021.
[30] Liu, J. and West, M., Combined parameter and state estimation in simulation-based
filtering, In A. Doucet, N. De Freitas and N. Gordon (Eds.), Sequential Monte Carlo
methods in practice, New York: Springer, 197–223.
[31] Lysy, M., and Pillai, N.S., (2013), Statistical inference for stochastic differential equations with memory, arXiv:1307.1164.
[32] Rogers, L. C. G. (1997) Arbitrage with fractional Brownian motion, Mathematical
Finance, 7 , 95–105.
[33] Van der Vaart, A.W. (1998), Asymptotic Statistics, Cambridge series in statistical and
probabilistic mathematics, Cambridge University Press.
[34] West (1993), M. Approximating posterior distributions by mixtures. Journal of the
Royal Statistical Society (series B), 55:409–422.

20

