arXiv:1804.00849v1 [math.ST] 3 Apr 2018

Robust estimation of continuous-time
ARMA models via indirect inference
Vicky Fasen-Hartmann †‡

Sebastian Kimmig §

In this paper we present a robust estimator for the parameters of a continuous-time
ARMA(p, q) (CARMA(p, q)) process sampled equidistantly which is not necessarily
Gaussian. Therefore, an indirect estimation procedure is used. It is an indirect estimation because we first estimate the parameters of the auxiliary AR(r) representation
(r ≥ 2p − 1) of the sampled CARMA process using a generalized M- (GM-)estimator.
Since the map which maps the parameters of the auxiliary AR(r) representation to the
parameters of the CARMA process is not given explicitly, a separate simulation part is
necessary where the parameters of the AR(r) representation are estimated from simulated
CARMA processes. Then, the parameter which takes the minimum distance between the
estimated AR parameters and the simulated AR parameters gives an estimator for the
CARMA parameters. First, we show that under some standard assumptions the GMestimator for the AR(r) parameters is consistent and asymptotically normally distributed.
Next, we prove that the indirect estimator is consistent and asymptotically normally distributed as well using in the simulation part the asymptotically normally distributed LSestimator. The indirect estimator satisfies several important robustness properties such as
weak resistance, πdn -robustness and it has a bounded influence functional. The practical
applicability of our method is demonstrated through a simulation study with replacement
outliers and compared to the non-robust quasi-maximum-likelihood estimation method.
AMS Subject Classification 2010: Primary: 62F10, 62F12, 62M10
Secondary: 60G10, 60G51
Keywords: AR process, CARMA process, indirect estimator, influence functional, GM-estimator,
LS-estimator, outlier, resistance, robustness

1 Introduction
The paper presents a robust estimator for the parameters of a discretely observed continuous-time
ARMA (CARMA) process. A weak ARMA(p, q) process in discrete-time is a weakly stationary
† Institute

of Stochastics, Englerstraße 2, D-76131 Karlsruhe, Germany.
by the Deutsche Forschungsgemeinschaft through the research grant FA 809/2-2 is gratefully acknowl-

‡ Financial support

edged.
§ Württembergische

Versicherung AG, Gutenbergstraße 30 D-70176 Stuttgart, Germany.

1

solution of the stochastic difference equation

φ (B)Xm = θ (B)Zm ,

m ∈ Z,

(1.1)

where B denotes the backward shift operator (i.e. BXm = Xm−1 ),

φ (z) = 1 − φ1 z − . . . − φ p z p

and

θ (z) = 1 + θ1 z + . . . + θq zq

are the autoregressive and the moving average polynomial, respectively, with φ1 , . . . , φ p , θ1 , . . . , θq ∈
R, φ p , θq 6= 0 and (Zm )m∈Z is a weak white noise. If (Zm )m∈Z is even an independent and identically
distributed (i.i.d.) sequence then we call (Xm )m∈Z a strong ARMA process. A natural continuoustime analog of this difference equation with i.i.d. noise (Zm )m∈Z is the formal p-th order stochastic
differential equation
a(D)Yt = c(D)DLt ,

t ∈ R,

(1.2)

where D denotes the differential operator with respect to t,
a(z) = z p + a1 z p−1 + . . . + a p

and

c(z) = c0 zq + c1 zq−1 + . . . + cq

are the autoregressive and the moving average polynomial, respectively, with p > q, and
a1 , . . . , a p , c0 , . . . , cq ∈ R, a p , c0 6= 0. The process (Lt )t∈R is a Lévy process, i.e., a stochastic process
with L0 = 0 almost surely, independent and stationary increments and almost surely càdlàg sample
paths. However, this is not the formal definition of a CARMA(p, q) process because a Lévy process is not differentiable. Going back to Brockwell (2001) a CARMA process can be defined via
a controller canonical state space representation. Necessary and sufficient conditions for the existence of strictly stationary CARMA processes are given in Brockwell and Lindner (2009). From
Brockwell and Lindner (2009) it is well known as well that a discretely sampled CARMA process
(Ymh )m∈Z (h > 0 fixed) admits a weak ARMA representation, but unfortunately this is in general for
Lévy driven models not a strong ARMA representation. For an overview and a comprehensive list of
references on CARMA processes we refer to Brockwell (2014).
In statistics the most fundamental question when considering robustness of an estimator is how the
estimator behaves when the data does not satisfy the model assumptions (cf. Huber and Ronchetti
(2009); Maronna et al. (2006); Olive (2017)). In the case of small deviations from the model assumptions a robust estimator should give estimations not too far away from the estimations of the original
model. The most common and best understood robustness property is distributional robustness where
the shape of the true underlying distribution deviates slightly from the assumed model. The amount of
measures for robustness is huge, e.g., qualitative robustness, quantitative robustness, optimal robustness, efficiency robustness and the breakdown point, to mention only a few. In contrast to the case of
i.i.d. random variables, in the case of time series, there exist several types of possible contamination
of the data which makes it more difficult to characterize robustness. In particular, for AR processes
it is well-known that the GM-estimator (cf. Boente et al. (1987); Künsch (1984); Martin (1980)) and
the RA-estimator (cf. Ben et al. (1999)) satisfy different robustness properties in contrast to M- or
LS-estimators which are sensitive to the presence of additive outliers (cf. Denby and Martin (1979)).
However, for general ARMA models the GM-estimator and the RA-estimator are again sensitive to
outliers and hence, non-robust (cf. Bustos and Yohai (1986)). Muler et al. (2009) develop a robust
estimation procedure for ARMA models by calculating the residuals of the ARMA models with the
help of BIP-ARMA models. Important for their results is that they have a strong ARMA model.

2

Unfortunately the results can not easily be extended to weak ARMA models which we have in our
context.
There exist a few papers dealing with the estimation of the parameters of a discretely sampled
CARMA process as Schlemm and Stelzer (2012); Brockwell et al. (2011); Fasen and Kimmig (2017);
Fasen-Hartmann and Scholz (2017). The papers have in common that they use a quasi maximum
likelihood estimator (QMLE). However, it is well known that a QMLE is sensitive to outliers and
irregularities in the data. Hence, we are looking for an alternative robust approach.
In this paper we use the indirect inference method, originally proposed by Smith (1993) and extended by Gouriéroux et al. (1993) (see also the overview in Gouriéroux and Monfort (1997)) for
models with intractable likelihood function. The core idea of the indirect estimation method is to avoid
estimating the parameter of interest directly and instead fit an auxiliary model to the data, estimate
the parameters of this auxiliary model and then use this estimates with simulated data to construct
an estimator for the original parameter of interest (see de Luna and Genton (2001) for a schematic
overview over the indirect estimation method). This method has been successfully used in different
contexts, see, e.g., Gouriéroux and Monfort (1997); Jiang and Turnbull (2004); de Luna and Genton
(2001, 2000). The latter two papers recognized that it is possible to construct robust estimators via
this approach, even for model classes where direct robust estimation is difficult. The reason is that it is
sufficient if the parameters of the auxiliary model are estimated by a robust estimation method. Therefore, de Luna and Genton (2001) present an indirect estimation procedure for strong ARMA processes
(without detailed assumptions and rigorous proofs). They fit an AR(r) process to the ARMA model
and estimate the parameters of the AR(r) process with a GM-estimator. We present a similar ansatz
in our paper for the estimation of the CARMA parameters. Since the discretely sampled CARMA
process admits a weak ARMA representation instead of a strong ARMA representation several proofs
have to be added and identifiability issues have to be taken into account.
The paper is structured as follows. In Section 2 we first present our parametric family of CARMA
processes and our model assumptions. Furthermore, we motivate that for any r ≥ 2p − 1 any CARMA
process has an AR(r) representation. Then, in Section 3, we introduce the indirect estimation procedure and give sufficient criteria for indirect estimators to be consistent and asymptotically normally distributed independent of the model; we have to assume at least consistent and asymptotically
normally distributed estimators in the estimation part and in the simulation part of the indirect estimation. Since the auxiliary AR(r) parameters of the sampled CARMA process are estimated by a
GM-estimator we give an introduction into GM-estimators in Section 4 and derive consistency and
asymptotically normality of this estimator in our setup. Moreover, we see that the GM-estimator is
still asymptotically normally distributed for CARMA processes with outliers as additive outliers and
replacement outliers. Our conclusions extend the results of Bustos (1982). Finally, in Section 5, we
are able to show that the indirect estimator for the parameters of the discretely observed CARMA process is consistent and asymptotically normally distributed using in the estimation part a GM-estimator
and in the simulation part a LS-estimator. Robustness properties of this estimator are topic of Section 6. In particular, qualitative robustness properties and the influence functional are studied. After
all, the simulation part, in Section 7, shows the practical applicability of our indirect estimator and its
robustness properties. We compare our estimator with the non-robust QMLE. The paper ends with a
conclusion in Section 8.
The present paper is mainly based on the PhD-thesis Kimmig (2016). The same arguments for
the indirect estimation as in Kimmig (2016) were later applied in do Rêgo Sousa et al. (2017) for the
estimation of the COGARCH process using only the LS- and the Yule-Walker-estimator and without
investigating any robustness properties.

3

Notation
We use as norms the Euclidean norm k·k in Rd and its operator k·k in Rm×d which is submultiplicative.
For a matrix A ∈ Rm×d we denote by AT its transpose. For a matrix function f (ϑ ) in Rm×d with
f (ϑ ))
ϑ ∈ Rs the gradient with respect to the parameter vector ϑ is ∇ϑ f (ϑ ) = ∂ vec(
∈ Rdm×s and
∂ϑT
D

P

ϑ f (ϑ ))
∈ Rdms×s . Finally, we write −→ for weak convergence and → for
similarly ∇2ϑ f (ϑ ) = ∂ vec(∇
∂ϑT
convergence in probability. In general C denotes a constant which may change from line to line.

2 Preliminaries
2.1 The CARMA model
In this paper we consider a parametric family of CARMA processes. Let Θ ⊆ RN(Θ ) (N(Θ ) ∈ N)
be a parameter space, p ∈ N be fixed and for any ϑ ∈ Θ let a1 (ϑ ), . . . , a p (ϑ ), c0 (ϑ ), . . . , c p−1 (ϑ ) ∈
R, a p (ϑ ) 6= 0 and c j (ϑ ) 6= 0 for some j ∈ {0, . . . , p − 1}. Furthermore, define

...
0

.. 
..
 0
.
. 
0
1


 ∈ R p×p ,

..
.. ..
Aϑ :=  ..

.
.
0
.
.


 0
0
... 0
1 
−a p (ϑ ) −a p−1 (ϑ ) . . . . . . −a1 (ϑ )
q(ϑ ) = sup{ j ∈ {0, . . . , p − 1} : cl (ϑ ) = 0 ∀ l > j} with sup 0/ := p − 1,


cϑ

0

1

0

:= (cq(ϑ ) (ϑ ), cq(ϑ )−1 (ϑ ), . . . , c0 (ϑ ), 0, . . . , 0)T ∈ R p .

Due to Brockwell (2001) the CARMA process (Yt (ϑ ))t∈R is then defined via the controller canonical
state space representation: Let (Xt (ϑ ))t∈R be a strictly stationary solution to the stochastic differential
equation
dXt (ϑ ) = Aϑ Xt (ϑ ) dt + e p dLt , t ∈ R,
(2.1a)
where e p denotes the p-th unit vector in R p . Then the process
Yt (ϑ ) := cϑT Xt (ϑ ),

t ∈ R,

(2.1b)

is said to be a CARMA process of order (p, q(ϑ )). This means that in our parametric family of
CARMA processes the order of the autoregressive polynomial is fixed to p but the order of the moving
average polynomial q(ϑ ) may change.
Furthermore, we have the discrete-time observations Yh , . . . ,Ynh of the CARMA process (Yt )t∈R =
(Yt (ϑ0 ))t∈R with fixed grid distance h > 0. Hence, the true model parameter is ϑ0 . The aim is to
receive from the observations Yh , . . . ,Ynh an estimator for ϑ0 . Throughout the paper we will assume
that the following Assumption A holds.
Assumption A.
(A.1) The parameter space Θ is a compact subset of RN(Θ ) .
(A.2) The true parameter ϑ0 is an element of the interior of Θ .
(A.3) E[L1 ] = 0, 0 < E|L1 |2 = σL2 < ∞ and there exists a δ > 0 such that E|L1 |4+δ < ∞.

4

(A.4) The eigenvalues of Aϑ have strictly negative real parts.
(A.5) For all ϑ ∈ Θ the zeros of cϑ (z) = c0 (ϑ )zq(ϑ ) + c1 (ϑ )zq(ϑ )−1 + . . . + cq(ϑ ) are different from
the eigenvalues of Aϑ .
(A.6) For any ϑ , ϑ ′ ∈ Θ we have (cϑ , Aϑ ) 6= (cϑ ′ , Aϑ ′ ).
(A.7) For all ϑ ∈ Θ the spectrum of Aϑ is a subset of {z ∈ C : − πh < Im(z) < πh } where Im(z) denotes
the imaginary part of z.
(A.8) The maps ϑ 7→ Aϑ and ϑ 7→ cϑ are three times continuous differentiable.
Remark 2.1.
(i) (A.1) and (A.2) are standard assumptions in point estimation theory.
(ii) (A.4) guarantees that there exists a stationary solution of the state space model (2.1a) and hence,
of the CARMA process (Yt (ϑ ))t∈R (see Marquardt and Stelzer (2007)). For this reason we will
assume throughout the paper that (Yt (ϑ ))t∈R is stationary.
(iii) A consequence of (A.4), (A.8), the compactness of Θ and the fact that the eigenvalues of a matrix are continuous functions of its entries (cf. (Bernstein, 2009, Fact 10.11.2)) is supϑ ∈Θ max{|λ | :
λ is eigenvalue of eAϑ } < 1 and hence, supϑ ∈Θ k eAϑ u k ≤ C e−ρ u for some C, ρ > 0.
(iv) Due to (A.5) the state space representation (2.1) of the CARMA process is minimal (cf. (Bernstein,
2009, Proposition 12.9.3) and (Hannan and Deistler, 2012, Theorem 2.3.3)).
(v) A consequence of (A.5) and (A.6) is that the family of CARMA processes (Yt (ϑ ))t∈R is identifiable from the spectral density and in combination with (A.7) that the same is true for the
discrete-time process (Ymh (ϑ ))m∈Z (cf. (Schlemm and Stelzer, 2012, Theorem 3.13)).
In the following we denote the autocovariance function of (Yt (ϑ ))t∈R as (γϑ (t))t∈R which has by
(Schlemm and Stelzer, 2012, Proposition 3.1) the form

γϑ (t) = cTϑ eAϑ t Σϑ cϑ ,
R

t ∈ R,

(2.2)

with Σϑ = σL2 0∞ eAϑ u e p eTp eAϑ u du. Due to Assumption A the autocovariance function is three times
continuous differentiable as well.

2.2 The AR(r) representation of a CARMA process
First, we define the auxiliary AR(r) representation of the sampled CARMA process (Ymh (ϑ ))m∈Z .
Proposition 2.2. For every ϑ ∈ Θ and every r ≥ 2p − 1, there exists a unique

π (ϑ ) := (π1 (ϑ ), . . . , πr (ϑ ), σ (ϑ )) ∈ Rr × [0, ∞)
such that

r

Um (ϑ ) := Ymh (ϑ ) − ∑ πk (ϑ )Y(m−k)h (ϑ )
k=1

5

(2.3)

is stationary with E[U1 (ϑ )] = 0, Var(U1 (ϑ )) = σ 2 (ϑ ) and


E Um (ϑ )Y(m−k)h (ϑ ) = 0 for k = 1, . . . , r.

(2.4)

We call π (ϑ ) the auxiliary parameter of the AR(r) representation of (Ymh (ϑ ))m∈Z .
Proof. First, we need to show that for any r ∈ N the covariance matrix of (Yh (ϑ ), . . . ,Y(r+1)h (ϑ ))
is non–singular. To see this, note that the autocovariance function of (Ymh (ϑ ))m∈Z is γϑ (mh) =
cTϑ eAϑ hm Σϑ cϑ , m ∈ Z (see (2.2)). Since Σϑ is non–singular (cf. (Schlemm and Stelzer, 2012, Corollary 3.9)) and cϑ 6= 0 p we have that γϑ (0) > 0. Moreover, the eigenvalues of Aϑ have strictly negative
real parts by (A.4) and therefore, γϑ (mh) → 0 as m → ∞ holds. By (Brockwell and Davis, 1991,
Proposition 5.1.1), it follows that the covariance matrix of (Yh (ϑ ), . . . ,Y(r+1)h (ϑ )) is non–singular
for every r ∈ N. Thus, a conclusion of (Brockwell and Davis, 1991, §8.1) is that there exist unique
π1 (ϑ ), . . . , πr (ϑ ), σ 2 (ϑ ) which solve the set of r + 1 Yule–Walker equations, namely



π1 (ϑ )


π ∗ (ϑ ) :=  . . .  = 

πr (ϑ )

γϑ (0)
γϑ (h)
..
.

γϑ (h)
γϑ (0)
..
.

γϑ ((r − 1)h) γϑ ((r − 2)h) · · ·

=: Γ (r−1)(ϑ )−1 γ (r−1) (ϑ ),

2

σ (ϑ )

=

···
···

∗

T (r−1)

γϑ (0) − π (ϑ ) γ

−1


γϑ ((r − 1)h)
γϑ (h)

γϑ ((r − 2)h)  . 
  .. 
..

.
γϑ (rh)
γϑ (0)

(ϑ ).

(2.5a)
(2.5b)

Remark 2.3. Um (ϑ ) can be interpreted as the error of the best linear predictor of Yϑ (mh) in terms of
Y(m−1)h (ϑ ), . . . ,Y(m−r)h (ϑ ). Per construction, however, the sequence (Um (ϑ ))m∈Z is not an uncorrelated sequence, Um (ϑ ) is only uncorrelated with Y(m−1)h (ϑ ), . . . ,Y(m−r)h (ϑ ).
Definition 2.4. Let Π ⊆ Rr+1 be the parameter space containing all possible parameter vectors of
stationary AR(r) processes. The map π : Θ → Π with ϑ 7→ π (ϑ ) and π (ϑ ) as given in Theorem 2.2
is called the link function or binding function.
Lemma 2.5. Let r ≥ 2p − 1. Then, π (ϑ ) is injective and three times continuously differentiable.
Proof. We make use of the fact that the discretely observed CARMA(p, q(ϑ )) process (Ymh (ϑ ))m∈Z
admits a representation as a stationary ARMA(p, p − 1) process with weak white noise of the form

φ (B)Ymh (ϑ ) = θ (B)εm (ϑ ),

(2.6)

p
(1 − ehλi z) (the λi being the eigenvalues of Aϑ ), θ (z) is a monic, Schur–stable
where φ (z) = ∏i=1
polynomial and (εm (ϑ ))m∈Z is a weak white noise (see (Brockwell and Lindner, 2009, Lemma 2.1)).
We can now decompose the map π : Θ → Π into three separate maps for which we define the
following spaces:

M := {(φ1 , . . . , φ p , θ1 , . . . , θ p−1 , σ ) ∈ R2p : The coefficients define a weak ARMA(p, p − 1)
model as in (1.1) for which φ (z) and θ (z) have no common zeros} ⊆ R2p ,

G := {γ = (γ0 , . . . , γr ) ∈ Rr+1 : The coefficients define the autocovariances up to order
r of a stationary stochastic process where Γ (r−1) is non-singular} ⊆ Rr+1 ,

6

Π := {(π1 , . . . , πr , σ ) ∈ Rr × (0, ∞) : (π1 , . . . , πr ) are the coefficients of a stationary
AR(r) process and σ 2 is the variance of the noise} ⊆ Rr+1 ,

where Γ (r−1) is defined as Γ (r−1) (ϑ ) in (2.5a). Denote by π1 : Θ → M the map which maps the
parameters of a CARMA process to the coefficients of the weak ARMA(p, p − 1) representation of its
sampled version as in (2.6). Denote by π2 : M → G the map which maps the parameters of a weak
ARMA(p, p − 1) process to its autocovariances of lags 0, . . . , r. Lastly, denote by π3 : G → Π the map
which maps a vector of autocovariances (γ0 . . . , γr ) to the parameters of the auxiliary AR(r) model.
Then we have that π = π3 ◦ π2 ◦ π1 . We will show that πi is injective for i = 1, 2, 3 and receive from
this the injectivity of π . The three-times continuous-differentiability of the map π follows from the
representation (2.5) and the three-times continuous-differentiability of the autocovariance function γϑ .
Step 1: π1 is injective.
Due to Assumption A and (Schlemm and Stelzer, 2012, Theorem 3.13) the family of sampled processes {(Ymh (ϑ )m∈Z : ϑ ∈ Θ )} is identifiable from the spectral densities and hence, for any ϑ 6= ϑ ′ ∈ Θ
the parameters of the weak ARMA process in (2.6) differ.
Step 2: π2 is injective if r ≥ 2p − 1.
The reason is that by the method of (Brockwell and Davis, 1991, p. 93), the autocovariances of an
ARMA(p, p − 1) process are completely determined as the solution of a difference equation with p
boundary conditions which depend on the coefficient vector (φ1 , . . . , φ p , θ1 , . . . , θ p−1 , σ ). If r ≥ 2p − 1
the number of equations r is greater than or equal to the number of variables 2p − 1 which results in
the injectivity of π2 (see also (de Luna and Genton, 2001, Section 4.1)). To be more precise, let
e) ∈ M .
θ = (φ1 , . . . , φ p , θ1 , . . . , θ p−1 , σ ) ∈ M and θe = (φe1 , . . . , φep , θe1 , . . . , θep−1 , σ
(p)
(p+1)×(p+1)
Case 1. (φ1 , . . . , φ p ) 6= (φe1 , . . . , φep ). Define Γ (θ ) ∈ R
similarly to Γ (p−1) (ϑ ) in (2.5a).
Due to (Brockwell and Davis, 1991, (3.3.9))
(−φ p
(−φep

...
...

− φ1
− φe1

1)Γ (p) (θ ) = (0
1)Γ (p) (θe) = (0

0
0

···
···

0),
0).

But since the vectors (−φ p . . . − φ1 1) and (−φep . . . − φe1 1) are linear independent this is
only possible if Γ (p) (θ ) 6= Γ (p) (θe) which implies π2 (θ ) 6= π2 (θe).
Case 2. (φ1 , . . . , φ p ) = (φe1 , . . . , φep ). Assume that π2 (θ ) = π2 (θe). But then due to (Brockwell and Davis,
1991, (3.3.9)), (γθ (k))k∈N0 = (γθe (k))k∈N0 and hence, θ = θe.
Step 3: π3 is injective.
We can also rewrite the linear equations (2.5) as a linear system with r + 1 equations and the r + 1
unknown variables γ0 , . . . , γr which gives the injectivity of π3 .
Finally, due to Lemma 2.5 we suppose throughout the paper:
Assumption B. Let r ≥ 2p − 1.

3 Indirect estimation
For fixed r, denote by πbn an estimator of π (ϑ0 ) that is calculated from the observations Y n =
(Yh , . . . ,Ynh ). If we were able to analytically invert the link function π and calculate π −1 (πbn ), then
π −1 (πbn ) we would be an estimator for ϑ0 = π −1 (π (ϑ0 )). However, this is not possible in general
since no analytic representation of π −1 exists. To overcome this problem, we perform a second estimation, which is based on simulations and constitutes the other building block of indirect estimation.

7

We fix a number s ∈ N and simulate a sample path of length sn of a Lévy process (LtS )t∈R with
ELS1 = 0 and E(LS1 )2 = σL2 . Then, for a fixed parameter ϑ ∈ Θ we generate a sample path of the
associated CARMA process (YtS (ϑ ))t∈R using the simulated path (LtS )t∈R . This gives us a vector of
S (ϑ )) of length sn. From this observation Y sn (ϑ )
“pseudo–observations” YSsn (ϑ ) = (YhS (ϑ ), . . . ,Ysnh
S
S
we estimate again π (ϑ ) by an estimator πbsn (ϑ ). The idea is now to choose that value of ϑ as estimator for ϑ0 which minimizes a suitable distance between πbn and πbnS (ϑ ). The formal definition is as
follows.

S (ϑ ) be an
Definition 3.1. Let πbn be an estimator for π (ϑ0 ) calculated from the data Y n , let πbsn
estimator for π (ϑ ) calculated from the pseudo–observations YSsn (ϑ ) = (YhS (ϑ ), . . . ,YϑS (snh)) and let
Ω ∈ RN(Θ )×N(Θ ) be a symmetric positive definite weighting matrix. The function LInd : Θ → [0, ∞) is
defined as
S
S
LInd (ϑ , Y n ) := [πbn − πbsn
(ϑ )]T Ω [πbn − πbsn
(ϑ )].

Then, the indirect estimator for ϑ0 is

ϑbnInd = arg min LInd (ϑ , Y n ).
ϑ ∈Θ

We are able to present general conditions under which this indirect estimator is consistent and
asymptotically normally distributed.
Theorem 3.2.
(a) Suppose that the following assumptions are satisfied:
P

(C.1) πbn → π (ϑ0 ) as n → ∞.

P

(C.2) supϑ ∈Θ kπbnS (ϑ ) − π (ϑ )k → 0 as n → ∞.
Then

P

sup |LInd (ϑ , Y n ) − QInd (ϑ )| → 0

and

ϑ ∈Θ

P
ϑbnInd → ϑ0 .

If we replace in (C.1) and (C.2) convergence in probability by almost sure convergence then we
can replace in the statement convergence in probability by almost sure convergence as well.
(b) Assume additionally to (C.1) and (C.2):
√
D
(C.3) n(πbnS (ϑ ) − π (ϑ )) −→ N (0, ΞS (ϑ )) as n → ∞ for any ϑ ∈ Θ .
√
D
(C.4) n(πbn − π (ϑ0 )) −→ N (0, ΞD (ϑ0 )) as n → ∞.
P

(C.5) For any sequence (ϑ n )n∈N with ϑ n → ϑ0 as n → ∞ the asymptotic behaviors
P

∇ϑ πbnS (ϑ n ) → ∇ϑ π (ϑ0 ),
∇2ϑ πbnS (ϑ n ) = OP (1),

hold as n → ∞ and ∇ϑ π (ϑ0 ) has full column rank N(Θ ).

Then, as n → ∞,
where

√
D
n(ϑbnInd − ϑ0 ) −→ N (0, ΞInd (ϑ0 )),

ΞInd (ϑ0 ) = JInd (ϑ0 )−1 IInd (ϑ0 )JInd (ϑ0 )−1

8

with
JInd (ϑ0 ) = [∇ϑ π (ϑ0 )]T Ω [∇ϑ π (ϑ0 )] and


1
IInd (ϑ0 ) = [∇ϑ π (ϑ0 )]T Ω ΞD (ϑ0 ) + ΞS (ϑ0 ) Ω [∇ϑ π (ϑ0 )].
s
Gouriéroux et al. (1993) develop for a dynamic model as well the consistency and the asymptotic
normality of the indirect estimator but under different assumptions mainly based on LInd (ϑ , Y n ) (see
as well Smith (1993)). These results are again summarized in Gouriéroux and Monfort (1997). In
the context of indirect estimation of ARMA models, (de Luna and Genton, 2001, p.22) mention the
asymptotic normality of their indirect estimator but without stating any regularity conditions and only
refereeing to (Gouriéroux and Monfort, 1997, Proposition 4.2).
Proof.
(a) We first start by proving the consistency. To this end, we define the map
QInd : Θ → [0, ∞)

as

ϑ 7→ [π (ϑ ) − π (ϑ0 )]T Ω [π (ϑ ) − π (ϑ0 )].

(3.1)

With this, we obtain
sup |LInd (ϑ , Y n ) − QInd (ϑ )|

ϑ ∈Θ

=
≤

S
S
(ϑ )]T Ω [πbn − πbsn
(ϑ )] − [π (ϑ ) − π (ϑ0 )]T Ω [π (ϑ ) − π (ϑ0 )]|
sup |[πbn − πbsn
ϑ ∈Θ
S
|πbnT Ω πbn − π (ϑ0 )T Ω π (ϑ0 )| + sup |πbsn
(ϑ )T Ω πbn − π (ϑ )T Ω π (ϑ0 )|
ϑ ∈Θ
S
S
S
(ϑ ) − π (ϑ0 )T Ω π (ϑ )| + sup |πbsn
(ϑ )T Ω πbsn
(ϑ ) − π (ϑ )T Ω π (ϑ )|.
+ sup |πbnT Ω πbsn
ϑ ∈Θ
ϑ ∈Θ

The four summands on the right–hand side converge in probability to zero as n → ∞. For the first one,
this is a consequence of (C.1). For the remaining three terms, the arguments are similar, so that we
treat only the second one exemplary. We have
S
sup |πbsn
(ϑ )T Ω πbn − π (ϑ )T Ω π (ϑ0 )|

ϑ ∈Θ

≤

S
sup |πbsn
(ϑ )T Ω πbn − π (ϑ )T Ω πbn | + sup |π (ϑ )T Ω πbn − π (ϑ )T Ω π (ϑ0 )|

ϑ ∈Θ

≤ kΩ k sup

ϑ ∈Θ

ϑ ∈Θ

S
kπbsn
(ϑ ) − π (ϑ )kkπbn k + kΩ k

P

sup kπ (ϑ )kkπbn − π (ϑ0 )k → 0.

ϑ ∈Θ

Here, we used the fact that supϑ ∈Θ kπ (ϑ )k is finite due to the continuity of the map π and the compactness of Θ as well as both (C.1) and (C.2). Therefore, the function LInd (ϑ , Y n ) converges uniformly
in ϑ in probability to the limiting function QInd (ϑ ). Per construction, ϑbnInd minimizes LInd (ϑ , Y n )
and QInd (ϑ ) has a unique minimum at ϑ = ϑ0 . Therefore, weak consistency of ϑbnInd follows by arguing as in the proof of (Schlemm and Stelzer, 2012, Theorem 2.4); although in their proof convergence
in probability is replaced by almost sure convergence, this doesn’t matter because we can use the
subsequence criterion which says that a sequence converges in probability iff any subsequence has a
further subsequence which converges almost surely.
The proof of strong consistency goes similarly by replacing convergence in probability by almost
sure convergence.

9

(b) For the asymptotic normality, note that
√
√
√
S
S
(ϑ0 )) = n(πbn − π (ϑ0 )) + n(π (ϑ0 ) − πbsn
(ϑ0 )).
n(πbn − πbsn

Since both estimators are independent from each other, we obtain with (C.3) and (C.4) that


√
1
D
S
(ϑ0 )) −→ N 0, ΞD (ϑ0 ) + ΞS (ϑ0 ) .
n(πbn − πbsn
s

(3.2)

Moreover,

0 = ∇ϑ LInd (ϑ , Y n )

ϑ =ϑbnInd

S bInd T
S bInd
= 2[∇ϑ πbsn
(ϑn )] Ω [πbsn
(ϑn ) − πbn ].

We now use a Taylor expansion of order 1 around the true value ϑ0 to obtain
√
n
0 =
n∇ϑ LInd (ϑbInd
, Y n)
√
√
=
n∇ϑ LInd (ϑ0 , Y n ) + n∇2ϑ LInd (ϑ n , Y n )(ϑbnInd − ϑ0 )
√
√ S
S
S
S
= 2[∇ϑ πbsn
(ϑ0 )]T Ω n[πbsn
(ϑ0 ) − πbn ] + 2[∇2ϑ πbsn
(ϑ n )]T Ω [πbsn
(ϑ n ) − πbn ] n(ϑbnInd − ϑ0 )
√
+2[∇ϑ πbS (ϑ n )]T Ω [∇ϑ πbS (ϑ n )] n(ϑbInd − ϑ0 ).
sn

sn

n

P

Here, ϑ n is such that kϑ n − ϑ0 k ≤ kϑbnInd − ϑ0 k and hence, ϑ n → ϑ0 as n → ∞. Moreover,
n

P

S
S
S
S
[∇2ϑ πbsn
(ϑ )]T Ω [πbsn
(ϑ n ) − πbn ] + [∇ϑ πbsn
(ϑ n )]T Ω [∇ϑ πbsn
(ϑ n )] → [∇ϑ π (ϑ0 )]T Ω [∇ϑ π (ϑ0 )] (3.3)

due to (C.1), (C.2), (C.5) and the continuity of π (ϑ ). Furthermore, the right-hand side is non-singular
since ∇ϑ π (ϑ0 ) has full column rank and Ω is non-singular. Finally, we write
√
n(ϑbnInd − ϑ0 )

−1
√
n
S
S
S
S
S
S
(ϑ )]T Ω [πbsn
(ϑ n ) − πbn ] + [∇ϑ πbsn
(ϑ n )]T Ω [∇ϑ πbsn
(ϑ n )]
(ϑ0 )]T Ω n(πbsn
(ϑ0 ) − πbn )
= [∇ϑ2 πbsn
[∇ϑ πbsn
and use (3.2), (3.3) and (C.5) to obtain as n → ∞,

−1
√
D
n(ϑbnInd − ϑ0 ) −→ [∇ϑ π (ϑ0 )]T Ω [∇ϑ π (ϑ0 )]
[∇ϑ π (ϑ0 )]T Ω · N




1
0, ΞD (ϑ0 ) + ΞS (ϑ0 ) .
s

This completes the proof.
Remark 3.3.

(a) The asymptotic covariance matrix can be written as


1
ΞInd (ϑ0 ) = H (ϑ0 ) ΞD (ϑ0 ) + ΞS (ϑ0 ) H (ϑ0 )T ,
s
where H (ϑ0 ) = [∇ϑ π (ϑ0 )T Ω ∇ϑ π (ϑ0 )]−1 [∇ϑ π (ϑ0 )]T Ω . This is the analog form of
(de Luna and Genton, 2001, Eq. (4.4)).
(b) Note that the asymptotic results hold for any r ≥ 2p − 1. But increasing the auxiliary AR order
does not necessarily yield better results. On the other hand, increasing s increases the efficiency.
For s → ∞ we receive ΞInd (ϑ0 ) → H (ϑ0 )ΞD (ϑ0 )H (ϑ0 )T . The best efficiency is received for
−1
s→∞ 
Ω = [ΞD (ϑ0 )]−1 in which case ΞInd (ϑ0 ) → ∇ϑ π (ϑ0 )T ΞD (ϑ0 )−1 ∇ϑ π (ϑ0 ) .
10

4 Estimating the auxiliary AR(r) parameters of a CARMA
process with outliers
In order to apply the indirect estimator to a discretely sampled CARMA process we need strongly consistent and asymptotically normally distributed estimators for the parameters of the auxiliary AR(r)
representation. In this section we will study generalized M- (GM-) estimators. The GM-estimator
will be applied to a CARMA process afflicted by outliers because we want to study some robustness
properties of our estimator as well. Outliers can be thought as typical observations that do not arise because of the model structure but due to some external influence, e.g., measurement errors. Therefore,
a whole sample of observations which contains outliers does not come from the true model anymore
but it is still close to it as long as the total number of outliers is not overwhelmingly large.
Definition 4.1. Let g : [0, 1] → [0, 1] be a function that satisfies g(γ )− γ = o(γ ) for γ → 0. Let (Vm )m∈Z
be a stochastic process taking only the values 0 and 1 with
P(Vm = 1) = g(γ )
γ

and let (Zm )m∈Z be a real-valued stochastic process. The disturbed process (Ymh (ϑ ))m∈Z is defined as
γ

Ymh (ϑ ) = (1 −Vm )Ymh (ϑ ) +Vm Zm .

(4.1)

γ

The disturbed process (Ymh (ϑ ))m∈Z is in general not a sampled CARMA process anymore.
Remark 4.2.
(a) The interpretation of this model is that at each point m ∈ Z an outlier is observed with probability
g(γ ) while the true value Ymh (ϑ ) is observed with probability 1 − g(γ ). The model has the
advantage that one can obtain both additive and replacement outliers by choosing the processes
(Zm )m∈Z and (Vm )m∈Z adequately. Specifically, to model replacement outliers, one assumes
that (Zm )m∈Z , (Vm )m∈Z and (Ymh (ϑ ))m∈Z are jointly independent. Then, if the realization of Vm
is equal to 1, the value Ymh (ϑ ) will be replaced by the realization of Zm justifying the use of
the name replacement outliers. On the other hand, modeling additive outliers can be achieved
by taking Zm = Ymh (ϑ ) + Wm for some process (Wm )m∈Z and assuming that (Ymh (ϑ ))m∈Z is
γ
independent from (Vm )m∈Z . Then we have Ymh (ϑ ) = Ymh (ϑ ) + VmWm such that the realization
of Wm is added to the realization of Ymh (ϑ ) if Vm is 1.
(b) Another advantage of this general outlier model is that one can easily model the temporal structure of outliers. On the one hand, if (Vm )m∈Z is chosen as an i.i.d. sequence with P(Vm = 1) = γ ,
then outliers typically appear isolated, i.e., between two outliers there is usually a period of time
where no outliers are present. On the other hand, one can also model patchy outliers by letting
(Bm )m∈Z be an i.i.d. process of Bernoulli random variables with success probability ε and
setting Vm = max(Bm−l , . . . , Bm ) for a fixed l ∈ N. Then as ε → 0,
P(Vm = 1) = 1 − (1 − ε )l = l ε + o(ε ),
which results in γ = l ε . For ε sufficiently small, outliers then appear in a block of size l.
Assumption D.
(D.1) The processes (Vm )m∈Z and (Zm )m∈Z are strictly stationary with E|V1 | < ∞ and E|Z1 | < ∞ for
some δ > 0.

11

(D.2) Either we have the replacement model where the processes (Ymh (ϑ ))m∈Z , (Vm )m∈Z and (Zm )m∈Z
are jointly independent, and (Vm )m∈Z and (Zm )m∈Z are exponentially strong mixing, i.e., αV (m) ≤
Cρ m and αZ (m) ≤ Cρ m for some C > 0, ρ ∈ (0, 1) and every m ∈ N. Or we have the additive
model with Zm = Ymh (ϑ ) + Wm where the processes (Ymh (ϑ ))m∈Z , (Vm )m∈Z and (Wm )m∈Z are
jointly independent, and (Vm )m∈Z and (Wm )m∈Z are exponentially strong mixing.
(D.3) For all a ∈ R, π ∈ Rr with |a| + kπ k > 0:
γ

γ

γ

P(aY(r+1)h (ϑ ) + π1Yrh (ϑ ) + . . . + πrYh (ϑ ) = 0) = 0.
We largely follow the ideas of Bustos (1982) for the GM-estimation of AR(r) parameters, however our
model and our assumptions are slightly different. Assumption D corresponds to (Bustos, 1982, Assumption (M2),(M4),(M5)). The main difference is that the sampled CARMA process (Ymh )m∈Z is in
Bustos (1982) an infinite-order moving average process whose noise is Φ –mixing which is in general
not satisfied for a sampled CARMA process. However, we already know from (Marquardt and Stelzer,
2007, Proposition 3.34) that a CARMA process is exponentially strong mixing which is weaker than
Φ -mixing. Therefore, we assume that (Vm )m∈Z , (Zm )m∈Z and (Wm )m∈Z are exponentially strong mixing instead of Φ -mixing as in Bustos (1982).
In the following we define GM-estimators. Let two functions φ : Rr × R → R and χ : R → R be
γ
γ
γ
given. Moreover, assume that we have observations Y n,γ (ϑ ) = (Yh (ϑ ),Y2h (ϑ ), . . . ,Ynh (ϑ )) from the
disturbed process in (4.1). The parameter

π GM (ϑ γ ) = (π1GM (ϑ γ ), . . . , πrGM (ϑ γ ), σ GM (ϑ γ ))
is defined as the solution of the equations
 γ
  γ


Yh (ϑ )
Yh (ϑ )
γ
γ
γ
ϑ
)
−
π
ϑ
)
−
.
.
.
−
π
ϑ
)
Y
(
Y
(
Y
(
1 rh
r h
  .. 
 
 (r+1)h
E φ  ...  ,
  .  = 0,
σ
γ
γ
Yrh (ϑ )
Yrh (ϑ )
 
!2  
γ
γ
γ
Y(r+1)h (ϑ ) − π1Yrh (ϑ ) − . . . − πrYh (ϑ )
 = 0
E χ 
σ

(4.2a)

(4.2b)

for (π1 , . . . , πr , σ ) ∈ Rr × (0, ∞). The idea is again that these are the parameters of the auxiliary AR
γ
representation of (Ymh (ϑ ))m∈Z . Note that π GM (ϑ γ ) depends on the processes (Vm )m∈Z and (Zm )m∈Z
as well. We choose not to indicate this in the notation to make the exposition more readable. Now, the
GM (ϑ γ ), . . . , π
GM (ϑ γ ), σ
bn,r
bnGM (ϑ γ )) based on φ and χ is defined to satisfy
GM-estimator πbnGM (ϑ γ ) = (πbn,1



γ
γ
Ykh (ϑ )
Ykh (ϑ )

 γ



GM (ϑ γ )Y γ
GM (ϑ γ )Y γ (ϑ ) 
πn,r
1
(ϑ )−...−b
 Y(k+r)h (ϑ )−bπn,1


..
..
kh
(k+r−1)h
,



φ




.
.
bnGM (ϑ γ )
σ
n − r k=1  γ

 γ

Y(k+r−1)h (ϑ )
Y(k+r−1)h (ϑ )

! 
γ
γ)−π
GM (ϑ γ )Y γ
GM (ϑ γ )Y γ (ϑ ) 2
n−r
b
b
ϑ
ϑ
)
−
.
.
.
−
π
Y
(
(
n,r
1
n,1
kh
(k+r)h
(k+r−1)h

χ
n − r k=1
σbnGM (ϑ γ )
n−r

∑



∑

= 0, (4.3a)

= 0. (4.3b)

Throughout the paper we assume that there exists a solution of (4.3) although this is not always the

12

case in practice.
Example 4.3.
(a) There are two main classes of GM-estimators, the so–called Mallows estimators and the Hampel–
Krasker–Welsch estimators. More information on them can be found in Bustos (1982);
Denby and Martin (1979); Martin (1980); Martin and Yohai (1986). In the literature, this kind
of estimators sometimes appear under the name BIF (for bounded influence) estimators. The
class of Mallows estimators are defined as φ (y, u) = w(y)ψ (u), where w is a strictly positive
weight function and ψ is a suitably chosen robustifying function. The Hampel–Krasker–Welsch
estimators are of the form
ψ (w(y)u)
φ (y, u) =
,
w(y)
where w is a weight function and ψ is again a suitably chosen bounded function.
(b) Typical choices for ψ are the Huber ψk –functions (cf. (Maronna et al., 2006, Eq. (2.28))).
Those functions are defined as ψk (u) = sign(u) min{|u|, k} for a constant k > 0. A possibility
for w is, e.g., w(y) = ψk (|y|)/|y| for a Huber function ψk . Another choice for ψ is the so-called
Tukey bisquare (or biweight) function which is given by


u2
ψ (u) = u 1 − 2
k

2

1{|u|≤k} ,

where k is a tuning constant.
(c) For the function χ , a possibility is χ (x2 ) = ψ 2 (x) − EZ [ψ 2 (Z)] with the same ψ function as in
the definition of φ . The random variable Z is suitably distributed.
In order to develop an asymptotic theory and to obtain a robust estimator it is necessary to impose
assumptions on φ and χ which we will do next analogous to (Bustos, 1982, (E1) - (E6)):
Assumption E.

Suppose φ : Rr × R → R and χ : R → R satisfy the following assumptions:

(E.1) For each y ∈ Rr , the map u 7→ φ (y, u) is odd, uniformly continuous and φ (y, u) ≥ 0 for u ≥ 0.
(E.2) (y, u) 7→ φ (y, u)y is bounded and there exists a c > 0 such that
|φ (y, u)y − φ (z, u)z| ≤ cky − zk
(E.3) The map u 7→

φ (y,u)
u

for all u ∈ R.

is non-increasing for y ∈ Rr and there exists a u0 ∈ R such that

φ (y,u0 )
u0

> 0.

(E.4) φ (y, u) is differentiable with respect to u and the map u 7→ ∂ φ∂(y,u)
is continuous, while (y, u) 7→
u
∂ φ (y,u)
∂ u y is bounded.
 
   γ
 

 γ
Yh (ϑ )
Yh (ϑ )





 ∂  ..    ..  
(E.5) E sup u  φ  .  , u  .   < ∞.

∂u
u∈R 


γ
γ
Yrh (ϑ )
Yrh (ϑ )
(E.6) χ is bounded and increasing on {x : −a ≤ χ (x) < b} where b = supx∈R χ (x) and a = −χ (0).
Furthermore, χ is differentiable and x 7→ xχ ′ (x2 ) is continuous and bounded. Lastly, χ (u20 ) > 0.

13

In the remaining of this section we always assume that Assumption D and E are satisfied.
Remark 4.4. As pointed out in (Bustos, 1982, p. 497) one can deduce from (Maronna and Yohai,
1981, Theorem 2.1) that there exists a solution π GM (ϑ γ ) ∈ Rr × (0, ∞) of equation (4.2) if
Assumption E holds. Moreover, there exists a compact set K ⊂ Rr × (0, ∞) with π GM (ϑ γ ) ∈ K and
for any π ∈ K c equation (4.2) does not hold (see (Bustos, 1982, p. 500)).
In general it is not easy to verify that π GM (ϑ γ ) is unique. Additionally, one would like to have that
π GM (ϑ 0 ) = π (ϑ ) are the parameters of the auxiliary AR(r) model in the case that the GM-estimator is
applied to realizations of an uncontaminated sampled CARMA process (Ymh (ϑ ))m∈Z . The following
proposition gives a sufficient condition.
Proposition 4.5. Suppose that Ur+1 (ϑ ) as defined in equation (2.3) satisfies
D

(Ur+1 (ϑ ),Yrh (ϑ ), . . . ,Yh (ϑ )) = (−Ur+1 (ϑ ),Yrh (ϑ ), . . . ,Yh (ϑ )).

(4.4)

Assume further that the function u 7→ φ (y, u) is nondecreasing and strictly increasing for |u| ≤ u0 ,
where u0 satisfies Assumptions (E.3) and (E.6), and the function χ is chosen in such a way that
"
 !#

U1 (ϑ ) 2
E χ
= 0.
(4.5)
σ (ϑ )
γ

Finally, assume that γ = 0 so that (Ymh (ϑ ))m∈Z = (Ymh (ϑ ))m∈Z . Then the auxiliary parameter π (ϑ )
as defined in Theorem 2.2 is the unique solution of (4.2), i.e., π GM (ϑ 0 ) = π (ϑ ).
Proof. Using similar arguments as in (Maronna and Yohai, 1981, Lemma 2.1) (limx→0 χ (x) < 0,
limx→∞ χ (x) = ∞, the continuity and boundedness of χ and the Intermediate Value Theorem) we
can show that for each fixed (π1 , . . . , πr ) ∈ Rr there exists a unique solution σ of the equation
"

 !#
Y(r+1)h (ϑ ) − π1Yrh (ϑ ) − . . . − πrYh (ϑ ) 2
= 0.
E χ
σ
By assumption (4.5), the function χ is chosen in such a way that for (π1 (ϑ ), . . . , πr (ϑ )) this unique
solution is σ (ϑ ). Therefore, we have that π (ϑ ) is a solution of (4.2b). Next, we show that π (ϑ ) is a
solution of (4.2a) as well. Since the function φ (y, u) is odd in u by Assumption (E.1), it holds that


 






Yh (ϑ )
Yh (ϑ )
Yh (ϑ )
Yh (ϑ )
 .  Ur+1 (ϑ )   .. 

 
 Ur+1 (ϑ )   .. 
E φ  ...  ,
  .  = E −φ  ..  , −
  . 
σ (ϑ )
σ (ϑ )
Yrh (ϑ )
Yrh (ϑ )
Yrh (ϑ )
Yrh (ϑ )


 

Yh (ϑ )
Yh (ϑ )
  ..  Ur+1 (ϑ )   .. 
(4.6)
= −E φ  .  ,
  .  ,
σ (ϑ )
Yrh (ϑ )
Yrh (ϑ )
where the last equality follows from (4.4). From this equation we can conclude that
 



Yh (ϑ )
Yh (ϑ )
 
 Ur+1 (ϑ )   .. 
E φ  ...  ,
  .  = 0,
σ (ϑ )
Yrh (ϑ )
Yrh (ϑ )
14

and therefore, π (ϑ ) is a solution of equation (4.2a).
Next, we show similarly to (Maronna and Yohai, 1981, Theorem 2.2(a)) for regression models that
π (ϑ ) is the unique solution. Assume that another solution π ′ = (π1′ , . . . , πr′ , σ ′ ) of (4.2) exists. But
then (π1′ , . . . , πr′ ) 6= (π1 , . . . , πr ). Note that the arguments in the derivation of (4.6) still hold if we
replace σ (ϑ ) in the denominator of the second argument of φ by σ ′ . Thus, we obtain that
 



Yh (ϑ )
Yh (ϑ )
 
 Ur+1 (ϑ )   .. 
E φ  ...  ,
  .  = 0,
σ′
Yrh (ϑ )
Yrh (ϑ )
and therefore
 




Yh (ϑ )
Yh (ϑ )
′
′
 
 Y(r+1)h (ϑ ) − π1Yrh (ϑ ) − . . . − πrYh (ϑ )   .. 
E φ  ...  ,
 . 
σ′
Yrh (ϑ )
Yrh (ϑ )

 


Yh (ϑ )
Yh (ϑ )
 ..  Ur+1 (ϑ )   .. 
−φ  .  ,
  .  = 0.
σ′
Yrh (ϑ )
Yrh (ϑ )

(4.7)

Since P((Yh (ϑ ), . . . ,Yrh (ϑ )) = (0, . . . , 0)) = 0 and P(Y(r+1)h (ϑ )− π1′ Yrh (ϑ )−. . .− πr′Yh (ϑ ) = Ur+1 (ϑ )) =
0 due to (D.3) for γ = 0 and u 7→ φ (y, u) is strictly increasing on the interval (−u0 , u0 ) for every y ∈ Rr
we have that
Y(r+1)h (ϑ ) − π1′ Yrh (ϑ ) − . . . − πr′Yh (ϑ )
≥ u0 P-a.s.
(4.8)
σ′
because otherwise (4.7) cannot hold. Now, π ′ is by assumption also a solution of (4.2b) and hence,
we have due to (4.8) and (E.6)
"

 !#
Y(r+1)h (ϑ ) − π1′ Yrh (ϑ ) − . . . − πr′Yh (ϑ ) 2
0=E χ
≥ χ (u20 ) > 0
σ′
which is a contradiction.
Remark 4.6.
(a) Assumption (4.4) holds if the distribution of Ur+1 (ϑ ) is symmetric and Ur+1 (ϑ ) is independent
of (Yrh (ϑ ), . . . ,Yh (ϑ )). This again is satisfied if (Lt )t∈R is a Brownian motion.
(b) The monotonicity assumption on φ is valid, e.g., for both the Mallows and Hampel–Krasker–
Welsch estimators when the function ψ is chosen as a Huber ψk –function with u0 = k.
(c) The assumption on χ is fulfilled, e.g., if χ is chosen as in Example 4.3(c) with
p
D
Z = U1 (ϑ )/ Var(U1 (ϑ )). In the case that the driving Lévy process is a Brownian motion this
means that Z ∼ N (0, 1).
Theorem 4.7. Suppose that there exists a unique solution π GM (ϑ γ ) of
n→∞
πbnGM (ϑ γ ) → π GM (ϑ γ ) P-a.s.

Proof. The proof goes in the same vein as the proof of (Bustos, 1982, Theorem 2.1).

15

(4.2).

Then

Next, we would like to deduce the asymptotic normality of the GM-estimator. We need the following lemma which is the analog of (Bustos, 1982, Lemma 3.1) under our different model assumptions.
Lemma 4.8. Define the map Ψ : Rr × Rr × (0, ∞) → Rr+1 as
  
  
y1
y1
  ..  yr+1 −π1 yr −...−πr y1   .. 
φ  .  ,
  . 
σ
.
Ψ (y, π ) = 

yr
yr 



2
χ yr+1 −π1 yσr −...−πr y1
Furthermore,
define the stochastic process
γ
γ
= Ψ (Ykh (ϑ ), . . . ,Y(k+r−1)h (ϑ ), π GM (ϑ γ )). Then
√

Ψ (ϑ γ ) = (Ψk (ϑ γ ))k∈N

as

Ψk (ϑ γ )

1 n−r
D
∑ Ψk (ϑ γ ) → N (0, IGM (ϑ γ )),
n − r k=1

where the (i, j)-th component of IGM (ϑ γ ) is
∞


[IGM (ϑ γ )]i j = E [Ψ1,i (ϑ γ )Ψ1, j (ϑ γ )] + 2 ∑ E Ψ1,i (ϑ γ )Ψ1+k, j (ϑ γ )

(4.9)

k=1

and Ψk,i (ϑ γ ) denotes the i–th component of Ψk (ϑ γ ), i = 1, . . . , r + 1. Especially, each [IGM (ϑ γ )]i j is
finite for i, j ∈ {1, . . . , r + 1}.
Proof. By the Cramer–Wold device, the statement of the Lemma is equivalent to the assertion that
√ 1 xT ∑n−r Ψk (ϑ γ ) converges to a univariate normal distribution with mean 0 and variance
k=1
n−r
xT IGM (ϑ γ )x for every x ∈ Rr+1 . According to (Ibragimov, 1962, Theorem 1.7), this holds if we
can show that
E|xT Ψk (ϑ γ )|2+δ < ∞

(4.10)

and that (xT Ψk (ϑ γ ))k∈N is strongly mixing with mixing coefficients αxT Ψ (ϑ γ ) (m) satisfying
∞

δ /(2+δ )

∑ αx Ψ (ϑ ) (m) < ∞
T

γ

for some δ > 0.

(4.11)

m=1

The same theorem then also states that xT IGM (ϑ γ )x < ∞ from which we then deduce that for i, j ∈
{1, . . . , r + 1} the entry [IGM (ϑ γ )]i j is finite and therefore, IGM (ϑ γ ) is well–defined.
We start to show the existence of the (2 + δ )–th moment of xT Ψk (ϑ γ ) in (4.10). Therefore, note
that
E|xT Ψk (ϑ γ )|2+δ ≤ Ckxk2+δ

r+1

∑ EkΨk,i (ϑ γ )k2+δ < ∞,

(4.12)

i=1

where the last inequality holds since Ψk,i (ϑ γ ) is bounded by (E.2) and (E.6).
γ
Finally, the process (Ymh (ϑ ))m∈Z is strongly mixing and the mixing coefficients satisfy the above
condition (4.11) for the following reason. Either we have in the case of replacement outliers that
γ
Ymh (ϑ ) = G(Vm , Zm ,Ymh (ϑ )) for some measurable function G and the three processes (Vm ), (Zm ) and

16

γ

(Ymh (ϑ )) are independent, or in the case of additive outliers we have Ymh (ϑ ) = G(Vm ,Wm ,Ymh (ϑ ))
for some measurable function G and the three processes (Vm ), (Wm ) and (Ymh (ϑ )) are independent.
Hence, by (Bradley, 2007, Theorem 6.6(II)), Assumption (D.2) and (Marquardt and Stelzer, 2007,
Proposition 3.34) we receive

αY γ (ϑ ) (m) ≤ αV (m) + αZ (m) + αY (ϑ ) (m) ≤ Cρ m
respectively, αY γ (ϑ ) (m) ≤ αV (m) + αW (m) + αY (ϑ ) (m) ≤ Cρ m for some C > 0 and ρ ∈ (0, 1). Furγ
γ
thermore, Ψk (ϑ γ ) depends only on the finitely many values Ykh (ϑ ), . . . ,Y(k+r)h (ϑ ) and by (Bradley,
2007, Remark 1.8(b)) this ensures that αΨ (ϑ γ ) (m) ≤ αY γ (ϑ ) (m + r) ≤ Cρ m . Thus, the strong mixing coefficients αxT Ψ (ϑ γ ) (m) of xT Ψ (ϑ γ ) satisfy the summability condition (4.11) and the lemma is
proven.
Lemma 4.9. Let the set K be given as in Remark 4.4 and for π = (π1 , . . . , πr , σ ) ∈ K define
 γ
   γ


Yh (ϑ )
Yh (ϑ )
γ
γ
γ
   ..  Y(r+1)h (ϑ )−π1Yrh (ϑ )−...−πrYh (ϑ )   .. 
E φ  .  ,
  . 
σ


γ
γ
γ
.

QGM (π , ϑ ) = 
Yrh"(ϑ )
Yrh (ϑ )

!#




2
γ
γ
γ
Y(r+1)h (ϑ )−π1Yrh (ϑ )−...−πr Yh (ϑ )


E χ
σ

(4.13)

P

Then ∇π QGM (π , ϑ γ ) exists. Moreover, for any sequence (πn )n∈N with πn → π GM (ϑ γ ) as n → ∞ we
have as n → ∞,
P

∇π QGM (πn , ϑ γ ) → ∇π QGM (π GM (ϑ γ ), ϑ γ ).
Proof. Note, first that for i, j = 1, . . . , r,
 γ


Yh (ϑ )
γ
γ
γ
(
Y
(
Y
(
ϑ
)
−
π
ϑ
)
−
.
.
.
−
π
ϑ
)
Y
1 rh
r h
∂  .  (r+1)h
 γ
sup
φ  ..  ,
 Y jh (ϑ )
σ
π ∈K ∂ πi
γ
Yrh (ϑ )

 γ

Yh (ϑ )
γ
γ
γ
γ


(
Y
(
Y
(
Y(r+1−i)h (ϑ )
ϑ
)
−
π
ϑ
)
−
.
.
.
−
π
ϑ
)
Y
1 rh
r h
∂
 γ
 ..  (r+1)h
= sup
φ  .  ,
 Y jh (ϑ )
∂u
σ
σ
π ∈K
γ
Yrh (ϑ )
 γ
 γ
  γ
  γ


Yh (ϑ )
Yh (ϑ )
Yh (ϑ )
Yh (ϑ )


∂


 
 


φ  ...  , u  ...   ...  ≤ C  ... 
≤ sup C
∂u
u∈R
γ
γ
γ
γ
Yrh (ϑ )
Yrh (ϑ )
Yrh (ϑ )
Yrh (ϑ )

due to Assumption (E.4) and the boundedness of 1/σ on the compact set K. By Assumption (D.1)
and (A.3) the expectation on the right-hand side is finite. Similarly,
 γ
 γ


Yh (ϑ )
Yh (ϑ )
γ
γ
γ
∂  .  Y(r+1)h (ϑ ) − π1Yrh (ϑ ) − . . . − πrYh (ϑ )   . 
sup
φ  ..  ,
  .. 
σ
π ∈K ∂ σ
γ
γ
Yrh (ϑ )
Yrh (ϑ )

17



  γ

γ
Yh (ϑ )
Yh (ϑ )
∂

 

≤ C sup u
φ  ...  , u  ...  .
∂
u
u∈R
γ
γ
Yrh (ϑ )
Yrh (ϑ )




The expectation on the right-hand side is finite due to Assumption (E.5). Similar
arguments, us γ
2 !
γ
γ
Y
(
ϑ
)−
π
Y
ϑ
)−...−
π
Y
ϑ
)
(
(
1 rh
r h
(r+1)h
ing Assumption (E.6), also show that ∂∂πi χ
for i = 1, . . . , r and
σ

 !
γ

∂
∂σ χ

γ

γ

Y(r+1)h (ϑ )−π1Yrh (ϑ )−...−πr Yh (ϑ )
σ

2

are uniformly dominated by integrable random variables.

Therefore, by (Billingsley, 1999, Theorem 16.8(ii)) (that is an application of dominated convergence)
∇π QGM (π , ϑ γ ) exists on K and the order of differentiation and expectation can be changed.
Moreover, due (E.4), (E.6) and (Billingsley, 1999, Theorem 16.8(i)) the map π 7→ ∇π QGM (π , ϑ γ )
P

P

is continuous. Hence, if πn → π GM (ϑ γ ) ∈ K then ∇π QGM (πn , ϑ γ ) → ∇π QGM (π GM (ϑ γ ), ϑ γ ).

Lemma 4.10. Let QGM (π , ϑ γ ) be defined as in (4.13) and suppose that ∇π QGM (π , ϑ γ ) is nonP
singular. Furthermore, let IGM(ϑ γ ) be given as in (4.9) and suppose that πbnGM (ϑ γ ) → π GM (ϑ γ )
as n → ∞. Then, as n → ∞,
√

D
n − rQGM (πbnGM (ϑ γ ), ϑ γ ) −→ N (0, IGM (ϑ γ )).

Proof. We use the decomposition
√

1 n−r
1 n−r
[QGM (πbnGM (ϑ γ ), ϑ γ ) + Ψk (ϑ γ )] − √
n − pQGM (πbnGM (ϑ γ ), ϑ γ ) = √
∑
∑ Ψk (ϑ γ ).
n − r k=1
n − r k=1

The first term is of order oP (1) due to (Bustos, 1982, Lemma 3.5) (cf. (Kimmig, 2016, Lemma A.5)
in our setting). The second term converges to N (0, IGM (ϑ γ )) due to Lemma 4.8. Hence, we receive
the statement.
Finally, we obtain the following analog version of (Bustos, 1982, Theorem 2.2) in our setting.
Theorem 4.11. Let QGM (π , ϑ γ ) be defined as in (4.13) and suppose that JGM (ϑ γ ) := ∇π QGM (π , ϑ γ )
P
is non-singular. Furthermore, let IGM (ϑ γ ) be given as in (4.9) and suppose that πbnGM (ϑ γ ) →
π GM (ϑ γ ) as n → ∞. Then, as n → ∞,
where

√
D
n − r(πbnGM (ϑ γ ) − π GM (ϑ γ )) −→ N (0, ΞGM (ϑ γ )),

ΞGM (ϑ γ ) := [JGM (ϑ γ )]−1 IGM (ϑ γ )[JGM (ϑ γ )]−1 .

(4.14)

Proof. Due to (4.2) we have QGM (π (ϑ γ ), ϑ γ ) = 0. Next, a first-order Taylor expansion around
πbnGM (ϑ γ ) gives
√
0 = n − rQGM (π (ϑ γ ), ϑ γ )
√
√
γ
γ
GM
(ϑ γ ) − πbnGM (ϑ γ )),
= n − rQGM (πbnGM (ϑ γ ), ϑ γ ) + n − r∇π QGM (π GM
n (ϑ ), ϑ )(π

γ
GM (ϑ γ ) − π
bnGM (ϑ γ )k. The statement follows then from a combiwhere kπ GM (ϑ γ ) − π GM
n (ϑ )k ≤ kπ
nation of Lemma 4.9 and Lemma 4.10.

18

5 The indirect estimator for the CARMA parameters
In Section 3 we already introduced the indirect estimator and presented in Theorem 3.2 sufficient
criteria for the indirect estimator to be consistent and asymptotically normally distributed. In the
following we want to show that these assumptions are satisfied in the setting of discretely sampled
CARMA processes when we use as estimator πbnS (ϑ ) in the simulation part the least-squares- (LS-)
estimator πbnLS (ϑ ) and for πbn the GM-estimator πbnGM (ϑ0 ).
LS (ϑ ) =
S (ϑ )) the LS-estimator π
bsn
Definition 5.1. Based on the sample YSsn (ϑ ) = (YhS (ϑ ), . . . ,Ysnh
LS (ϑ ), . . . , π
LS (ϑ ), σ
LS (ϑ )) of π (ϑ ) minimizes
bsn,r
bsn
(πbsn,1

LLS (π , YSsn (ϑ )) :=

2
1 sn−r  S
S
S
ϑ
)
−
π
ϑ
)
−
.
.
.
−
π
ϑ
)
Y
(
Y
(
Y
(
1 (k+r−1)h
r kh
∑ (k+r)h
sn − r k=1

(5.1)

LS (ϑ ) is defined as
bsn
in Π ′ := π (Θ ) and σ
2
σbLS,sn
(ϑ ) =

2
1 sn−r  S
S
S
LS
LS
b
b
(
(
ϑ
)
−
π
ϑ
)Y
ϑ
)
−
.
.
.
−
π
ϑ
)Y
ϑ
)
.
(
Y
(
(
∑ (k+r)h
sn,1
sn,r
kh
(k+r−1)h
sn − r k=1

Remark 5.2. The quasi ML-function for the auxiliary AR(r) parameters of the discretely sampled
CARMA process is defined as
!
S
S
S (ϑ ))2
(ϑ ) − π1Y(k+r−1)h
(Y(k+r)h
(ϑ ) − . . . − πrYkh
1 sn−r
sn
2
LMLE (π , YS (ϑ )) =
∑ log(σ ) +
sn − r k=1
σ2
sn
MLE (ϑ ) = arg min
and the quasi ML-estimator as πbsn
π ∈Π ′ LMLE (π , YS (ϑ )). It is well known that for
the estimation of AR(r) parameters the ML-estimator and the LS-estimator are equivalent (this can
be seen by straightforward calculations taking the derivatives of the ML-function LMLE which are
proportional to the derivatives of LLS ).

Theorem 5.3. Let Assumption A, B, D and E hold. Suppose that the unique solution π GM (ϑ00 ) of (4.2)
for (Ymh )m∈Z is π (ϑ0 ), that ∇ϑ π (ϑ0 ) has full column rank N(Θ ) and that JGM (ϑ0 ) is non-singular.
∗
Further, assume that E|LS1 |2N for some N ∗ ∈ N with 2N ∗ > max(N(Θ ), 4 + δ ). If πbnS (ϑ ) = πbnLS (ϑ )
and πbn = πbnGM (ϑ0 ) then the indirect estimator ϑbnInd is weakly consistent and
√

where

D

n(ϑbnInd − ϑ0 ) −→ N (0, ΞInd (ϑ0 )),

ΞInd (ϑ0 ) = JInd (ϑ0 )−1 IInd (ϑ0 )JInd (ϑ0 )−1
with
JInd (ϑ0 ) = [∇ϑ π (ϑ0 )]T Ω [∇ϑ π (ϑ0 )] and


1
IInd (ϑ0 ) = [∇ϑ π (ϑ0 )]T Ω ΞGM (ϑ0 ) + ΞLS (ϑ0 ) Ω [∇ϑ π (ϑ0 )],
s
where the matrix ΞLS (ϑ ) is defined as in (4.14) with φ (y, u) = u and χ (x) = x − 1.

19

We have already proven that (C.1) and (C.4) of Theorem 3.2 are satisfied. To show the remaining
conditions on the LS-estimator πbnLS (ϑ ) we require several auxiliary lemmata. The remaining of this
section is devoted to that.
For the ease of notation we write in the following for the Lévy process (LtS )t∈R shortly (Lt )t∈R and
∗
hence, assume that E|L1 |2N for some 2N ∗ > max(N(Θ ), 4 + δ ); similarly (YtS )t∈R is (Yt )t∈R .
Lemma 5.4. Define for any ϑ ∈ Θ the function fϑ (u) = cTϑ eAϑ u e p 1[0,∞) (u) and

Gϑ (u) = fϑ (u), ∇ϑ fϑ (u), ∇2ϑ fϑ (u) .

Then P-a.s. we have



Ymh (ϑ ), ∇ϑ Ymh (ϑ ), ∇ϑ2 Ymh (ϑ ) m∈Z
which is strongly mixing and ergodic.

=

Z

mh
−∞

Gϑ (mh − u) dLu



m∈Z

The proof is moved to Appendix A.
Proposition 5.5. For j, l ∈ {0, . . . , r} define

γbϑ ,n (l, j) =

1 n−r
∑ Y(k+l)h (ϑ )Y(k+ j)h (ϑ ).
n − r k=1

Then for i, u ∈ {1, . . . , N(Θ )} the following statements hold.
P

(a) supϑ ∈Θ |γbϑ ,n (l, j) − γϑ (l − j)| → 0.
(b) supϑ ∈Θ
(c) supϑ ∈Θ

∂ b
∂ ϑi γϑ ,n (l,

P

j) − ∂∂ϑi γϑ (l − j) → 0.

∂2 b
∂ ϑi ∂ ϑu γϑ ,n (l,

P

j) − ∂ ϑi∂∂ ϑu γϑ (l − j) → 0.

Proof. (a) First, we prove the pointwise convergence of the sample autocovariance function and second, that γbϑ ,n (l, j) is locally Hölder-continuous which results in a stochastic equicontinuity condition.
Then we are able to apply (Pollard, 1990, Theorem 10.2) which gives the uniform convergence.
Step 1. Pointwise convergence. From Lemma 5.4 we already know that (Ymh (ϑ ))m∈Z is a stationary
and ergodic sequence with E|Ymh (ϑ )|2 < ∞ due to E|L1 |2 < ∞. Then Birkoff’s Ergodic Theorem gives
as n → ∞,
P
γbϑ ,n (l, j) → γϑ (l − j).

Step 2. γbϑ ,n (l, j) is locally Hölder-continuous. Let γ ∈ [0, 1 − N(Θ )/(2N ∗ )) and
Uk :=

sup

0<kϑ1 −ϑ2 k<1
ϑ1 ,ϑ2 ∈Θ

|Ykh (ϑ1 ) −Ykh (ϑ2 )|
.
kϑ1 − ϑ2 kγ
d

∗

Since ((Ymh (ϑ ))ϑ ∈Θ )m∈Z is a stationary sequence, Uk = U1 and due to Lemma A.3, EU12N < ∞. In
particular, for any ϑ1 , ϑ2 ∈ Θ with kϑ1 − ϑ2 k < 1 the upper bound
|Ykh (ϑ1 ) −Ykh (ϑ2 )| ≤ Uk kϑ1 − ϑ2 kγ

20

and hence,
|Y(k+l)h (ϑ1 )Y(k+ j)h (ϑ1 ) −Y(k+l)h (ϑ2 )Y(k+ j)h (ϑ2 )|


≤ sup |Y(k+l)h (ϑ )| + sup |Y(k+ j)h (ϑ )| (Uk+l +Uk+ j ) kϑ1 − ϑ2 kγ
ϑ ∈Θ
ϑ ∈Θ
{z
}
|
∗
=:Uk+l,k+
j

hold. Finally,

with

|γbϑ1 ,n (l, j) − γbϑ2 ,n (l, j)| ≤
∗
E(Uk+l,k+
j)

1 n−r ∗
∑ Uk+l,k+ j kϑ1 − ϑ2kγ
n − r k=1

∗
= E(U1+l,1+
j)

for kϑ1 − ϑ2 k < 1

(5.2)


1/2
 
2
2
<∞
≤ 4 E sup Yh (ϑ ) EU1
ϑ ∈Θ

where we used Lemma A.3 to get the finite expectation.
1/γ . Then (5.2) and Markov’s inequality
∗
Step 3. Let ε , η > 0. Take 0 < δ < min{1, ηε /E(U1+l,1+
j )}
give


δγ
∗
< ε.
P  sup |γbϑ1 ,n (l, j) − γbϑ2 ,n (l, j)| > η  ≤ E(U1+l,1+
)
j
η
0<kϑ1 −ϑ2 k<δ
ϑ1 ,ϑ2 ∈Θ

A conclusion of this stochastic equicontinuity condition, the pointwise convergence in Step 1 and
(Pollard, 1990, Theorem 10.2) is the uniform convergence.
The proof of (b,c) goes in the same vein as the proof of (a).
Proposition 5.6.
P
(a) supϑ ∈Θ |πbnLS (ϑ ) − π (ϑ )| → 0.

P

(b) supϑ ∈Θ |∇ϑ πbnLS (ϑ ) − ∇ϑ π (ϑ )| → 0.
P

(c) supϑ ∈Θ |∇ϑ2 πbnLS (ϑ ) − ∇ϑ2 π (ϑ )| → 0.

Proof. Define


γbϑ ,n (r, r − 1)


(r−1)
..
γbn
(ϑ ) = 

.
γbϑ ,n (r, 0)

Then

and



γbϑ ,n (r − 1, r − 1) · · ·

(r−1)
..
Γbn
(ϑ ) = 
.
···
γbϑ ,n (r − 1, 0)


 LS
πbn,1 (ϑ )


(r−1)
(r−1)
πbn∗ (ϑ ) :=  ...  = [Γbn
(ϑ )]−1 γbn
(ϑ ),
LS (ϑ )
πbn,r
(r−1)

2
σLS,n
(ϑ ) =γbϑ ,n (r, r) − [πbn∗ (ϑ )]T γbn

21

(ϑ ).


γbϑ ,n (0, r − 1)

..
.
.

γbϑ ,n (0, 0)

(5.3)

A conclusion of Proposition 5.5(a) and the definition of Γ (r−1)(ϑ ) and γ (r−1) (ϑ ) in (2.5) is that
(r−1)

sup kΓbn

ϑ ∈Θ

P

(ϑ ) − Γ (r−1) (ϑ )k → 0

(r−1)

sup kγbn

and

ϑ ∈Θ

P

(ϑ ) − γ (r−1) (ϑ )k → 0.

(5.4)

Due to the continuity and the positive definiteness of Γ (r−1) (ϑ ) (cf. proof of Proposition 2.2), and the
compactness of Θ we receive supϑ ∈Θ kΓ (r−1) (ϑ )−1 k < ∞. Hence, statement (a) is a consequence of
(5.3)-(5.4) and (2.5).
(b) Note that
h
h
i
i
(r−1)
(r−1)
(r−1)
(r−1)
(r−1)
∂ b∗
−1 ∂ Γ
−1 γb(r−1) (ϑ ) + [Γ
−1 ∂ γ
b
b
b
b
b
π
ϑ
)
=
−[
Γ
ϑ
)]
ϑ
)
[
Γ
ϑ
)]
ϑ
)]
ϑ
)
,
(
(
(
(
(
(
n
n
n
n
∂ ϑi n
∂ ϑi n
∂ ϑi n
h
h
i
i
∂
∗
(r−1) (ϑ )]−1 ∂ Γ (r−1) (ϑ ) [Γ (r−1) (ϑ )]−1 γ (r−1) (ϑ ) + [Γ (r−1) (ϑ )]−1 ∂ γ (r−1) (ϑ ) .
∂ ϑi π (ϑ ) = −[Γ
∂ ϑi
∂ ϑi
(5.5)

Applying Proposition 5.5(b) we receive that
sup
ϑ ∈Θ

∂ b(r−1)
∂ (r−1)
∂ (r−1)
∂ (r−1)
P
P
Γn
Γ
γbn
γ
(ϑ ) −
(ϑ ) −
(ϑ ) → 0 and sup
(ϑ ) → 0. (5.6)
∂ ϑi
∂ ϑi
∂ ϑi
ϑ ∈Θ ∂ ϑi

Then the same arguments as in (a) and (5.4)-(5.6) lead to statement (b).
(c) The proof goes in analog lines as in (a) and (b).
P

Corollary 5.7. Let ϑ n be a sequence in Θ with ϑ n → ϑ0 . Then the following statements hold:
P

(a) πbnLS (ϑ n ) → π (ϑ0 ).
P

(b) ∇ϑ πbnLS (ϑ n ) → ∇ϑ π (ϑ0 ).
P

(c) ∇ϑ2 πbnLS (ϑ n ) → ∇ϑ2 π (ϑ0 ).

Proof. (a) We use the upper bound
kπbnLS (ϑ n ) − π (ϑ0 )k ≤ sup kπbnLS (ϑ ) − π (ϑ )k + kπ (ϑ n ) − π (ϑ0 )k.
ϑ ∈Θ

The first term converges in probability to 0 due Proposition 5.6(a) and the second term because π (ϑ )
P

is continuous (see Lemma 2.5) and ϑ n → ϑ0 . The proof of (b,c) goes on the same way.
Theorem 5.8. For any ϑ ∈ Θ the LS-estimator πbnLS (ϑ ) is strongly consistent and as n → ∞,
√

D

n(πbnLS (ϑ ) − π (ϑ )) −→ N (0, ΞLS (ϑ )) .

Proof. Due to Proposition 5.6(a) we already know that the LS-estimator πbnLS (ϑ ) is consistent. The
asymptotic normality of πbnLS (ϑ ) follows in principle from Theorem 4.11 by interpreting the least
squares estimator as a particular GM-estimator with φ (y, u) = u and χ (x) = x − 1. An assumption of
Theorem 4.11 is that the Jacobian JGM (ϑ ) = ∇π QGM (π (ϑ ), ϑ ) is non–singular. This we can verify

22

by direct calculations for the LS-estimator because

γϑ (0)
γϑ (h) . . . γϑ ((r − 1)h)
 γϑ (h)
γϑ (0) . . . γϑ ((r − 2)h)
1 

..
..
.
..
..
∇π QLS (π , ϑ ) = JLS (ϑ ) = −

.
.
.
σ (ϑ ) 
γϑ ((r − 1)h)
...
...
γϑ (0)
0
0
...
0


0
0

..  .
.

0
2

Hence, JLS (ϑ ) is non–singular if and only if the upper left r × r block is non-singular. However, the upper left block is up to a positive factor the covariance matrix of the random vector
(Yh (ϑ ), . . . ,Yrh (ϑ )) which is non–singular (cf. proof of Proposition 2.2).
Still, we need to be careful because the function φ and χ do not satisfy Assumptions (E.2), (E.4)
and (E.6) with respect to boundedness. However, a close inspection of the proof of Theorem 4.11
reveals that the boundedness is only used at two points. First, in Lemma 4.8 where we deduce the
finiteness of the expectation in (4.12). However, for the LS-estimator


Ψk,i (ϑ ) = Y(k+r)h (ϑ ) − π1 (ϑ0 )Y(k+r−1)h (ϑ ) − . . . − πr (ϑ0 )Ykh (ϑ ) Y(k+i−1)h (ϑ )
for i = 1, . . . , r and

Ψk,r+1 (ϑ ) =



Y(k+r)h (ϑ ) − π1 (ϑ )Y(k+r−1)h (ϑ ) − . . . − πr (ϑ )Ykh (ϑ )
σ (ϑ )

2

− 1.

Therefore, inequality (4.12) follows since the Lévy process (Lt )∈R has finite (4+ δ )–th moment which
then transfers to (Yt (ϑ ))t∈R by (Marquardt and Stelzer, 2007, Proposition 3.30) and subsequently the
2 + δ /2-moment of Ψk,i (ϑ ).
Second, the boundedness assumptions are used in the proof of Lemma 4.9 to deduce the existence
of ∇π QLS (π , ϑ ) and its continuity. But by the above calculations ∇π QLS (π , ϑ ) exists obviously and
is continuous.
Proof of Theorem 5.3. (C.1) and (C.4) follow from Theorem 4.7 and Theorem 4.11. (C.2) is proven
in Proposition 5.6. (C.3) is a consequence of Theorem 5.8. Finally, (C.5) is derived in Corollary 5.7.

6 Robustness properties of the indirect estimator
In this section we study the robustness properties of the indirect estimator for the CARMA parameters
of (Ymh )m∈Z where we assume that πbn = πbnGM (ϑ0 ) is the GM-estimator satisfying Assumptions A, B,
E, (D.3) for (Ymh )m∈Z and that π (ϑ0 ) is the unique solution of (4.2) for (Ymh )m∈Z (a sufficient criterium
for this is given in Proposition 4.5). Moreover, we require similarly to (E.2):
(E2’) (y, u) 7→ φ (y, u) is bounded and there exists c > 0 such that
kφ (y1 , u)y1 − φ (y2 , u)y2 k ≤ cky1 − y2 k/ min(ky1 k, ky2 k),
kφ (y, u1 )y − φ (y, u2 )yk ≤ c|u1 − u2 |/ min(|u1 |, |u2 |).

Under some mild assumptions on ψ (u) and w(y) both the Mallows estimator and the Hampel–
Krasker–Welsch estimator satisfy these conditions (cf. (Boente et al., 1987, p.1305)). For the simula-

23

tion part we take some estimator πbnS (ϑ ), e.g., the LS-estimator, such that as n → ∞,
P

sup kπbnS (ϑ ) − π (ϑ )k → 0

ϑ ∈Θ

holds.

6.1 Resistance and qualitative robustness
Roughly speaking an estimator is robust when small deviations from the nominal model have not
much effect on the estimator. This property is known as qualitative robustness or resistance of the
estimator and was originally introduced in Hampel (1971) for i.i.d. sequences. The same article also
gives a slight extension to the case of data that are generated by permutation–invariant distributions,
introducing the term π –robustness (Hampel, 1971, p.1893). Of course, time series do not satisfy
the assumption of permutation invariance in general. Therefore, there have been various attempts
to generalize the concept of qualitative robustness to the time series setting. (Boente et al., 1987,
Theorem 3.1) prove that their πdn –robustness for time series is equivalent to Hampel’s π –robustness
for i.i.d. random variables and therefore, extends Hampel’s π –robustness. They go ahead and define
the term of resistance as well. The concept of resistance has the intuitive appeal of making a statement
about changes in the values of the estimator when comparing two deterministic samples while πdn –
robustness is only a statement concerning the distribution of the estimator, which is in general not
easily tractable. We apply the definitions in the sense of Boente et al. (1987) to study the robustness
properties for our indirect estimator.
To this end, let y be a (infinite-length) realization of the discretely sampled CARMA process
(Ymh )m∈Z . Formally, we can write that y = (ymh )m∈N ∈ R∞ , where R∞ denotes the infinite cartesian product of R with itself. On this space, equipped with the Borel σ -field B ∞ we denote the
set of all probability measures by P(R∞ ). In the following, we denote for y ∈ R∞ as above by
yn = (yh , y2h , . . . , ynh ) the vector of the first n coordinates. Finally, PY (h) denotes the probability measure of the discrete-sampled CARMA process (Ymh )m∈Z .
Definition 6.1. Let y ∈ R∞ and let (ϑbn )n∈N be a sequence of estimators. Denote by ϑbn (zn ) the value
of ϑbn when it is calculated using the deterministic realization zn ∈ Rn .
(a) (ϑbn )n∈N is called resistant at y if for every ε > 0 there exists a δ > 0 such that
o
n
sup kϑbn (zn ) − ϑbn (wn )k : zn , wn ∈ Bδ (yn ) ≤ ε ∀n ∈ N,

(6.1)

where Bδ (x) denotes an open ball with center x and radius δ with respect to the metric


#{i ∈ {1, . . . , n} : |zni − wni | ≥ ε }
n
n
dn (z , w ) = inf ε :
≤ε .
n

(b) (ϑbn )n∈N is called asymptotically resistant at y if for any ε > 0 there exists a δ > 0 and N0 (ε , y) ∈
N such that (6.1) holds for n ≥ N0 (ε , y).
(c) For Q ∈ P(R∞ ) we say that (ϑbn )n∈N is strongly resistant at Q if
o
n
Q y ∈ R∞ : (ϑbn )n∈N is resistant at y = 1.
24

(d) (ϑbn )n∈N is called asymptotically strongly resistant at Q if
o
n
Q y ∈ R∞ : (ϑbn )n∈N is asymptotically resistant at y = 1.

(e) (ϑbn )n∈N is called weakly resistant at Q if for any ε > 0 there exists a δ > 0 such that
n
n
o
o
Q y ∈ R∞ : sup kϑbn (zn ) − ϑbn (wn )k : zn , wn ∈ Bδ (yn ) ≤ ε
≥ 1 − ε ∀ n ∈ N.

(f) (ϑbn )n∈N is called asymptotically weakly resistant at Q if for any ε > 0 there exist a δ > 0 and
N(ε ) ∈ N such that
n
n
o
o
∞
n
n
n
n
n
b
b
Q y ∈ R : sup kϑn (z ) − ϑn (w )k : z , w ∈ Bδ (y ) ≤ ε
≥ 1 − ε ∀ n ≥ N(ε ).

With this definition at hand, we want to study the question whether our indirect estimator for the parameters of a CARMA processes is resistant. We will make use of the fact that the indirect estimator
consists of two independent parts, the GM-estimator for the parameters of the auxiliary AR representation, which deals with possible outliers in the observations, and the outlier–free estimator of the AR
parameters based on simulated data.
Theorem 6.2. The GM-estimator (πbnGM (ϑ0 ))n∈N is strongly resistant at PY (h) .

Proof. First of all, (πbnGM (ϑ0 ))n∈N is asymptotically strongly resistant at PY (h) . This follows from
(Boente et al., 1987, Theorem 5.1). The theorem requires that φ and χ fulfill Assumption E, (E2’)
and that the limiting equation has a unique solution, which we assumed. Moreover, (Ymh )m∈Z is
stationary and ergodic due to (Marquardt and Stelzer, 2007, Proposition 3.34) and fulfills (D.3).
Next, by (Cox, 1981, Lemma 5) πbnGM (ϑ0 ) is a continuous function of Y n for every n ∈ N. This in
combination with the asymptotically strongly resistance of (πbnGM (ϑ0 ))n∈N at PY (h) implies the strong
resistance due to (Boente et al., 1987, Proposition 4.2).
Theorem 6.3. The indirect estimator (ϑbnInd )n∈N is weakly resistant and asymptotically weakly resistant at PY (h) .
Proof. Let ε > 0. Since QInd has a unique minimum in ϑ0

η := 3−1

inf

kϑ −ϑ0 k>ε /4

QInd (ϑ ) > 0.

Moreover, the map x 7→ xT Ω x is continuous and hence, uniformly continuous on the compact set
Π ′ = π (Θ ). Thus, there exists an ε ′ > 0 such that
T

sup |π T Ω π − π ′ Ω π ′ | ≤

π ,π ′ ∈Π ′
kπ −π ′ k≤ε ′

η
.
8

Define

Ω0 (ε , δ ) :=

\

n∈N

{y ∈ R∞ : sup {kπbn (zn ) − πbn (wn )k : zn , wn ∈ Bδ (yn )} ≤ ε /2} .

25

Due to the strong resistance of the GM-estimator πbn = πbnGM (ϑ0 ) at PY (h) given in Theorem 6.2 and
(Boente et al., 1987,
Proposition 4.1) there exists an δ > 0 such that for
e
ε := min{ε , ε ′ , η /(8kΩ k sup π ∈Π ′ kπ k)}
P(Ω0 (e
ε , δ )) ≥ 1 −

e
ε
ε
≥ 1− .
2
2

Let y ∈ Ω0 (e
ε , δ ) and zn ∈ Bδ (yn ). Then for any n ∈ N,
sup |LInd (ϑ , yn ) − LInd (ϑ , zn )|

ϑ ∈Θ

=

sup | − 2[πbn (yn ) − πbn (zn )]T Ω πbnS (ϑ ) + πbn (yn )T Ω πbn (yn ) − πbn (zn )T Ω πbn (zn )|

ϑ ∈Θ

≤

sup 2kπ kkΩ kkπbn (yn ) − πbn (zn )k + |πbn (yn )T Ω πbn (yn ) − πbn (zn )T Ω πbn (zn )| ≤

π ∈Π ′

η
.
4

(6.2)

P

On the other hand, from Theorem 3.2(a) we know that supϑ ∈Θ |LInd (ϑ , Y n ) − QInd (ϑ )| → 0. Hence,
there exists an N(ε ) ∈ N such that


ε
n
P sup kLInd (ϑ , Y ) − QInd (ϑ )| ≤ η > 1 −
∀ n ≥ N(ε ).
2
ϑ ∈Θ

ε , δ ) and n ≥ N(ε ).
Define Ωn := {y ∈ R∞ : supϑ ∈Θ |LInd (ϑ , yn )− QInd (ϑ )| ≤ η } and let y ∈ Ωn ∩ Ω0 (e
Then
|LInd (ϑ0 , yn )| = |LInd (ϑ0 , yn ) − QInd (ϑ0 )| ≤ η ,

(6.3)

and with the definition of η we receive
inf

|ϑ −ϑ0 |≥ ε4

|LInd (ϑ , yn )| ≥

inf

|ϑ −ϑ0 |≥ ε4

|QInd (ϑ )| − sup |LInd (ϑ , yn ) − QInd (ϑ )| ≥ 3η − η = 2η . (6.4)
ϑ ∈Θ

Since ϑbnInd (yn ) minimizes LInd (ϑ , yn ) we can deduce from (6.3) and (6.4) that

ε
kϑbnInd (yn ) − ϑ0 k < .
4

(6.5)

Let zn ∈ Bδ (yn ). Due to (6.2) and (6.3) we obtain

|LInd (ϑ0 , zn )| ≤ sup |LInd (ϑ , zn ) − LInd (ϑ , yn )| + |LInd (ϑ0 , yn )| ≤
ϑ ∈Θ

η
5η
+η =
.
4
4

(6.6)

Likewise, (6.2) and (6.4) give us that
inf

|ϑ −ϑ0 |≥ ε4

|LInd (ϑ , zn )| ≥

inf

|ϑ −ϑ0 |≥ ε4

|LInd (ϑ , yn )| − sup |LInd (ϑ , zn ) − LInd (ϑ , yn )| ≥
ϑ ∈Θ

Since ϑbnInd (zn ) minimizes LInd (ϑ , zn ) we can conclude from (6.6) and (6.7) that

ε
kϑbnInd (zn ) − ϑ0 k < .
4

26

7η
.
4

(6.7)

(6.8)

Finally, (6.5) and (6.8) result in

ε
kϑbnInd (zn ) − ϑbnInd (yn )k ≤ kϑbnInd (zn ) − ϑ0 k + kϑbnInd (yn ) − ϑ0 k < .
2

ε , δ ) ∩ Ωn ) ≥ 1 − ε , and for y ∈ Ω0 (e
ε , δ ) ∩ Ωn and
To summarize, for n ≥ N(ε ) we have PY (h) (Ω0 (e
n
n
n
z , w ∈ Bδ (y ) we have
kϑbnInd (zn ) − ϑbnInd (wn )k ≤ kϑbnInd (zn ) − ϑbnInd (yn )k + kϑbnInd (yn ) − ϑbnInd (wn )k < ε .

This gives the asymptotically weakly resistance of (ϑbnInd )n∈N at PY (h) .
By definition, ϑbnInd depends on Y n through a continuous function applied to πbnGM (ϑ0 ) and therefore,
Ind
b
ϑn is a continuous function in Y n . This and the asymptotically weakly resistance of (ϑbnInd )n∈N at
PY (h) imply the weakly resistance at PY (h) by (Boente et al., 1987, Proposition 4.2).
Remark 6.4. If the stronger version supϑ ∈Θ kπbnS (ϑ ) − π (ϑ )k → 0 P-a.s. holds then it is possible to
show on a similar way that the indirect estimator (ϑbnInd )n∈N is even strongly resistant.

As already mentioned, one could also define qualitative robustness of a sequence of estimators by
demanding that the distribution of the estimator does not change too much when the data is changed
slightly. To make this notion explicit, we first define a pseudometric for measures on metric spaces.

Definition 6.5. For a metric space (M, d) with Borel sets B(M), the Prokhorov distance πd between
two measures µ , ν on B(M) with respect to d is defined as

πd (µ , ν ) := inf{ε > 0 : µ (A) ≤ ν ({x ∈ M : d(x, A) < ε }) + ε ∀ A ∈ B(M)}.
This pseudometric is a key component of the definition of qualitative robustness.
Definition 6.6. Let dΘ be a metric on Θ and let ρn be a pseudometric on P(Rn ), n ∈ N. For P ∈
P(R∞ ) denote by Pn the n–th order marginal of P. Finally, Pϑbn ∈ P(Θ ) is the distribution of the
estimator ϑbn under Pn . Then the sequence of estimators (ϑbn )n∈N is called ρn –robust at P if for every
ε > 0 there exists a δ > 0 such that for every Qn ∈ P(Rn ) with ρn (Pn , Qn ) < δ :

πdΘ (Pϑbn , Qϑbn ) ≤ ε .

As shown in (Boente et al., 1987, Theorem 3.1), this is a direct generalization of the definition of
π –robustness given by Hampel (1971) for i.i.d. processes.
Theorem 6.7. On Rn we define the metric
dn (xn , yn ) = inf {ε : #{i : |xi − yi | ≥ ε }/n ≤ ε }
and use the Prokhorov distance with respect to dn ,

πdn (Pn , Qn ) = inf{ε > 0 : Pn (A) ≤ Qn ({xn ∈ Rn : dn (xn , A) < ε }) + ε ∀ A ∈ B(Rn )}
on B(Rn ). Then the sequence of estimators (ϑbnInd )n∈N is πdn –robust at PY (h) .

Regarding the metric dn two points are close if all coordinates except a small fraction are close.

27

Proof. Due to Theorem 6.3 the estimator (ϑbnInd )n∈N is weakly resistant at PY (h) . A conclusion of
(Boente et al., 1987, Theorem 4.2(i)) is then the πdn –robustness of (ϑbnInd )n∈N at PY (h) .

In summary, we can say that our indirect estimator is weakly resistant at PY (h) as well as πdn –robust.
This is in contrast to, e.g., M–estimators, which are not qualitatively robust even in the case of linear
regression (cf. (Maronna and Yohai, 1981, p.8)).

6.2 The influence functional
We continue our investigation of robustness of ϑbnInd with the study of the influence functional of the indirect estimator. Intuitively speaking, the influence functional measures the change in the asymptotic
bias of an estimator caused by an infinitesimal amount of contamination in the data. This measure
of robustness was originally introduced as the influence curve by Hampel (1974) for i.i.d. processes.
It was later generalized to the time series context by Künsch (1984) who explicitly studies the estimation of autoregressive processes. However, in the paper of Künsch only estimators which depend on a finite–dimensional marginal distribution of the data–generating process and a very specific
form of contamination are considered. To remedy this, a further generalization was then made by
Martin and Yohai (1986) who consider the influence functional and explicitly allow for the estimators to depend on the measure of the process which makes more sense in the time series setup (cf.
(Martin and Yohai, 1986, Section 4)). We work with their definition in the following.
In addition to the assumptions given at the beginning of Section 6 we assume throughout this section
γ
γ
γ
that there exists a unique solution π GM (ϑ0 ) of (4.2) for (Ymh )m∈Z = (Ymh (ϑ0 ))m∈Z for any 0 ≤ γ ≤ 1,
and JGM (ϑ0 ) is non-singular. We denote the probability measure associated to the distribution of
γ
γ
(Ymh )m∈Z by PY (h) for 0 ≤ γ ≤ 1. Note that γ = 0 corresponds to the case where there are no outliers,
i.e., we can observe the nominal process without error and then write PY0 (h) = PY (h) . Similarly PZ is the
γ
γ
distribution of (Zm )m∈Z and PW is the distribution of (Wm )m∈Z . We write {PY (h) } := {PY (h) , 0 ≤ γ ≤
1} ⊆ P(R∞ ) and introduce the statistical functional
γ

TGM : {PY (h) } → Π

as

γ

γ

PY (h) 7→ π GM (ϑ0 ).

Then, the definition of the influence functional for the GM-estimator is
γ

γ

IFGM (PZ , {PY (h) }) := lim

γ →0

TGM (PY (h) ) − TGM (PY (h) )

γ

π GM (ϑ0γ ) − π (ϑ0 )
γ →0
γ

= lim

(6.9)

whenever this limit is well–defined. Note that the influence functional depends on the whole “arc”
γ
of contaminated measures {PY (h) , 0 ≤ γ ≤ 1}. This is the most important difference to the definition
γ
used by Künsch (1984), because in that paper the approximation PY (h) = (1 − γ )PY (h) + γν for some
∞
fixed ν ∈ P(R ) is used (Künsch, 1984, Eq. (1.11)). The influence functional measures the effect
of an infinitesimal contamination of the true process by the process (Zm ) on the asymptotic estimate
defined via the functional TGM .
In a similar vein, we can define the influence functional for the estimation of the parameter ϑ0 of
our CARMA process. Analogous to TGM , we first define a suitable statistical functional
γ

TInd : {PY (h) } → Θ

as

γ

γ

γ

PY (h) 7→ ϑ0Ind (γ ) := arg min[π (ϑ ) − π GM (ϑ0 )]T Ω [π (ϑ ) − π GM (ϑ0 )].
ϑ ∈Θ

This is the analog of ϑ0 = arg minϑ ∈Θ [π (ϑ ) − π (ϑ0 )]T Ω [π (ϑ ) − π (ϑ0 )] in the uncontaminated case
(cf. (3.1)). With this and ϑ0Ind (0) = ϑ0 due to π GM (ϑ0 ) = π (ϑ0 ) the definition of the influence

28

functional of the indirect estimator is
γ

γ

IFInd (PZ , {PY (h) }) = lim

γ →0

TInd (PY (h) ) − TInd (PY (h) )

γ

ϑ0Ind (γ ) − ϑ0
.
γ →0
γ

= lim

We are interested in properties of this functional, in particular, if it is bounded. Boundedness of the
influence functional implies that the estimate arising from the contaminated process cannot move too
far away from the one in the uncontaminated case if the rate of contamination is infinitesimal. This
property is well–known for the influence functional for GM-estimators of AR processes. Since these
estimators are an integral building block of the indirect estimator, one can hope that it carries over to
our scenario and indeed it does, since the two functionals are proportional.
Proposition 6.8. Assume that ∇ϑ π (ϑ0 ) has full column rank N(Θ ). Then if IFGM exists then IFInd
exists and
γ
γ
IFInd (PZ , {PY (h) }) = H (TInd (PY (h) ))IFGM (PZ , {PY (h) }),
where
H (TInd (PY (h) )) = H (ϑ0 ) = [∇ϑ π (ϑ0 )T Ω ∇ϑ π (ϑ0 )]−1 ∇ϑ π (ϑ0 )T Ω .
Proof. This follows from (de Luna and Genton, 2000, Theorem 1) as special case.
From this theorem we see that the question of the boundedness of the influence functional for the
indirect estimator of a discretely sampled CARMA process reduces to the question of the boundedness
of the influence functional for the GM-estimatior of the auxiliary AR representation of the sampled
CARMA process.
Theorem 6.9. Let the additive outlier model be given and let Assumption D hold. Then there exists a
constant K > 0 such that
γ

kIFGM (PW , {PY (h) })k ≤ 2(r + 1)K JGM (ϑ0 )−1 .
Proof. The plan is to apply (Martin and Yohai, 1986, Theorem 4.3). Therefore, we have to check
that (Martin and Yohai, 1986, Eq. (4.6)) hold. Sufficient conditions for this equation are given
in (Martin and Yohai, 1986, Theorem 4.2) which are obviously satisfied in our case due to (E.2),
γ
(E.6) and due to the fact that π GM (ϑ0 ) depends only on the distribution of the finite random vector
γ
γ
(Yh , . . . ,Y(r+1)h ). For the same reasons and with our assumption that J (ϑ0 ) is non-singular the other
conditions in (Martin and Yohai, 1986, Theorem 4.3) are satisfied as well.
As well due to (Martin and Yohai, 1986, Theorem 4.3) it is possible to give an explicit (but not a very
γ
γ
handy) representation of IFGM (PW , {PY (h) }) and hence, though Proposition 6.8 for IFInd (PW , {PY (h) }).

6.3 The breakdown point
The breakdown point is (for a sample of data with fixed length n) the maximum percentage of outliers which can be contained in the data without ”ruining” the estimator. In this sense, it measures
how much the observed data can deviate from the nominal model before catastrophic effects in the
estimation procedure happen. However, the formal definition depends on the model and the estimator. Maronna and Yohai (1991) and Maronna et al. (1979) deal explicitly with the breakdown point of
GM-estimators in regression models and Martin and Yohai (1985) and Martin (1980) study it in the
time series context. A very general definition of the breakdown point is given in (Genton and Lucas,

29

2003, Definition 1 and Definition 2). Heuristically speaking, the fundamental idea of that definition
is that the breakdown point is the smallest amount of outlier contamination with the property that
the performance of the estimator does not get worse anymore if the contamination is increased further. As already mentioned in (Martin, 1980, p. 239) (the proof is given in the unpublished paper
of Martin and Jong (1977)), and later in (de Luna and Genton, 2001, p. 377) and (Genton and Lucas,
2003, p. 89), the breakdown point of the GM-estimator applied to estimate the parameters of an
AR(r) process is 1/(r + 1). Hence, the breakdown point of our indirect estimator is as well 1/(r + 1)
since the other building block of the indirect estimator, the estimator πbnS (ϑ ) is applied to a simulated
outlier–free sample.

7 Simulation study
We simulate a CARMA process on the interval [0, 1000] and choose a sampling distance of h = 1,
resulting in n = 1000 observations of the discrete–time process. The simulated processes are driven
either by a standard Brownian motion or by a univariate NIG Lévy process. For the NIG Lévy process
we use the parameters α = 3, β = 1, ∆ = 1, δ = 2.5145 and µ = −0.8890. These parameters result in
a zero–mean Lévy process with variance approximately 1 which allows for comparison of the results
to the Brownian motion case.
The indirect estimator is defined as in Section 5. We take πbn as GM-estimator πbnGM (ϑ0 ) using
the S–Plus software which provides a pre–built function ar.gm for applying GM-estimators to AR
processes. This function uses a Mallows estimator as in Example 4.3(a). The weight function w(y) is
the Tukey bisquare function from Example 4.3(b) applied to kyk, for the function ψ (u) the user can
choose between the Huber ψk –function and the bisquare function. The function is implemented as an
iterative least squares procedure as described by (Martin, 1980, p. 231ff.) and therefore also allows
to use first some iteration steps with the Huber ψk –function and then some steps with the bisquare
function. As advocated by Martin (1980), our experience is that doing 6 iterations using the Huber
function and then 2 steps with the bisquare function works well. In our experiments we use k = 4
for the tuning constant of the ψk –function. In general, we set s = 75 to obtain the simulation–based
S (ϑ )) in the simulation part of the indirect procedure. The
observations YSsn (ϑ ) = (YhS (ϑ ), . . . ,Ysnh
type of Lévy process used for the simulation is of the same type as the Lévy process driving the
CARMA process. For the estimator πbnS (ϑ ) we apply the least squares estimator and as weighting
matrix Ω we take the identity matrix for convenience reasons. Some experiments in which we first
estimated the asymptotic covariance matrix of the GM-estimator by the empirical covariance matrix
of a suitable number of independent realizations of πbn and set Ω to be the inverse of that estimate did
not significantly affect the procedure positively or negatively so that the use of the convenient identity
matrix seems justified. For the outlier model we choose the process (Vm )m∈Z as i.i.d. Bernoulli
random variables where P(V1 = 1) = γ varies. In all but one of the studies, the process (Zm )m∈Z
is chosen to be constant Zm = ξ for m ∈ Z. We use varying values of ξ . In the experiment where
(Zm )m∈Z has a different structure, this will be mentioned explicitly. In each experiment, we calculate
the indirect estimator and, for comparison purposes, the QMLE as defined in Schlemm and Stelzer
(2012) in 50 independent replications and report on the average estimated value, the bias relative to
the true parameters and the empirical variance of the parameter estimates.

30

7.1 CARMA processes driven by a Brownian motion
(1)

In a first experiment, we use as true process a CARMA(1,0) process with parameter ϑ0 = −2. This
process is of particular interest because its discretely sampled version admits an AR(1) representation.
For this reason, one would expect the procedure to work very well here as the auxiliary representation
is actually exact. Naturally, we use r = 1 in this case. We consider three different scenarios of outlier
contamination. In the first case, we set ξ = 5 and γ = 0.1. In the second, we set ξ = 10 and γ = 0.1,
while for the last one we choose ξ = 5 and γ = 0.15. Note that already ξ = 5 represents quite large
outliers, since for a sample path in this situation we typically observe that the values of the discretely
sampled process lie between −3.5 and 3.5. The results of the simulation studies are given in Table 1.
ξ = 5, γ = 0.1

ϑ1

Mean
−2.4497

MLE
Bias
−0.4497

ϑ1

Mean
−5.9097

MLE
Bias
−3.9097

ϑ1

Mean
−2.7873

MLE
Bias
−0.7873

Variance
Mean
0.0559
−2.0768
ξ = 10, γ = 0.1
Variance
Mean
0.2360
−2.0245
ξ = 5, γ = 0.15
Variance
0.3860

Mean
−2.1703

Indirect
Bias
−0.0768

Variance
0.0513

Indirect
Bias
−0.0245

Variance
0.0663

Indirect
Bias
−0.1703

Variance
0.0973

(1)

Table 1: Results for ϑ0 .
As we can see, already in the case ξ = 5 and γ = 0.1, the indirect estimator performs vastly better
than the QMLE giving a much less biased estimate at a similar variance. In the situation where ξ is
increased to 10 the QMLE has lost all its information about the true parameter and provides no useful
estimate anymore. On the other hand, the indirect estimator stays close to the true value (we explain
the even smaller bias in comparison to the situation with ξ = 5 as caused by the relatively small
number of 50 iterations. Of course one should not expect the bias to systematically decrease when
ξ increases), while the variance has increased only slightly. Increasing γ to 0.15 but keeping ξ = 5
shows that both estimators perform worse than in the situation with γ = 0.1, which is to be expected.
But once again, the indirect estimator deals much better with the higher outlier percentage than the
QMLE. Comparing these studies, we see that for the indirect estimator the percentage of outliers has
a bigger effect on the estimates than the actual size of the outliers.
In our next study we investigate a CARMA(3,1) process. This especially means that the sampled
process is not a weak AR process anymore. The true parameter is

(2)
ϑ0 = −1 −2 −2 0 1 .

For this process, we choose r = 5, which is also the minimum order of the auxiliary AR representation
(2)
to satisfy Assumption B. We do four experiments in this setup. We estimate ϑ0 for each of the
following contamination configurations: ξ = 5 and γ = 0.1, ξ = 10 and γ = 0.1, ξ = 5 and γ = 1/6,
and ξ = 5 and γ = 0.25. Remember that in this situation, the breakdown point has an upper bound
of 1/6 since we have r = 5. Hence, γ = 0.25 lies above the breakdown point and we would expect to
encounter problems in the estimation procedure, while for γ ≤ 1/6 these problems should not occur.

31

We will see that this is indeed the case. The results of the four experiments are given in Table 2.
ξ = 5, γ = 0.1
MLE
Bias
0.3124
−0.6307
−1.3831
2.5467
−0.4379

ϑ1
ϑ2
ϑ3
ϑ4
ϑ5

Mean
−0.6876
−2.6307
−3.3831
2.5467
0.5621

ϑ1
ϑ2
ϑ3
ϑ4
ϑ5

Mean
−0.4614
−1.4955
−1.8424
3.3178
2.6317

ϑ1
ϑ2
ϑ3
ϑ4
ϑ5

Mean
−0.4445
−2.3450
−3.5119
3.1446
0.7903

MLE
Bias
0.5555
−0.3450
−1.5119
3.1446
−0.2097

ϑ1
ϑ2
ϑ3
ϑ4
ϑ5

Mean
−0.1022
−1.5663
−4.2751
4.2080
0.8238

MLE
Bias
0.8978
0.4337
−2.2751
4.2080
−0.1762

MLE
Bias
0.5386
0.5045
0.1576
3.3178
1.6317

Variance
Mean
0.0294
−1.0174
0.2550
−1.9930
0.0573
−1.9954
0.0190
0.0048
0.0498
1.0007
ξ = 10, γ = 0.1
Variance
Mean
0.0478
−1.0202
0.1616
−2.0063
0.0953
−1.9802
0.1280
0.0047
0.0477
0.9987
ξ = 5, γ = 61
Variance
Mean
0.0029
−1.0052
0.1135
−1.9815
0.1823
−2.0210
0.0376
0.0106
0.1317
1.0043
ξ = 5, γ = 0.25
Variance
0.0029
0.1135
0.1823
0.0376
0.1317

Mean
−0.0638
−2.9026
−2.1834
1.1326
1.8766

Indirect
Bias
−0.0174
0.0070
0.0046
0.0048
0.0007
Indirect
Bias
−0.0202
−0.0063
0.0198
0.0047
−0.0013

Variance
0.0114
0.0068
0.0194
0.0040
0.0064

Variance
0.0110
0.0107
0.0259
0.0066
0.0074

Indirect
Bias
−0.0052
0.0185
−0.0210
0.0106
0.0043

Variance
0.0071
0.0136
0.0276
0.0057
0.0054

Indirect
Bias
0.9362
−0.9026
−0.1834
1.1326
0.8766

Variance
0.0024
0.0008
0.0094
0.0018
0.0021

(2)

Table 2: Results for ϑ0 .
For the first two experiments where γ = 0.1, we immediately recognize the maximum likelihood
estimate as basically worthless being severely biased and very far from the true parameter value.
Especially the inclusion of a zero component in the true parameter seems to pose a major problem
since this component is affected by the most bias. On the other hand, the indirect estimator is still
very close to the true parameter value in all components including the zero component. Increasing ξ
to 10 while keeping γ = 0.1 results in a slightly worse performance of the indirect estimator. However,
the increase in the bias is not very substantial and the estimates are still reasonably close to the true
values. This reaffirms that the indirect estimation procedure works in practical scenarios and the
results in the former experiments were not due to the use of the CARMA(1,0) process.
The increase of γ from 0.1 to 1/6 also affects the performance of the indirect estimator. For all
(2)
components of ϑ0 , except the first one, the bias and the variance of the indirect estimator increase.
However, the loss in quality of the indirect estimator is manageable and the calculated estimates still

32

resemble the true parameter. This means that even at the breakdown point of 1/6, the performance
of the indirect estimator is satisfying, although of course not as good as for lower contamination
probabilities.
The situation is vastly different in the experiment with γ = 0.25 where γ is above the breakdown
point. Here, we see not surprisingly that the indirect estimator, too, gives estimates which are severely
biased and are quite far away from the true parameters. We also observe that the numerical procedure
used to obtain the parameter estimates quite often fails to deliver a result at all because the algorithm
terminates with an error. The error occurs when the estimated value of ϑ0 is not an element of Θ
anymore. The results in the table are averaged over experiments in which the algorithm did deliver
a result, the failed attempts were discarded. The ratio of successful to unsuccessful experiments was
roughly equal to 1:2, i.e., the algorithm failed about twice as often as it succeeded. In this sense, we
can say that the estimator has broken down: for a given outlier–contaminated sample, it either does
not return an admissible estimate at all, or, if it does, the estimate is far away from the true parameter.
The latter statement is also evident from the fact that the variances of the indirect estimates are far
smaller in this case than in the other experiments which intuitively means that the algorithm typically
returns very similar bad estimates if it returns a result at all.
Lastly, we use a CARMA(2,1) process with true parameter

(3)
ϑ0 = −0.5 −1 −2 .

We report on two simulation studies, which are a bit different than the ones conducted before. In
the first study, we take γ = 0 and compare the performance of the indirect estimator to that of the
QMLE in the uncontaminated case. In the next study, the process (Zm )m∈Z is a sequence of i.i.d.
random variables with P(Z1 = ξ ) = P(Z1 = −ξ ) = 21 , i.e., every time an outlier appears the sign is
chosen randomly with equal probability. This change in model serves to study the performance of
the procedure in circumstances that are a bit more complicated than the simple additive outlier model
with fixed size ξ . We choose γ = 0.1 and ξ = 10 in this study. In both studies r = 3. The results are
given in Table 3.

ϑ1
ϑ2
ϑ3

Mean
−0.5112
−0.9969
−2.0427

ϑ1
ϑ2
ϑ3

Mean
−0.3796
−2.0699
−4.0399

γ =0
MLE
Indirect
Bias
Variance
Mean
Bias
−0.0112
0.0071
−0.5080 −0.0080
0.0031
0.0056
−1.0019 −0.0019
−0.0427
0.0246
−2.0156 −0.0156
P(Z1 = 10) = P(Z1 = −10) = 0.5, γ = 0.1
MLE
Indirect
Bias
Variance
Mean
Bias
0.1204
0.0585
−0.4908
0.0092
−1.0699
0.2484
−1.0479 −0.0479
−2.0399
0.1408
−2.0015 −0.0015

Variance
0.0086
0.0050
0.0112

Variance
0.0134
0.0132
0.0091

(3)

Table 3: Results for ϑ0 .
In the situation without outliers, both estimators are very close to the true parameter values. The
differences in the bias and the variance are not systematic and probably due to the approximations in
the numerical procedure and only 50 iterations. We see that both estimators can be used to achieve sat(1)
(2)
isfying results. For the parameters ϑ0 and ϑ0 we conducted the same study (i.e. in the outlier–free

33

situation) and observed basically the same results. Using the indirect estimator also yields satisfying results under the more complicated outlier model. This is not surprising since the GM-estimator,
which controls the robustness properties of ϑbnInd , is per construction not sensitive to the sign of the,
outlier, only to the absolute size.

7.2 CARMA processes driven by a NIG Lévy process
In the following experiments, the driving Brownian motion is replaced by the NIG Lévy process.
We repeat some of the experiments of Table 1 and Table 2 using the same outlier configurations and
auxiliary AR orders as in those experiments. The results are given in Table 4 and Table 5.
Mean
−2.3851

ϑ1

MLE
Bias
−0.3851

Variance
0.0650

Mean
−2.0536

Indirect
Bias
−0.0536

Variance
0.0860

(1)

Table 4: Results for ϑ0 , ξ = 5, γ = 0.1, driving NIG Lévy process
ξ = 5, γ = 0.1

ϑ1
ϑ2
ϑ3
ϑ4
ϑ5

Mean
−0.6669
−2.5668
−3.3917
0.5077
2.5503

MLE
Bias
0.3331
−0.5668
−1.3917
−0.4923
2.5503

ϑ1
ϑ2
ϑ3
ϑ4
ϑ5

Mean
−0.4691
−1.5286
−1.8365
2.6651
3.3668

MLE
Bias
0.5309
0.4714
0.1635
1.6651
3.3668

Variance
Mean
0.0311
−1.0113
0.2676
−1.9936
0.0626
−1.9967
0.0535
1.0027
0.0236
−0.0070
ξ = 10, γ = 0.1
Variance
0.0689
0.1200
0.0668
0.0316
0.0981

Mean
−1.0030
−1.9712
−2.0321
1.0056
0.0127

Indirect
Bias
−0.0113
0.0064
0.0033
0.0027
−0.0070
Indirect
Bias
−0.0030
0.0288
−0.0321
0.0056
0.0127

Variance
0.0056
0.0042
0.0096
0.0043
0.0035

Variance
0.0045
0.0112
0.0178
0.0038
0.0056

(2)

Table 5: Results for ϑ0 , γ = 0.1, driving NIG Lévy process
We see that the results do not substantially differ from those in the Brownian motion case. In all
three experiments with the NIG Lévy process, the QMLE, just like in the Brownian motion case,
ceases to be a meaningful estimator in the presence of outliers in the data. The indirect estimator, on
the other hand, continues to provide good estimates.

8 Conclusion
In this paper we presented an indirect estimation procedure for the parameters of a discretely observed
CARMA process by estimating the parameters of its auxiliary AR(r) representation using a GMestimator. Since there does not exist an explicit form of the map between the CARMA parameters
and the AR parameters, an additional simulation step to get back from the AR parameters to the
CARMA parameters was necessary. Sufficient conditions were given such that the indirect estimator

34

is consistent and asymptotically normally distributed, on the one hand, in a general context, but on
the other hand as well for the special case where πbn = πbnGM (ϑ0 ) and πbnS (ϑ ) = πbnLS (ϑ ). Moreover, we
investigated different robustness properties and showed that the indirect estimator is weakly resistant,
πdn -robust and has a bounded influence functional.
Summarizing the simulation studies, the indirect estimator performs convincingly for various orders
p and q of the CARMA process, for different driving Lévy processes and for a variety of outlier
configurations. As soon as outliers are present the use of the indirect estimator instead of the QMLE
seems advisable, since there was no situation in which the performance of the QMLE came close to
being satisfying in contrast to the indirect estimator. Of course, it is clear that the indirect estimator
has its bounds as well. We especially saw that both an increase of γ , the proportion of outliers, and
the size of the outliers affect the performance. Increasing γ too far eventually causes the estimator to
break down because we get above the breakdown point. Unfortunately, the breakdown point 1/(r + 1)
is very low if r is large, which is, however, necessary if the order p of the CARMA(p, q) process is
large.

A Appendix: Proof of Lemma 5.4
Before we state the proof we require some auxiliary results.
Lemma A.1. For any i, j, l ∈ {1, . . . , N(Θ )} the following conditions hold:
(a) supϑ ∈Θ | fϑ (u)| ≤ C e−ρ u and

(b) supϑ ∈Θ
(c) supϑ ∈Θ
(d) supϑ ∈Θ

∂ fϑ (u)
∂ϑj
∂ 2 fϑ (u)
∂ ϑ j ∂ ϑi

R∞

supϑ ∈Θ fϑ (u)2 du < ∞.

2
R
du < ∞.
≤ Cu e−ρ u and 0∞ supϑ ∈Θ ∂ ∂fϑϑ(u)
j
0

≤ Cu2 e−ρ u and

∂ 3 fϑ (u)
∂ ϑ j ∂ ϑi ∂ ϑl

R∞
0

≤ Cu3 e−ρ u and

supϑ ∈Θ

R∞



∂ 2 fϑ (u)
∂ ϑ j ∂ ϑi

supϑ ∈Θ

0



2

du < ∞.

∂ 3 fϑ (u)
∂ ϑ j ∂ ϑi ∂ ϑl

2

du < ∞.

Proof. (a) Due to Remark 2.1(iii) we have that supϑ ∈Θ k eAϑ u k ≤ C e−ρ u and hence, supϑ ∈Θ | fϑ (u)| ≤
kc k supϑ ∈Θ k eAϑ u k ≤ C e−ρ u using the continuity of cϑ on the compact set Θ . Finally,
sup
R ∞ ϑ ∈Θ ϑ
2
0 supϑ ∈Θ fϑ (u) du < ∞.
(b) A consequence of Wilcox (1967) is that


Z u
∂ e Aϑ u
∂
=
eAϑ (u−s)
Aϑ eAϑ s ds,
∂ϑj
∂ϑj
0
A u

and hence, supϑ ∈Θ k ∂∂eϑϑj k ≤ Cu e−ρ u . From this and

∂ fϑ (u)
=
∂ϑj



∂ cϑ
∂ϑj



Aϑ u

e

e p + cϑ




∂ eAϑ u
ep
∂ϑj

we receive
Z ∞

sup

0 ϑ ∈Θ



∂ fϑ (u)
∂ϑj

2

du ≤ C

(c,d) can be proven similarly to (a) and (b).

35

Z ∞
0

u2 e−2ρ u du < ∞.

∗

Lemma A.2. Let (Lt )t∈R be a Lévy process with E|L1 |2N < ∞ for some N ∗ R∈ N and E(L1 ) = 0.
T ∗
∗
∞
φ (u) dLu |2N < ∞
Furthermore, let φ : R → R be a measurable map with φ ∈ Nj=1 L2 j (R). Then E| −∞
and there exist finite constants c j1 ,..., jk such that
E

Z

∞
−∞

φ (u) dLu

2N ∗

N∗

=

∑

∑

k=1

c j1 ,..., jk

j1 +...+ jk =N ∗
j1 ,..., jk ∈N0

Z

∞

φ

2 j1

−∞

 Z
(u) du · · ·

∞

φ

2 jk

−∞


(u) du .

Proof. For N ∗ = 2 the result was already derived in (Cohen and Lindner, 2013, Lemma 3.2), the
proof for general N ∗ uses the same ideas. Let ν be the Lévy measure of (Lt )t∈R and V its Gaussian
parameter. Define
1
ψ (s) = − V 2 s2
2

Z ∞

−∞

φ 2 (u) du +

Z ∞Z ∞

−∞ −∞

[eisφ (u)x −1 − isφ (u)x]ν (dx) ds,

and write ψ (i) (s) for the i-th derivative of ψ (s). Due to the Lévy-Khintchine formula we get


 Z ∞
= exp (ψ (s)) .
ξ (s) := E exp is
φ (u) dLu
−∞

2N ∗
R∞
∗
φ (u) dLu
is obtained as (2N ∗ )-th derivative of ξ (s) at s = 0 times (−1)N . StraightThen E −∞
∗
forward calculations yield that the (2N ∗ )-th derivative ξ (2N ) (s) of ξ (s) has the form



ξ

(2N ∗ )

2N ∗


(s) =  ∑

k=1

∑

i1 +...+ik =2N ∗
i1 ,...,ik ∈N0


cei1 ,...,ik ψ (i1 ) (s) · · · ψ (ik ) (s) exp(ψ (s)).

Plugging in for s = 0 and taking into account that ψ (0) = 1 and ψ (1) (0) = 0 gives
E

Z

∞

−∞

φ (u) dLu

2N ∗

N ∗ (2N ∗ )

= (−1) ξ

N∗

(0) =

∑
k=1

∑

j1 +...+ jk =N ∗
j1 ,..., jk ∈N0

c j1 ,..., jk ψ (2 j1 ) (0) · · · ψ (2 jk ) (0).

Finally, with

 Z ∞
2
x ν (dx)
ψ (0) = −V −
φ (u) du ,
−∞
−∞

Z ∞
 Z ∞
(2i)
i1
2i
2i
x ν (dx)
ψ (0) = (−1)
φ (u) du ,
(2)



2

Z ∞

2

−∞

−∞

i ≥ 2,

we obtain the result.
∗

Lemma A.3. Let N ∗ ∈ N be such that 2N ∗ > N(Θ ) and E|L1 |2N < ∞. For any i, j, l ∈ {1, . . . , N(Θ )}
the following maps are P-a.s. Hölder-continuous of order γ ∈ [0, 1 − N(Θ )/(2N ∗ )):
(a) ϑ 7→
(b) ϑ 7→
(c) ϑ 7→

R∞
0

fϑ (u) dLu =: Z(ϑ ),

R∞ ∂
( j)
0 ∂ ϑ fϑ (u) dLu =: Z (ϑ ),
j

R∞

∂
0 ∂ ϑ j ∂ ϑi fϑ (u) dLu

=: Z ( j,i) (ϑ ).

36

∗

Moreover, E(supϑ ∈Θ |Z(ϑ )|2N ) < ∞ and for U := sup 0<kϑ1 −ϑ2 k<1
ϑ1 ,ϑ2 ∈Θ

The same is true for Z ( j) (ϑ ) and Z ( j,i) (ϑ ).

|Z(ϑ1 )−Z(ϑ2 )|
kϑ1 −ϑ2 kγ

∗

we have EU 2N < ∞.

Proof.
(a) Let ϑ1 , ϑ2 ∈ Θ and define φ (u) := fϑ1 (u) − fϑ2 (u). Due to a Taylor expansion we obtain

φ (u) = fϑ1 (u) − fϑ2 (u) = ∇ϑ fϑe(u) (u)(ϑ1 − ϑ2 )

for some ϑe(u) ∈ Θ with kϑe(u) − ϑ2 k ≤ kϑ1 − ϑ2 k. Hence, we receive by Lemma A.1
|φ (u)| ≤ sup k∇ϑ fϑ (u)kkϑ1 − ϑ2 k ≤ Cu e−ρ u kϑ1 − ϑ2 k.
ϑ ∈Θ

Plugging this into Lemma A.2 gives
E(Z(ϑ1 ) − Z(ϑ2 ))2N
N∗

≤ C∑

k=1

∗

∑

j1 +...+ jk =N ∗
j1 ,..., jk ∈N0

|c j1 ,..., jk |

Z

∞

−∞

2 j1 2 j1 −2 j1 ρ

kϑ1 − ϑ2 k

e

u

 Z
du · · ·

∞
−∞

2 jk 2 jk −2 jk ρ

kϑ1 − ϑ2 k

u

e


du

∗

≤ Ckϑ1 − ϑ2 k2N .
Then, an application of Kolmogorov-Chentsov Theorem (cf. (Schilling and Partzsch, 2014, Theo∗
rem 10.1)) yield the Hölder continuity and EU 2N < ∞. Since Θ is compact, some straightforward
∗
calculations yield E(supϑ ∈Θ |Z(ϑ )|2N ) < ∞.
The proofs of (b)-(c) are similarly to the proof of (a) and thus, skipped.
Lemma A.4. Let [a, b] ⊆ R be a bounded interval, (Lt )t∈R be a Lévy process with finite second moments and let g : [a, b] × R → R be differentiable in the first component with derivative ∂ g(∂ϑϑ,u) . Moreover, assume
(a)

∂ g(ϑ , u)
is bounded B([a, b]) ⊗ B([−u1 , u2 ])-measurable for all u1 , u2 > 0,
∂ϑ

(b) lim sup

Z

∂ g(ϑ , u)
∂ g(ϑ , u)
du = 0 and lim sup
= 0,
∂ϑ
∂ϑ
|u|→∞ ϑ ∈[a,b]

M→∞ ϑ ∈[a,b] |u|>M

(c) ϑ 7→
(d) ϑ 7→

Z ∞

−∞

g(ϑ , u) dLu is P-a.s. continuous,

Z ∞
∂ g(ϑ , u)
−∞

∂ϑ

dLu is P-a.s. continuous.

Then, outside a P-zero set, Z(ϑ , ω ) :=
val (a, b) and

∂
∂ϑ

Z ∞

−∞

R∞

−∞ g(ϑ , u) dLu (ω )

g(ϑ , u) dLu (ω ) =

is continuous differentiable over the inter-

Z ∞
∂ g(ϑ , u)
−∞

37

∂ϑ

dLu (ω ).

Proof. An application of Fubini’s Theorem for Lévy-integrals (see (Brockwell and Schlemm, 2013,
Theorem 2.4)) gives that for ϑe ∈ [a, b]
Z ϑe Z ∞
∂ g(ϑ , u)
a

−∞

∂ϑ

dLu d ϑ =

Z ∞ Z ϑe
∂ g(ϑ , u)
−∞ a

∂ϑ

d ϑ dLu

P-a.s.

The remaining of the proof follows the same line as (Hutton and Nelson, 1984, Theorem 2.2).
Proof of Lemma 5.4. A combination of Lemma A.1-Lemma A.4 result in the MA representation. A
conclusion of the MA representation and (Fuchs and Stelzer, 2013, Theorem 3.5) is that the process
is mixing and therefore, in particular, ergodic.

References
Ben, M., Martinez, E. and Yohai, V. (1999). Robust estimation in vector autoregressive movingaverage models. J. Time Ser. Anal., 20(4), 381–399.
Bernstein, D. (2009). Matrix Mathematics: Theory, Facts, and Formulas. Princeton Univ. Press.
Billingsley, P. (1999). Convergence of probability measures. 2nd edn. John Wiley & Sons, Inc., New
York.
Boente, G., Fraiman, R. and Yohai, V. (1987). Qualitative robustness for stochastic processes. Ann.
Statist., 15(3), 1293–1312.
Bradley, R. (2007). Introduction to Strong Mixing Conditions. Volume 1. Kendrick Press, Heber City.
Brockwell, P. (2001). Lévy-driven CARMA processes. Ann. Inst. Statist. Math., 53, 113–124.
Brockwell, P. (2014). Recent results in the theory and applications of CARMA processes. Ann. Inst.
Statist. Math., 66(4), 647–685.
Brockwell, P. and Davis, R. (1991). Time Series: Theory and Methods. 2nd edn. Springer, New York.
Brockwell, P., Davis, R. and Yang, Y. (2011). Estimation for non-negative Lévy-driven CARMA
processes. J. Bus. Econom. Statist., 29(2), 250–259.
Brockwell, P. and Lindner, A. (2009). Existence and uniqueness of stationary Lévy-driven CARMA
processes. Stochastic Process. Appl., 119(8), 2660–2681.
Brockwell, P. and Schlemm, E. (2013). Parametric estimation of the driving Lévy process of multivariate CARMA processes from discrete observations. J. Multivariate Anal., 115, 217–251.
Bustos, O. (1982). General M-estimates for contaminated pth-order autoregressive processes: consistency and asymptotic normality. Z. Wahrscheinlichkeit, 59(4), 491–504.
Bustos, O. and Yohai, V. (1986). Robust estimates for ARMA models. J. Amer. Statist. Assoc.,
81(393), 155–168.
Cohen, S. and Lindner, A. (2013). A central limit theorem for the sample autocorrelations of a Lévy
driven continuous time moving average process. J. Statist. Plann. Inference, 143(8), 1295–1306.

38

Cox, D. (1981). Metrics on stochastic processes and qualitative robustness. University of Washington,
Dept. of Statistics, Techn. Rep, 3.
de Luna, X. and Genton, M. (2000). Robust simulation-based estimation. Statist. Probab. Lett., 48(3),
253–259.
de Luna, X. and Genton, M. (2001). Robust simulation-based estimation of ARMA models. J.
Comput. Graph. Stat., 10(2), 370–387.
Denby, L. and Martin, R. (1979). Robust estimation of the first-order autoregressive parameter. J.
Amer. Statis. Assoc., 74(365), 140–146.
do Rêgo Sousa, T., Haug, S. and Klüppelberg, C. (2017). Indirect Inference for Lévy-driven
continuous-time GARCH models. ArXiv https://arxiv.org/abs/1712.09870.
Fasen, V. and Kimmig, S. (2017). Information criteria for multivariate CARMA processes. Bernoulli,
23, 2860–2886.
Fasen-Hartmann, V. and Scholz, M. (2017). Quasi-maximum likelihood estimation for cointegrated solutions of continuous-time state space models observed at discrete time points. ArXiv
https://arxiv.org/abs/1712.08665.
Fuchs, F. and Stelzer, R. (2013). Mixing conditions for multivariate infinitely divisible processes
with an application to mixed moving averages and the supOU stochastic volatility model. ESAIM:
Probability and Statistics, 17, 455–471.
Genton, M. and Lucas, A. (2003). Comprehensive definitions of breakdown points for independent
and dependent observations. J.R. Stat. Soc. Ser. B Stat. Methodol., 65(1), 81–94.
Gouriéroux, C. and Monfort, A. (1997). Simulation-based Econometric Methods. Core Lectures.
Oxford University Press, Oxford.
Gouriéroux, C., Monfort, A. and Renault, E. (1993). Indirect inference. J. Appl. Econometr., 8,
85–118.
Hampel, F. (1971). A general qualitative definition of robustness. Ann. Math. Statist., pages 1887–
1896.
Hampel, F. (1974). The influence curve and its role in robust estimation. J. Amer. Statist. Assoc.,
69(346), 383–393.
Hannan, E. and Deistler, M. (2012). The Statistical Theory of Linear Systems. Society for Industrial
and Applied Mathematics. Reprint of the 1988 original.
Huber, P. and Ronchetti, E. (2009). Robust Statistics. 2nd edn. Wiley.
Hutton, J. and Nelson, P. (1984). Interchanging the order of differentiation and stochastic integration.
Stochastic Process. Appl., 18(2), 371–377.
Ibragimov, I. (1962). Some limit theorems for stationary processes. Theory Probab. Appl., 7(4),
349–382.

39

Jiang, W. and Turnbull, B. (2004). The indirect method: inference based on intermediate statistics – a
synthesis and examples. Statist. Sci., 19(2), 239–263.
Kimmig, S. (2016). Statistical Inference for MCARMA Processes. Ph.D. thesis, Karlsruher Institut of
Technology, Karlsruhe.
Künsch, H. (1984). Infinitesimal robustness for autoregressive processes. Ann. Statist., 12(3), 843–
863.
Maronna, R., Bustos, O. and Yohai, V. (1979). Bias-and efficiency-robustness of general M-estimators
for regression with random carriers. In: Smoothing techniques for curve estimation, pages 91–116.
Springer.
Maronna, R., Martin, D. and Yohai, V. (2006). Robust Statistics. Wiley.
Maronna, R. and Yohai, V. (1981). Asymptotic behavior of general M-estimates for regression and
scale with random carriers. Z. Wahrscheinlichkeit, 58(1), 7–20.
Maronna, R. and Yohai, V. (1991). The breakdown point of simultaneous general M estimates of
regression and scale. J. Amer. Statist. Assoc., 86(415), 699–703.
Marquardt, T. and Stelzer, R. (2007). Multivariate CARMA processes. Stochastic Process. Appl.,
117(1), 96–120.
Martin, D. (1980). Robust estimation of autoregressive models. In: D. Brillinger and G. Tiao (Eds.),
Directions in time series, vol. 1, pages 228–262. Institute of Mathematical Statistics, Hayward, CA.
Martin, D. and Jong, J. (1977). Asymptotic properties of robust generalized M-estimates for the first
order autoregressive parameter. Unpublished Bell Labs Technical Memo, Murray Hill, NJ.
Martin, D. and Yohai, V. (1985). Robustness in time series and estimating ARMA models. In: E. Hannan, P. Krishnaiah and M. Rao (Eds.), Handbook of statistics, vol. 5, pages 119–155. Elsevier.
Martin, D. and Yohai, V. (1986). Influence functionals for time series. Ann. Statist., 14(3), 781–818.
Muler, N., Peña, D. and Yohai, V. (2009). Robust estimation for ARMA models. Ann. Statist., 37(2),
816–840.
Olive, D. (2017). Robust Multivariate Analysis. Springer, Cham.
Pollard, D. (1990). Empirical processes: theory and applications, vol. 2. Institute of Mathematical
Statistics, Hayward, CA.
Schilling, R. and Partzsch, L. (2014). Brownian motion. 2nd edn. De Gruyter, Berlin.
Schlemm, E. and Stelzer, R. (2012). Quasi maximum likelihood estimation for strongly mixing state
space models and multivariate Lévy-driven CARMA processes. Electron. J. Stat., 6, 2185–2234.
Smith, A. (1993). Estimating nonlinear time-series models using simulated vector autoregressions. J.
Appl. Econometr., 8(S1), S63–S84.
Wilcox, R. (1967). Exponential operators and parameter differentiation in quantum physics. J. Math.
Physics, 8, 962–982.

40

