Maximum Likelihood Estimation for Wishart processes
Aurélien Alfonsi1 , Ahmed Kebaier2 , and Clément Rey1

arXiv:1508.03323v2 [] 15 Apr 2016

1

CERMICS, École des Ponts, UPE, Projet MathRisk ENPC-INRIA-UMLV,
Champs-sur-Marne, France
2
Université Paris 13, Sorbonne Paris Cité, LAGA, CNRS (UMR 7539)
April 18, 2016

Abstract
In the last decade, there has been a growing interest to use Wishart processes for
modelling, especially for financial applications. However, there are still few studies on
the estimation of its parameters. Here, we study the Maximum Likelihood Estimator
(MLE) in order to estimate the drift parameters of a Wishart process. We obtain precise
convergence rates and limits for this estimator in the ergodic case and in some nonergodic
cases. We check that the MLE achieves the optimal convergence rate in each case.
Motivated by this study, we also present new results on the Laplace transform that
extend the recent findings of Gnoatto and Grasselli [17] and are of independent interest.

Keywords : Wishart processes, Laplace transform, parameter inference, maximum likelihood, limit theorems, local asymptotic properties.
AMS MSC 2010: 62F12, 44A10, 60F05, 91B70.

1

Introduction and preliminary results

The goal of this paper is to study the maximum likelihood estimation of the parameters of
Wishart processes. These processes have been introduced by Bru [7] and take values in the
set of positive semidefinite matrices. Let d ∈ N∗ denote the dimension, Md be the set of
real d-square matrices, Sd+ (resp. Sd+,∗ ) be the subset of positive semidefinite (resp. definite)
matrices, Sd (resp. Ad ) the subset of symmetric (resp. antisymmetric) matrices. Wishart
processes are defined by the following SDE
√
√



dXt = αa> a + bXt + Xt b> dt + X t dWt a + a> dWt> X t , t > 0
(1)
X0 = x ∈ Sd+ ,
where α > d−1, a ∈ Md , b ∈ Md and (Wt )t>0 denotes a d-square matrix made of independent
√
Brownian motions. We recall that for x ∈ Sd+ , x is the unique matrix in Sd+ such that
e-mails : alfonsi@cermics.enpc.fr, kebaier@math.univ-paris13.fr, reyc@cermics.enpc.fr. This research
benefited from the support of the “Chaire Risques Financiers”, Fondation du Risque and the Laboratories of
Excellence Bézout http://bezout.univ-paris-est.fr/ and MME-DII http://labex-mme-dii.u-cergy.fr/ .

1

1 INTRODUCTION AND PRELIMINARY RESULTS
√

2

2

x = x. It is shown by Bru [7] and Cuchiero et al. [8] in a more general affine setting that
the SDE (1) has a unique strong solution when α > d + 1 and a unique weak solution when
α > d − 1. Besides, we have Xt ∈ Sd+,∗ for any t > 0 when x ∈ Sd+,∗ and α > d + 1. In this
paper, we will denote by W ISd (x, α, b, a) the law of (Xt , t > 0) and W ISd (x, α, b, a; t) the
law of Xt . In dimension d = 1, Wishart processes are known as Cox-Ingersoll-Ross processes
in the literature. It is worth recalling that the law of X only depends on a through a> a since
we have
√
W ISd (x, α, b, a) = W ISd (x, α, b, a> a),
law

see e.g. equation (12) in [1]. Therefore, the parameters to estimate are α, b and a> a.
Wishart processes have been originally considered by Bru [6] to model some biological
data. Recently, they have been widely used in financial models in order to describe the
evolution of the dependence between assets. Namely, Gourieroux and Sufana [19] and Da
Fonseca et al. [10] have proposed a stochastic volatility model for a basket of assets that
assumes that the instantaneous covariance between the assets follows a Wishart process.
This extends the well-known Heston model [21] to many assets. Wishart processes have also
been used for interest rates models. Affine term structure models involving these processes
have been proposed for example by Gourieroux and Sufana [20], Gnoatto [16] and Ahdida
et al. [2]. For these models, the question of estimating the parameters of the underlying
Wishart process may be important for practical purposes and should be possible thanks to
the profusion of financial data. This issue has been considered by Da Fonseca et al. [9] for the
model presented in [10]. However, there is no dedicated study on the Maximum Likelihood
Estimator (MLE) for Wishart processes. For the Cox-Ingersoll-Ross process, the estimation
of parameters has been studied earlier, motivated in particular by its use for interest rates
(see Fournié and Talay [14]). Later on, the MLE has been studied by Overbeck [31] including
some nonergodic cases, and more recently by Ben Alaya and Kebaier [4, 5]. This paper
completes the literature by studying the MLE for Wishart processes.
In this paper, we will follow the theory developed in the books by Lipster and Shiryaev [27]
and Kutoyants [23] and assume that we observe the full path (Xt , t ∈ [0, T ]) up to time T > 0.
This choice will be convenient from a mathematical point of view to study the convergence
of the MLE. Of course, in practice it can be relevant to study precisely the estimation when
we only observe the process on a discrete time-grid. This is left for further research, but
we already observe in our numerical experiments that the discrete approximation of the
MLE gives a satisfactory estimation of Wishart parameters (see Section 6). It is worth
noticing that once we observe the path (Xt , t ∈ [0, T ]), the parameter a> a is known. In
fact, we can calculate the quadratic covariation (see for example Lemma 2 in [1]) and get for
i, j, k, l ∈ {1, · · · , d}
Z
hXi,j , Xk,l iT =

T

(a> a)j,l (Xs )i,k + (a> a)j,k (Xs )i,l + (a> a)i,l (Xs )j,k + (a> a)i,k (Xs )j,l ds.

0

(2)
This leads to
>

(a a)i,i
(a> a)i,j

Z T
−1
1
= hXi,i iT
,
(Xs )i,i ds
4
0

Z T
Z T
−1
1
=
hXi,j , Xi,i iT − (a> a)i,i
(Xs )i,j ds
(Xs )i,i ds
,
2
0
0

(3)

1 INTRODUCTION AND PRELIMINARY RESULTS

3

for 1 6 i, j 6 d and j 6= i. We note that these quantities are well defined as soon as the path
(Xt , t ∈ [0, T ]) has a finite quadratic variation and is such that Xt ∈ Sd+,∗ dt-a.e., which is
satisfied by the paths of Wishart processes (see Proposition 4 in [7]). We will assume that
a> a ∈ Sd+,∗ and denote by a ∈ Md an invertible matrix that matches the observed value of
a> a: a can be for example the square root of a> a or the Cholesky decomposition of a> a.
Then, we know that Yt = (a> )−1 Xt a−1 follows the law W ISd ((a> )−1 xa−1 , α, (a> )−1 ba> , Id ),
see e.g. equation (13) in [1]. It is therefore sufficient to focus on the estimation of the
parameters α and b when a = Id , which we consider now.
We first present the MLE of θ = (b, α), and we denote by Pθ the original probability
measure under which X satisfies
h
i
√
√
dXt = αId + bXt + Xt b> dt + X t dWt + dWt> X t .
(4)
When no confusion is possible, we also denote P this probability. We consider α0 > d + 1
and set θ0 = (α0 , 0). We will assume for the joint estimation of α and b that
α > d + 1 and x ∈ Sd+,∗ .

(5)

The latter assumption is not restrictive in practice since the condition α > d + 1 ensures that
Xt ∈ Sd+,∗ for any t > 0. Due to this assumption, we know by Theorem 4.1 in Mayerhofer [29]
that
Z T

Z
p
dPθ0 ,T
1 T
α0 − α p
>
:= exp
Tr[Hs dWs ] −
Tr[Hs Hs ]ds , with Ht =
( Xt )−1 − b Xt
dPθ,T
2 0
2
0
Rt
defines a probability measure under which W̃t = Wt − 0 Hs> ds is a d × d-Brownian motion,
where Pθ,T is the restriction of Pθ to the σ-algebra σ(Ws , s ∈ [0, T ]). We have
√
√
dXt = α0 Id dt + X t dW̃t + dW̃t> X t ,
and the likelihood is then defined by (see Lipster and Shiryaev [27], Chapter 7)
0
Lθ,θ
=
T

1

,

R
RT
T
E exp 0 Tr[Hs dWs ] − 12 0 Tr[Hs Hs> ]ds FTX

(6)

where (FtX )t>0 denote the filtration generated by the process X.
Proposition 1.1. For X ∈ Sd+,∗ , let LX : Sd → Sd be the linear application defined by
LX (Y ) = XY + Y X. It is invertible, and the likelihood of (Xt , t ∈ [0, T ]) is given by


Z
α − α
 T
det[XT ]
α − α0 α + α0
αT
0
Tr[Xs−1 ]ds −
= exp
log
−
−1−d
Tr[b]
4
det[x]
4
2
2
0
Z
h


i 1Z T
h


i 
1 T
−1
>
>
>
+
Tr LXt bXt + Xt b dXt −
Tr L−1
bX
+
X
b
(bX
+
X
b
)
dt .
t
t
t
t
Xt
2 0
4 0
(7)

0
Lθ,θ
T

Lemmas B.1 and B.2 states some properties of LX , and the proof of Proposition 1.1 is given
dP 0 ,T
∈ FTX if, and only if b ∈ Sd ,
in Appendix A. In particular, we see from this proof that dPθθ,T

1 INTRODUCTION AND PRELIMINARY RESULTS

4

in which case the likelihood has the following simpler form
0
Lθ,θ
T



Z
α − α
Tr[bXT ] − Tr[bx] 1 T
det[XT ]
0
= exp
+
Tr[b2 Xs ]ds
log
−
4
det[x]
2
2 0
Z

 T
αT
α − α0 α + α0
−
−1−d
Tr[Xs−1 ]ds −
Tr[b] , (8)
4
2
2
0

since L−1
Xt (bXt + Xt b) = b.
Now, we want to maximize the likelihood and observe that the quantity in the exponential (7) is quadratic with respect to (b, α) and goes almost surely
to −∞ when k(b, α)k → +∞.

> ] by Lemma B.1. Then, CauchyTo do so, we first remark that Tr[b] = Tr[L−1
bX
+
X
b
t
t
Xt
Schwarz inequality yields to


Z
q
p
√ −1 
α
1 T
−1
>
Tr
2LXs bXs + Xs b
Xs √
Xs ds
| Tr[αb]| =
T 0
2
Z
Z
h


i


1 T
α2 1 T
−1
> 2
Tr (LXs bXs + Xs b ) Xs ds +
Tr Xs−1 ds
(9)
6
T 0
4 T 0
Z T
Z

i

h


α2 1 T
1
−1
>
>
=
Tr LXs bXs + Xs b (bXs + Xs b ) ds +
Tr Xs−1 ds,
2T 0
4 T 0
and it is strict almost surely, which gives that the quadratic form in the exponential (7) is
negative definite. There is thus a unique global maximum of (7) on R × Md . We know from
Lemma B.2 that L−1
Xs is self-adjoint, and we get with straightforward calculations that the
MLE θ̂T = (b̂T , α̂T ) is characterized by the following equations:


(
RT
det[XT ]
1
T
−1
− α̂T −1−d
log
4
4
0 Tr[Xs ]ds − 2 Tr[b̂T ] = 0,
det[x]
(10)
R T −1
R T −1
α̂T T
>
0 LXs (dXs )Xs − 0 LXs (b̂T Xs + Xs b̂T )Xs ds − 2 Id = 0.
Unless in the ergodic case, we will not be able to obtain convergence results for this estimator.
Instead, we will mostly work with the MLE estimator when b is known to be symmetric. This
enables us to work with more tractable formulas, even if the calculations are already quite
involved in case. Analyzing the general case would require development of further arguments.
Besides, we can consider that Wishart processes with b symmetric already form an interesting
family of processes that may be rich enough in many applications. When b ∈ Sd , the unique
global maximum θ̂T = (b̂T , α̂T ) of (8) on R × Sd is characterized by the following equations:


(
RT
det[XT ]
1
T
−1
log
− α̂T −1−d
4
4
0 Tr[Xs ]ds − 2 Tr[b̂T ] = 0,
det[x]
(11)
R
XT −x
α̂T T
1 T
−
(
b̂
X
+
X
b̂
)ds
−
I
=
0.
s
s
T
T
d
2
2 0
2
To get more explicit formulas, we have to invert this linear system. For X ∈ Sd and a ∈ R,
we define the linear applications
LX : Sd → Sd
Y 7→ Y X + XY

and LX,a : Sd → Sd
Y 7→ Y X + XY − 2a Tr[Y ]Id .

We introduce the following shorthand notation
Z T
−1
Z T
RT :=
Xs ds, QT :=
Tr[Xs−1 ]ds
,
0

0


ZT := log

det[XT ]
det[x]

(12)


,

(13)

1 INTRODUCTION AND PRELIMINARY RESULTS

5

and note that QT and ZT are defined only for α > d + 1 while RT is defined for α > d − 1
and belongs almost surely to Sd+,∗ .1 By using the convexity property of the inverse, see e.g.
Mond and Pecaric [30], we have when α > d + 1
"
 #
Q−1
RT −1
Tr
< T , a.s.
(14)
T
T

We get α̂T = 1 + d + QT ZT − 2T Tr[b̂T ] and LRT ,T 2 QT (b̂T ) = XT − x − T (QT ZT + 1 + d) Id .
By (14) and Lemma B.1, the latter equation can be inverted, which leads to


 −1


α̂
=
1
+
d
+
Q
Z
−
2T
Tr
L
(X
−
x
−
T
[Q
Z
+
1
+
d]
I
)
T
T
T
T T
d
 T
RT ,T 2 QT
(15)


−1
b̂T = LRT ,T 2 QT (XT − x − T [QT ZT + 1 + d] Id ) .
The estimator of α when α ∈ [d − 1, d + 1) given by the MLE is no longer well defined.
The same thing already occurs in dimension d = 1 for the CIR process, see Ben Alaya and
Kebaier [4]. However, it is still possible to estimate the parameter b ∈ Md when α > d − 1
is known. In this case, we denote θ = (b, α) and θ0 = (0, α) and get by repeating the same
arguments that

i

h
1 Z T
−1
>
0
dX
bX
+
X
b
Tr
L
Lθ,θ
=
exp
t
t
t
Xt
T
2 0
Z T


i
h

1
αT
>
>
−
bX
+
X
b
(bX
+
X
b
)
dt
−
Tr L−1
Tr[b]
,
t
t
t
t
Xt
4 0
2
and the MLE is characterized by
Z
Z T
−1
LXs (dXs )Xs −
0

0

T

>
L−1
Xs (b̂T Xs + Xs b̂T )Xs ds −

αT
Id = 0.
2

(16)

When b is known a priori to be symmetric, the likelihood and the MLE are then given by
 Tr[bX ] − Tr[bx] 1 Z T

αT
T
θ,θ0
LT = exp
−
Tr[b2 Xs ]ds −
Tr[b] ,
(17)
2
2 0
2
b̂T = L−1
RT (XT − x − αT Id ) .

(18)

The goal of the paper is to study the convergence of the MLE under the original probability Pθ . To do so, we first consider the case where the Wishart process is ergodic. By
Lemma C.1, this holds if −(b + b> ) ∈ Sd+,∗ when b ∈ Md , and the ergodicity is equivalent
to −b ∈ Sd+,∗ when b ∈ Sd . Then, we can use Birkhoff’s ergodic theorem to determine the
convergence of the MLE. Section 2 presents these results for (15) when α > d+1, for (7) when
α > d + 1 and for both (18) and (16) when α > d − 1. Section 3 studies the convergence of
the MLE in some nonergodic cases, namely when b = λ0 Id with λ0 > 0 and when b is known
to be symmetric. More precisely, when b = 0, we obtain convergence results for (15) when
α > d + 1 and for (18) when α > d − 1. When λ0 > 0, we only obtain convergence results
1
This is obvious when α > d − 1 since Xt ∈ Sd+,∗ a.s. by Proposition 4 in [7]. For α = d − 1, we would have
by contradiction the existence of vT ∈ FTX such that ∀t ∈ [0, T ], vT> Xt vT = 0. This is clearly not possible by
using the connection with matrix-valued Ornstein-Uhlenbeck in this case, see eq. (5.7) in [7].

2 STATISTICAL INFERENCE: THE ERGODIC CASE

6

for (18) when α > d−1. In all these cases, we analyse the convergence by the mean of Laplace
transforms. Though limited to some nonergodic cases, we however recover and extend the
recent convergence results obtained by Ben Alaya and Kebaier for the one-dimensional CIR
process [4, 5]. In Section 4, we check that the MLE achieves the optimal rate of convergence
in the different cases by proving local asymptotic properties. Last, we study in Section 5 the
Laplace transform of (XT , RT ). This study can be of independent interest and improves the
recent results of Gnoatto and Grasselli [17].

2

Statistical Inference of the Wishart process: the ergodic case

When −(b + b> ) ∈ Sd+,∗ , the Wishart process Xt converges in law when t → +∞ to the
R∞
√
>
stationary law X∞ ∼ W ISd (0, α, 0, 2q∞ ; 1/2) with q∞ = 0 esb esb ds for any starting
point x ∈ Sd+ by Lemma C.1. Therefore this is the unique stationary law which is thus
extremal, and we know by Stroock ([35], Theorem 7.4.8) that it is then ergodic, see also
Pagès [32], Annex A. We introduce the following quantity
R∞ := Eθ (X∞ ).
From the ergodic Birkhoff’s theorem, we have
RT a.s.
−→ R∞ ,
T
Besides, when α > d + 1, Q∞ =

1
−1
Eθ (Tr[X∞
])

as T → +∞.

(19)

is finite and satisfies
−1

Q∞ Tr[R∞ ] < 1,

(20)

due to the convexity property of the inverse, see e.g. Mond and Pecaric [30]. Again, the
ergodic Birkhoff’s theorem gives
a.s.

T QT −→ Q∞ =

1
−1 ,
Eθ (Tr[X∞
])

as T → +∞.

(21)

This section is organized as follows. First, we study the MLE (15) when b is known to be
symmetric in the cases α > d + 1 and α = d + 1. Then, we focus on the MLE (10) when
b ∈ Md and α > d + 1. The analysis follows the same steps and reuses some calculations
made in the symmetric case. Last, we study the convergence of the MLE when α > d − 1 is
known, in both symmetric and general cases.

2.1

The global MLE estimator of θ = (b, α) when b is known to be symmetric

When b ∈ Sd , the ergodicity is by Lemma (C.1) equivalent
to −b ∈ Sd+,∗ , which we assume
√
in this subsection. We have X∞ ∼ W ISd (0, α, 0, −b−1 ; 1/2) and it is easy to get from (4)
that αId + bR∞ + R∞ b = 0, which gives R∞ = − α2 b−1 ∈ Sd+,∗ . We will also show in the
proof of Theorem 2.1 that
α − (1 + d)
Q∞ :=
.
(22)
2 Tr[−b]

2 STATISTICAL INFERENCE: THE ERGODIC CASE

7

We consider the convergence of the MLE given by (15) when α > d + 1. We introduce
the following martingales:
Z tp
Z t
√
dWs> X s ,
Mt :=
Xs dWs +
(23)
0
0
Z t
p
Tr[( Xs )−1 dWs ].
Nt :=
(24)
0

We use the dynamics of (Xt )t>0 under Pθ and Itô’s formula for (Zt )t>0 (see e.g. Bru [7],
equation (2.6)) to get on the one hand
ZT = (α − 1 − d)Q−1
T + 2 Tr[b]T + 2NT .

XT = x + αT Id + LRT (b) + MT ,

(25)

On the other hand, we obtain from (11) and (13) that XT = x + α̂T T Id + LRT (b̂T ) and
ZT = (α̂T − 1 − d)Q−1
T + 2T Tr[b̂T ], which yields to




α̂T − α = 2T QT Tr[b − b̂T ] + 2QT NT
(26)
LRT (b̂T − b) = (α − α̂T )T Id + MT = 2T 2 QT Tr[b̂T − b]Id + MT − 2T QT NT Id .


√
T (b̂T − b, α̂T − α)
Theorem 2.1. Assume that −b ∈ Sd+,∗ and α > d + 1. Under Pθ ,
converges in law when T → +∞ to the centered Gaussian vector (G, H) that takes values in
Sd × R and has the following Laplace transform: for c, λ ∈ Sd × R,
2Q∞ λ2

Eθ [exp (Tr[cG] + λH)] = exp

−1

1 − Q∞ Tr[R∞ ]

2Q∞ λ

−
1−

−1
−1 Tr[cR∞ ]
Q∞ Tr[R∞ ]

!
+

Tr[cL−1
(c)]
R∞ ,Q∞

Proof. By (14) and Lemma B.1, we can rewrite the system (26) as follows




√
N
M
N

T
T
T
−1
 T (α̂T − α) = 2T QT √ − 2T QT Tr L R
√ − 2T QT Id √


T ,T Q

T
T
T
T
T


√


T (b̂T − b) = L−1

RT
T


,T QT

MT
N
√ − 2T QT Id √ T
T
T


.

Note that, for i, j, k, l ∈ {1, . . . , d} we have
hMi,j , Mk,l it = [δjl (Rt )i,k + δjk (Rt )i,l + δil (Rt )j,k + δik (Rt )j,l ] ,
hMi,j , N it = 2tδij

and

hN it = Q−1
t ,

(27)

where δij stands for the Kronecker symbol.
So, it follows from the central limit theorem for martingales (see e.g., Kutoyants [23],
NT
√T , √
Proposition 1.21), that ( M
) converges in law under Pθ towards a centered Gaussian
T
T
vector (G̃, H̃) taking values in Sd × R such that


Eθ (G̃i,j G̃k,l ) = δjl (R∞ )i,k + δjk (R∞ )i,l + δil (R∞ )j,k + δik (R∞ )j,l ,
−1

Eθ (G̃i,j H̃) = 2δi,j and Eθ (H̃ 2 ) = Q∞ .

(28)

.

2 STATISTICAL INFERENCE: THE ERGODIC CASE

8

From (25) and (21), we obtain (22). From Lemma B.1, the function (X, Y, a) 7→ L−1
X,a (Y ) is
√
√
continuous, and we get by Slutsky’s theorem that ( T (b̂T − b), T (α̂T − α)) converges in
law to the Gaussian vector




h

i
−1
(G, H) = L−1
G̃
−
2Q
G̃
−
2Q
H̃I
H̃
−
Tr
L
H̃I
,
2Q
.
d
d
∞
∞
∞
R ,Q
R ,Q
∞

∞

∞

∞

We are interested to calculate the Laplace transform of this law. First, we calculate the
Laplace transform of (G̃, H̃):
 

h

i
1 2 −1
2
λ Q∞ + 4λ Tr[c] + 4 Tr[c R∞ ] .
∀c ∈ Sd , λ ∈ R, Eθ exp Tr[cG̃] + λH̃ = exp
2
(29)
We want to calculate for c ∈ Sd and λ ∈ R,
h

i
Eθ [exp (Tr[cG] + λH)] = Eθ exp Tr[(c − 2λQ∞ Id )G] + 2λQ∞ H̃ .
Due to (20) and Lemma B.1, we can introduce c̃ = L−1
R

∞ ,Q∞

(c − 2λQ∞ Id ). We have

R∞ c̃ + c̃R∞ − 2Q∞ Tr[c̃]Id = c − 2λQ∞ Id ,
and thus
Tr[(c − 2λQ∞ Id )G] = Tr[(R∞ c̃ + c̃R∞ − 2Q∞ Tr[c̃]Id )G]
= Tr[c̃(R∞ G + GR∞ − 2Q∞ Tr[G]Id )] = Tr[c̃(G̃ − 2Q∞ H̃Id )].
We therefore obtain from (29)
Eθ [exp (Tr[cG] + λH)]
h

i
= Eθ exp Tr[c̃(G̃ − 2Q∞ H̃Id )] + 2λQ∞ H̃
i
h

= Eθ exp Tr[c̃G̃] + 2Q∞ (λ − Tr[c̃])H̃


= exp 2 (λ − Tr[c̃])2 Q∞ + 2(λ − Tr[c̃]) Tr[c̃]Q∞ + Tr[c̃2 R∞ ] .
Since 2 Tr[c̃2 R∞ ] = Tr[c̃(c̃R∞ + R∞ c̃)] = Tr[c̃c] + 2Q∞ (Tr[c̃] − λ) Tr[c̃], we get

Eθ [exp (Tr[cG] + λH)] = exp 2λ(λ − Tr[c̃])Q∞ + Tr[c̃c] .
We now use that L−1
R

(Id ) =

Since we have Tr[L−1
R

(c)] =

∞ ,Q∞

∞ ,Q∞

result.

−1
1
−1 R
2(1−Q∞ Tr[R∞ ]) ∞
−1
Tr[R∞ c]
by
−1
2(1−Q∞ Tr[R∞ ])

to get c̃ = L−1
R

∞ ,Q∞

−1

(c) − λ

Q∞ R ∞

−1

1−Q∞ Tr[R∞ ]

.

Lemma B.1, this yields to the claimed

When α = d + 1, the rate of convergence of the MLE of α is even better as stated by the
following theorem.

√
T (b̂T − b), T (α̂T − α)
Theorem 2.2. Assume −b ∈ Sd+,∗ and α = d+1. Then, under Pθ ,


converges in law when T → +∞ to G, −2τ−−1Tr[b] Tr[b] , where τa = inf{t ≥ 0, Bt = a} with
(Bt )t≥0 a given one-dimensional standard Brownian
 motion and
 G is a Gaussian vector
−1
(c)] , c ∈ Sd .
independent of B such that Eθ [exp (Tr[cG])] = exp Tr[cLR
∞

2 STATISTICAL INFERENCE: THE ERGODIC CASE

9

Proof. By (14) and Lemma B.1, we can rewrite the system (26) as follows





NT
MT
NT

−1
1
2
3/2

√
√ − 2T QT Id
T (α̂T − α) = 2T QT
− T Tr L RT


,T QT

T
T
T
T


√


 T (b̂T − b) = L−1
RT
T



MT
N
√ − 2T 3/2 QT Id T
T
T

,T QT

(30)


.

From (25), we have
NT
1
=
log
T
2T



det[XT ]
det[x]


− Tr[b].

As for −b ∈ Sd+,∗ the Wishart process (Xt )t≥0 is stationary with invariant limit distribution
X∞ we easily deduce that NTT converges in probability to − Tr[b] when T → ∞. Then, it
follows from (19) that
(T −1 RT , T −1 NT ) →θ (R∞ , − Tr[b]),
P

as T → ∞.

(31)

Hence, we only need to study the asymptotic behavior of the couple (T −1/2 MT , T 2 QT ).
According to Theorem 4.1 in Mayerhofer [29], we have for λ > 0 and Γ ∈ Sd



Z
1
2 T
λ2 −1
2λ
λ
2
= 1.
NT −
Q + √ Tr[ΓMT ] −
Tr[Γ Xs ]ds − √ Tr[Γ]
Eθ exp
T
2T 2 T
T 0
T
T

(32)

Now, let us introduce the quantity




NT
1
λ2 −1
√
AT = Eθ exp λ
+ λ Tr[b] exp − 2 QT +
Tr[ΓMT ]
T
2T
T


Z
2 T
× exp −
Tr[Γ2 Xs ]ds + 2 Tr[Γ2 R∞ ] .
T 0


Then, by (32) we easily get AT = exp λ Tr[b] + 2 Tr[Γ2 R∞ ] + √2λT Tr[Γ] . We now write
h

i
−1
λ2
1
√
+
AT = ÃT + Eθ exp − 2T
Q
Tr[ΓM
]
with
T
2
T
T




1
λ2 −1
ÃT = Eθ (exp (ξT ) − 1) exp − 2 QT + √ Tr[ΓMT ]
2T
T
Z T
NT
2
ξT = λ
+ λ Tr[b] −
Tr[Γ2 Xs ]ds + 2 Tr[Γ2 R∞ ].
T
T 0


Cauchy-Schwarz inequality and Q−1
T > 0 give
|ÃT | 6

1/2
Eθ [exp (2ξT )

− 2 exp (ξT ) +

1/2
1]Eθ




2
exp √ Tr[ΓMT ] .
T

On the one hand, Proposition 5.1 with m = −b ∈ Sd+,∗ gives






2
2
2
6 Eθ exp
Tr[Γ RT ]
< ∞.
Eθ exp √ Tr[ΓMT ]
T
T

2 STATISTICAL INFERENCE: THE ERGODIC CASE

10

On the other hand, we have for any r > 0,



λr
Eθ [exp (rξT )] 6 Eθ exp
exp(2r Tr[Γ2 R∞ ]).
NT
T
From (25), we have
"

 λr #


det[XT ] 2T
λr
= exp(−λr Tr[b])Eθ
NT
.
Eθ exp
T
det[x]
The sublinear growth of the
of the Wishart SDE and the convergence to a sta coefficients
λ̃ 
T]
tionary law gives that Eθ det[X
is uniformly bounded in T > 0, λ̃ < 1 and theredet[x]


 λr
T ] 2T
< ∞. This gives the uniform integrability of the family
fore supT > λr Eθ det[X
det[x]
2

(exp (2ξT ) , T > λ). Then, we deduce from (31) that Eθ [exp (2ξT ) − 2 exp (ξT ) + 1]
and thus ÃT

→

T →+∞

→

T →+∞

0

0.

Hence, we obtain




λ2 −1
1
lim Eθ exp − 2 QT + √ Tr[ΓMT ]
= lim AT = exp λ Tr[b] + 2 Tr[Γ2 R∞ ] .
T →∞
T →∞
2T
T
Therefore, we deduce by Lemma B.4 the following convergence in law
!


q
q
Q−1
M
T
>
T
,√
⇒ τ− Tr[b] , R∞ G̃ + G̃
R∞ as T → ∞,
T2
T
where G̃i,j 1 6 i, j 6 d are independent standard normal variables. Together with (31), we
obtain that
q
q
−1
2
−1
−1/2
>
(T RT , T QT , T NT , T
MT ) ⇒ (R∞ , 1/τ− Tr[b] , − Tr[b], R∞ G̃ + G̃
R∞ ), (33)
which gives the claim by (30) and Lemma B.4.

2.2

The global MLE estimator of θ = (b, α) when b ∈ Md

We define the linear operators L̄X , L̄X,a : Md → Md by
>
L̄X (Y ) = L−1
X (Y X + XY )X, L̄X,a (Y ) = L̄X (Y ) − a Tr[Y ]Id .

From (4), we get ZT = (α − 1 − d)Q−1
T + 2 Tr[b]T + 2NT . This yields with (10) to



α̂T − α = 2T QT Tr[b − b̂T ] + 2QT NT
(34)

 RT
0

L̄Xs (b̂T − b)ds −

T 2 QT

Tr[b̂T − b]Id =

We now define
1
L̂T (Y ) =
T

Z

RT

−1
0 LXs (dMs )Xs

T

L̄Xs (Y )ds − T QT Tr[Y ]Id ,
0

− T QT NT Id .

2 STATISTICAL INFERENCE: THE ERGODIC CASE

11

which is a linear operator on Md . By using the convexity of the inverse
function, there
γ RT
1
exists γ ∈ (0, 1) that depends on (Xs , s ∈ [0, T ]) such that T QT = T 0 Tr[X −1 ] ds. We get
s
RT
L̂T (Y ) = T1 0 L̄Xs , γ−1 (Y )ds. By Lemma B.3, L̂T is self adjoint and positive. It is even
Tr[Xs

]

positive definite since Tr[L̂T (Y )> Y ] = 0 implies by using Lemmas B.3 that Y Xs + Xs Y > = 0
a.s. on [0, T ] under Pθ , and therefore the quadratic variation of Tr[uY Xs ] is equal to zero for
RT
any u ∈ Sd . This gives 0 Tr[(uY + Y > u)Xs (uY + Y > u)]ds = 0 and thus uY + Y > u = 0 for
all u ∈ Sd , which necessarily implies Y = 0. Thus, we rewrite (34) as
 √
√
NT
 T (α̂T − α) = −2T QT Tr[ T (b̂T − b)] + 2T QT √
T


√
R T −1
(35)
NT
√1
√

T (b̂T − b) = L̂−1
L
I
(dM
)X
−
T
Q
.
s
s
T
d
T
Xs
T 0
T
We will assume −(b + b> ) ∈ Sd+,∗ and know from Lemma C.1 that XT converges in law under
√
Pθ to the stationary law X∞ ∼ W ISd (0, α, 0, 2q∞ ; 1/2). We define
L̂∞ (Y ) = Eθ [L̄X∞ (Y )] − Q∞ Tr[Y ]Id .

(36)

Note that for Y ∈ Sd , L̂∞ (Y ) = Y R∞ − Q∞ Tr[Y ]Id . From the convexity of the inverse func−1 ]] with γ ∈ (0, 1), and thus L̂ (Y ) = E [L̄
γ
tion, Q∞ = γEθ [1/ Tr[X∞
(Y )] is a self∞
θ X∞ ,
−1
Tr[X∞ ]

adjoint positive operator by Lemma B.3. It is even positive definite since Tr[L̂∞ (Y )> Y ] = 0
implies by Lemma B.3 that Y X∞ + X∞ Y > = 0 almost surely. Since the law X∞ has a
positive density on Sd+,∗ , this gives Y u + uY > = 0 for any u ∈ Sd and thus Y = 0.
√

Theorem 2.3. Assume −(b+b> ) ∈ Sd+,∗ and α > d+1. Then, under Pθ ,
T (b̂T − b), T (α̂T − α)
converges in law when T → +∞ to the centered Gaussian vector (G, H) that takes values in
Md × R and has the following Laplace transform: for c, λ ∈ Md × R,

2Q∞ λ2
2Q∞ λ
−1
>
Eθ [exp(Tr[c G] + λH)] = exp
(37)
−1 −
−1 Tr[cR∞ ]
1 − Q∞ Tr[R∞ ] 1 − Q∞ Tr[R∞ ]
1
−1
−1
>
−1
−1
>
+ Eθ [Tr[L−1
X∞ (L̂∞ (c)X∞ + X∞ L̂∞ (c) )(L̂∞ (c)X∞ + X∞ L̂∞ (c) )]]
4
!
−1
−1
1
Tr[cR∞ ]Q∞
Tr[cR∞ ]
−1
+
−1
−1 − Tr[L̂∞ (c)]
1 − Q Tr[R ] 2 1 − Q Tr[R ]
∞

∞

∞

∞

Proof. From the ergodic Birkhoff’s theorem, L̂T converges almost surely
R to L̂∞ , and thus
−1 . We define the martingale M̂ = T L−1 (dM )X . We
L̂−1
converges
almost
surely
to
L̂
s
s
T
∞
T
Xs
0
have
1
−1
>
d(M̂s )i,j = Tr[e>
Tr[L−1
i,j LXs (dMs )Xs ] =
Xs (ei,j Xs + Xs ei,j )dMs ]
2
1 X
−1
=
(LX
(ei,j Xs + Xs e>
i,j ))k,l (dMs )k,l .
s
2
1≤k,l≤d

We get from (27) that
i
1 h −1
>
>
0
0
hd(M̂s )i,j , d(M̂s ) i = Tr LXs (ei,j Xs + Xs ei,j )(ei ,j Xs + Xs ei0 ,j 0 ) ds,
2
>
hd(M̂s )i,j , dNs i = Tr[L−1
Xs (ei,j Xs + Xs ei,j )]ds = δi,j ds.
i0 ,j 0

2 STATISTICAL INFERENCE: THE ERGODIC CASE

12

NT
√T , √
From the central limit theorem for martingales, ( M̂
) converges in law under Pθ towards
T
T

a centered Gaussian vector (Ĝ, Ĥ) taking values in Md × R such that
ii
1 h h
>
>
0 ,j 0 X∞ + X∞ e 0 0 )
Eθ (Ĝi,j Ĝi0 ,j 0 ) = Eθ Tr L−1
(e
X
+
X
e
)(e
i,j
∞
∞
i
i,j
i ,j
X∞
2
−1
2
Eθ (Ĝi,j Ĥ) = δi,j and Eθ (Ĥ ) = Q∞ .

(38)

We thus have the following Laplace transform for c ∈ Md , λ ∈ R,


1
1 2 −1
−1
>
>
>
Eθ [exp(Tr[c Ĝ] + λĤ)] = exp
Eθ [Tr[LX∞ (cX∞ + X∞ c )(cX∞ + X∞ c )]] + λ Tr[c] + λ Q∞ .
4
2
(39)
√

By using Slutsky’s Theorem, we get from (35) that
T (b̂T − b, α̂T − α) converges in law
under Pθ when T → +∞ to the centered Gaussian vector



−1
(G, H) = L̂−1
.
∞ (Ĝ − Q∞ ĤId ), 2Q∞ Ĥ − Tr[L̂∞ (Ĝ − Q∞ ĤId )]
Now, we use that L̂−1
∞ (Id ) =

−1
1
−1 R∞
1−Q∞ Tr[R∞ ]

and that L̂∞ is self-adjoint to get for c ∈ Md ,

−1

>

Tr[c G] =

>
Tr[L̂−1
∞ (c) Ĝ]

−

Q∞ Ĥ Tr[c> R∞ ]
1−

−1
Q∞ Tr[R∞ ]

2Q∞

, H=
1−

−1 (Ĥ
Q∞ Tr[R∞ ]

−1

− Tr[R∞ Ĝ])].

From (39), we obtain after some calculations (37), using in particular that for m ∈ Md , s ∈ Sd ,
>
>
Eθ [Tr[L−1
X∞ ((m + s)X∞ + X∞ (m + s) )((m + s)X∞ + X∞ (m + s) )]]

(40)

>
>
>
2
= Eθ [Tr[L−1
X∞ (mX∞ + X∞ m )(mX∞ + X∞ m )]] + 2 Tr[s(mR∞ + R∞ m )] + 2 Tr[s R∞ ]

and taking m = L̂−1
∞ (c) and s = −λ

2Q∞

−1
1−Q∞ Tr[R∞ ]

−1

R∞ .

Remark 2.1.
√ It is interesting to compare Theorems 2.1 and 2.3 and see that the asymptotic
variance of T (α̂T −α) is the same in both cases. Instead, for the estimation of the symmetric
part of b, we can check that the asymptotic variance is greater when we do not know a priori
−1
−1
>
that b is symmetric. For c ∈ Sd , we have c = 21 [L̂−1
∞ (c)R∞ + R∞ L̂∞ (c) ] − Q∞ Tr[L̂∞ (c)]Id
−1

−1
−1
−1
and c = LR
(c)R∞ + R∞ LR
(c) − 2Q∞ Tr[LR
(c)]Id . Multiplying by R∞ , we
∞ ,Q∞
∞ ,Q∞
∞ ,Q∞
get
−1
Tr[cR∞ ]
−1
−1
Tr[L̂∞ (c)] =
−1 = 2 Tr[LR∞ ,Q∞ (c)],
1 − Q∞ Tr[R∞ ]
−1
and then L̂−1
∞ (c) = 2LR

∞ ,Q∞

(c) + ∆ with ∆R∞ + R∞ ∆> = 0. This gives from (40)

1
−1
−1
>
−1
−1
>
Eθ [Tr[L−1
X∞ (L̂∞ (c)X∞ + X∞ L̂∞ (c) )(L̂∞ (c)X∞ + X∞ L̂∞ (c) )]]
4
!
−1
−1
1
Tr[cR∞ ]Q∞
Tr[cR∞ ]
−1
+
−1
−1 − Tr[L̂∞ (c)]
1 − Q∞ Tr[R∞ ] 2 1 − Q∞ Tr[R∞ ]
1
−1
−1
>
>
2
=2 Tr[LR
(c)2 R∞ ] + Eθ [Tr[L−1
X∞ (∆X∞ + X∞ ∆ )(∆X∞ + X∞ ∆ )]] − 2Q∞ Tr[LR∞ ,Q∞ (c)]
∞ ,Q∞
4
1
−1
−1
>
>
= Tr[LR
(c)c] + Eθ [Tr[L−1
X∞ (∆X∞ + X∞ ∆ )(∆X∞ + X∞ ∆ )]] ≥ Tr[LR∞ ,Q∞ (c)c],
,Q
∞
∞
4
since L−1
X∞ is a self-adjoint positive operator.

3 STATISTICAL INFERENCE: SOME NONERGODIC CASES

2.3

13

The MLE estimator of b

When α ∈ [d − 1, d + 1), we are no longer able to study the convergence of the MLE of α. It
is however still possible to get the speed of convergence of the MLE of b.
Theorem 2.4. Assume that −b ∈ Sd+,∗ and α > d − 1. For T > 0, we consider b̂T defined
√
by (18). Then, under Pθ , T (b̂T − b) converges in law to a centered
Gaussian

 vector G on
−1
Sd with the following Laplace transform Eθ [exp (Tr[cG])] = exp Tr[cLR (c)] , c ∈ Sd .
∞

Assume that −(b + b> ) ∈ Sd+,∗ and α > d − 1. For T > 0, we consider b̂T defined by (16).
√
Then, under Pθ , T (b̂T − b) converges in law to a centered Gaussian
vector G on Md with

the following Laplace transform: for c ∈ Md , Eθ exp Tr[c> G]


1
−1
−1
−1
>
−1
−1
>
= exp
Eθ [Tr[LX∞ (Ľ∞ (c)X∞ + X∞ Ľ∞ (c) )(Ľ∞ (c)X∞ + X∞ Ľ∞ (c) )]] ,
4
with Ľ∞ (c) = Eθ [L̄X∞ (Y )].
Proof. We could prove the result for (18) by using the explicit Laplace transform Proposition 5.1. Here, we use the same arguments as before based on the ergodic property. From (18),
we have




√
1
MT
−1
−1
.
T (b̂T − b) = L RT √ (XT − x − bRT − RT b − αT Id ) = L RT √
T
T
T
T
M
√T
T

converges in law to the centered Gaussian vector G̃
√
defined by(28).
Slutsky’s
theorem
and
(19)
give
then
the
convergence
of
T (b̂T − b) to

−1
G = LR
G̃ , whose Laplace transform is given by Lemma B.4.
∞
√
√
R
1 T
Similarly, we get from (16) that T (b̂T −b) = Ľ−1
T (M̂T / T ) with ĽT (Y ) = T 0 L̄Xs (Y )ds.
We easily check that Tr[ĽT (Y )Y ] ≥ Tr[L̂T (Y )Y ] and Tr[Ľ∞ (Y )Y ] ≥ Tr[L̂∞ (Y )Y ] for
Y ∈ Md . Therefore, ĽT and Ľ∞ are positive definite and thus invertible. Besides, the
ergodic theorem gives that Ľ−1
T converges almost surely to Ľ∞ . The result follows from (39)
and the self-adjoint property of Ľ∞ .
As in the proof of Theorem 2.1,

3

Statistical Inference of the Wishart process: some nonergodic cases

This section studies the convergence of the MLE in the case b = b0 Id with b0 > 0. When
b0 = 0 and α > d + 1, we are able to describe the rate of convergence of the MLE of (b, α)
given by (15), when b is known to be symmetric. When b0 > 0 and α > d − 1, we can
also obtain the rate of convergence of the MLE of b given by (18). Last, when b is known
a priori to be diagonal, the MLE of b has a simpler form and we can describe precisely its
convergence.

3.1

The global MLE of θ = (b, α) when b = 0

The following result provides the asymptotic behavior of the estimator of the couple when
α > d + 1 and b = 0 in (4).

3 STATISTICAL INFERENCE: SOME NONERGODIC CASES

14

Theorem 3.1. Assume
that b = 0 and α > d + 1. Let (b̂T , α̂T ) be the MLE defined by (15).
p
Then, (T (b̂T − b), log(T )(α̂T − α)) converges in law under Pθ when T → +∞ to
!
r

α
−
(d
+
1)
L−1
X10 − αId , 2
G ,
R10
d
p
Rtp
where Xt0 = αtId + 0 Xs0 dWs + dWs> Xs0 is a Wishart process with the same parameters
Rt
but starting from 0, Rt0 = 0 Xs0 ds and G ∼ N (0, 1) is an independent standard Normal
variable.
Proof. From (15) and (26), we obtain
 p
T Tr[b̂T ]

log(T )QT + 2 log(T )QT √ NT
log(T )(α̂T − α) = −2 √


log(T )
log(T )






T b̂T

= L−1
RT

T2


,QT

XT
T

−

x
T


− (QT ZT + 1 + d)Id ,


and we are interested in studying the convergence in law of

√ NT , XTT
log(T )

,

RT
T2


. By Theo-

rem 4.1 in [29], for µ > 0 and T > 1,
dP
= exp
dPθ

µ2
µN
p T −
log(T ) 2QT log(T )

!

defines a change of probability and (Xt )t∈[0,T ] is a Wishart process with degree α + √

µ
log(T )

under P. Let λ1 , λ2 ∈ Sd+,∗ and
"
!




#
µNT
µ2
λ2
λ1
AT = Eθ exp p
−
exp − Tr
XT − Tr
RT
.
T
T2
log(T ) 2QT log(T )
By Proposition 5.1, we have






λ2
λ1
AT = E exp − Tr
XT − Tr
RT
T
T2

 0 −1 
1
1
1
√
=
exp
−
Tr
V V x
→
α ,
α+µ/ log(T )
T →+∞ det[V ] 2
2T
2
det[V ]

(41)

where
p
p
p
p
p
p
V = ( 2λ1 )−1 sinh( 2λ1 )2λ2 + cosh( 2λ1 ), V 0 = 2 cosh( 2λ1 )λ2 + 2λ1 sinh( 2λ1 ).
(42)
0
0
We note that this limit does not depend on µ and is the Laplace transform of (X1 , R1 ) by
Proposition 5.1.
1
d
We now use that QT log(T
→ α−(d+1)
a.s., see Lemma C.2 and we define
)
T →+∞





λ2
λ1
µNT
µ2 d
ÃT = Eθ [exp(ξT )] and ξT = p
−
− Tr
XT − Tr
RT
T
T2
log(T ) 2(α − (d + 1))

3 STATISTICAL INFERENCE: SOME NONERGODIC CASES

15

that is finite by using equation (73) of Lemma C.2 since ξT 6 √µNT . We have
log(T )


AT = ÃT + Eθ


exp

µ2 d
µ2
−
2(α − (d + 1)) 2QT log(T )

The Cauchy-Schwarz inequality gives
"

2

(AT − ÃT ) 6 Eθ

exp

µ2 d
µ2
−
2(α − (d + 1)) 2QT log(T )





− 1 exp (ξT ) .

2 #


−1

Eθ [exp (2ξT )] .

Since QT log(T ) is positive for T > 1 and converges a.s. to α−(d+1)
, the first expectation
d
goes to 0 while the second one is bounded by using again (73). Therefore, AT − ÃT → 0,
T →+∞

and we get
"



Eθ exp


√ NT





λ
λ
µN
p T − Tr 2 XT − Tr 12 RT
T
T
log(T )

Thus,
, XT
log(T ) T
independent of X 0 .

,

RT
T2


converges in law to (

q

!#

exp
→

T →+∞

d
0
0
α−(d+1) G, X1 , R1 ),



µ2 d
2(α−(d+1))
α

det[V ] 2


.

where G ∼ N (0, 1) is

From (25), we have
NT
1
log(T )QT p
+ α,
QT Z T + 1 + d = 2 p
log(T )
log(T )

and therefore QT ZT + 1 + d converges in probability to α. Slutsky’s theorem gives then the
following convergence in law: as T → +∞,
s
!
!
NT
XT RT
d
α
−
(d
+
1)
p
,
,
, QT ZT + 1 + d, QT log(T ) ⇒
G, X10 , R10 , α,
.
α − (d + 1)
d
log(T ) T T 2
(43)
This gives the claimed convergence for (α̂T , b̂T ) due to the continuity property given in
Lemma B.1.
Theorem 3.2. Assume that b = 0 and α = d + 1. Let (b̂T , α̂T ) be the MLE defined by (15).
Then, (T (b̂T − b), log(T )(α̂T − α)) converges in law under Pθ when T → +∞ to


 4
−1
0
LR0 X1 − αId ,
,
1
dτ1
p
Rtp
where Xt0 = αtId + 0 Xs0 dWs + dWs> Xs0 is a Wishart process with the same parameters
Rt
but starting from 0, Rt0 = 0 Xs0 ds and τ1 = inf{t > 0, Bt = 1} where B is a standard
Brownian motion independent from W .
Proof. The proof follows the same line as the one of Theorem 3.1, but we now write
log(T )(α̂T − α) = −2

T Tr[b̂T ]
NT
log(T )2 QT + 2 log(T )2 QT
,
log(T )
log(T )

3 STATISTICAL INFERENCE: SOME NONERGODIC CASES
while we still have T b̂T = L−1
RT

T2

dP
dP
Sd+,∗ ,

for µ > 0 and T > 1,

= exp


,Q

T


− (QT ZT + 1 + d)Id . By Theorem 4.1 in [29],

µ2
− 2QT log(T
defines a change of probability, and we
2
)

XT
T

µNT
log(T )

16

−

x
T

define for λ1 , λ2 ∈








λ2
λ1
µNT
µ2
.
exp − Tr
RT
AT = Eθ exp
−
XT − Tr
log(T ) 2QT log(T )2
T
T2
By Proposition 5.1, we have
1

AT =
det[V ]

α+µ/log(T )
2




1
1
exp −
Tr V 0 V −1 x
→
α ,
T →+∞ det[V ] 2
2T

where V and V 0 are defined by (42).
NT
d
We now use that log(T
) → 2 in probability, see Lemma C.2, and define




λ1
µd
µ2
λ2
ÃT = Eθ [exp(ξT )] , ξT =
−
− Tr
XT − Tr
RT
2
2QT log(T )2
T
T2
and have






µNT
µd
exp
−
− 1 exp (ξT ) .
log(T )
2
 
We note that exp(ξT ) 6 exp µd
2 . By using Lemma C.2 and the uniform integrability (74),
AT = ÃT + Eθ

we get that AT − ÃT

→

T →+∞

0 and therefore





λ2
λ1
exp(−µd/2)
µ2
.
− Tr
XT − Tr
RT
→
α
2
2
T →+∞
2QT log(T )
T
T
det[V ] 2



2 
Therefore, XTT , RT T2 , QT log(T )2 converges in law to X10 , R10 , d2 τ11 , where τ1 is indepen

Eθ exp −

dent of X 0 . We observe that QT ZT =
theorem gives


1
log(T ) QT

NT
XT RT
,
,
, QT ZT + 1 + d, QT log(T )2
log(T ) T T 2

ZT
log(T )2 log(T
) . Lemma C.2 and Slutsky’s


⇒

d 0 0
, X , R , d + 1,
2 1 1

 2 !
2
1
,
d
τ1

(44)

which gives the claim by using the formulas for T b̂T and log(T )(α̂T − α).

3.2

The MLE of b

Until the end of this section we consider that α > d − 1 is known and study the speed of
convergence of the estimator of b defined by (18).
3.2.1

Case b = 0.

Theorem 3.3. Assume that b = 0 and α > d − 1. For T > 0, let b̂T be defined by (18).
When T → +∞, T (b̂T − b) converges in law under Pθ to L−1
X10 − αId , where (Xt0 )t>0 is
R10
p
Rtp
Rt
the solution to Xt0 = αtId + 0 Xs0 dWs + dWs> Xs0 and Rt0 = 0 Xs0 ds.

3 STATISTICAL INFERENCE: SOME NONERGODIC CASES
Proof. From (18), we have T b̂T = L−1
RT

T2



XT
T

−

x
T

17


− αId . Let V and V 0 be defined by (42).

Similarly to (41), we have by Proposition 5.1 for λ1 , λ2 ∈ Sd+,∗






 1
 0 −1 
λ1
1
1
λ2
=
XT − Tr
R
Tr
V V x
→
exp
−
Eθ exp − Tr
α
α .
T
T
→+∞
T
T2
2T
2
det[V ]
det[V ] 2


This gives the convergence in law of XTT , RT T2 to (X10 , R10 ) and then the claimed result.
3.2.2

Case b = b0 Id , b0 > 0.

In this case b = b0 Id with b0 > 0. In order to identify the speed of convergence and the limit
law, we use the Laplace transform approach. We have the following result,
Theorem 3.4. Assume that b = b0 Id , b0 > 0, and α > d − 1. For T > 0let b̂T defined by
√ 
√
(18). When T → +∞, exp(b0 T )(b̂T − b) converges in law under Pθ to L−1
X G̃ + G̃ X
X


where X ∼ W ISd 2bx0 , α, 0, Id ; 4b12 and G̃ is an independent d-square matrix whose elements
0
are independent standard Normal variables.
The proof of this results relies on the explicit calculation of the Laplace transform of
(XT , RT ) and is postponed to Subsection 5.2.
Obviously, the case b = b0 Id is very particular. One would like to consider more general
nonergodic cases or ideally to be able to state a general convergence results of b̂T towards b
for any b ∈ Sd . Despite our efforts, we have not been able to get such a result. The reason
why we can handle the ergodic case and the nonergodic case b = b0 Id with√b0 > 0 is that the
convergence of all the matrix terms occurs at the same speed, namely 1/ T for the ergodic
case, 1/T for b = 0 and e−b0 T when b0 > 0. In the other cases, there is no such a simple
scalar rescaling. Heuristically, there may be different speeds of convergence that are difficult
to disentangle because of the different matrix products. To get an idea of this, we present
now the case of the estimation of b when b is known to be a diagonal matrix. In this case,
we obtain different speed of convergence for each diagonal terms.
3.2.3

The MLE of b when b is known a priori to be diagonal.

We assume that α > d − 1 is known and that b is a diagonal matrix, i.e. b = diag(b1 , · · · , bd ).
We want to estimate the diagonal elements by maximizing the likelihood. We denote θ0 = (0, α).
As in (17), we have
 Tr[bX ] − Tr[bx] 1 Z T

αT
T
θ,θ0
LT = exp
−
Tr[b] .
Tr[b2 Xs ]ds −
2
2 0
2
By differentiating this with respect to bi , 1 6 i 6 d, we get
0
∂bi Lθ,θ
T
0
Lθ,θ
T

1
(XT )i,i − xi,i − αtId − 2bi
=
2

Z

T


(Xs )i,i ds ,

0

and therefore the MLE of b is given by
(b̂T )i =

(XT )i,i − xi,i − αT
.
2(RT )i,i

(45)

4 OPTIMALITY OF THE MLE

18

We therefore obtain
(b̂T )i − bi =

(XT )i,i − xi,i − αT − 2bi (RT )i,i
.
2(RT )i,i

(46)

Let us observe that this estimator is precisely the one obtained by Ben Alaya and Kebaier [5]
for the CIR process. This is not very surprising since we know from (4), (2) and b diagonal
that there exists independent Brownian motions β i , 1 6 i 6 d such that
q
d(Xt )i,i = (α + 2bi (Xt )i,i )dt + 2 (Xt )i,i dβti .
Thus, the diagonal elements follow independent CIR processes, and the observation of the
non diagonal elements does not improve the ML estimation. We can obtain the asymptotic convergence by applying Theorem 1 in [4], up to a small correction in the nonergodic
case which is given by our Theorem 3.4 in dimension d = 1. This yields to the following
proposition.
Proposition 3.1. Let α > d − 1 and b a diagonal matrix. Let t = diag(1t , · · · , dt ) be a
diagonal matrix with

1
t− 2
if bi < 0

it =
t−1
if bi = 0

exp(−bi t)
if bi > 0
Then, under Pθ , −1
T diag((b̂T )1 − b1 , · · · , (b̂T )d − bd ) converges in law to a diagonal matrix D
made with independent elements. Each diagonal element Di is distributed as follows:

q
−2bi


if bi < 0

α G


X10 −α
if bi = 0
∀i ∈ {1, · · · , d}, Di =
2R10
law 
G

if bi > 0
 r xi,i /(2bi ) ,


X
2
1/(4b )
i

where Xtx = x + αt + 2

4

Rtp
0

Xsx dWs , Rt0 =

Rt
0

Xs0 ds, and G ∼ N (0, 1) is independent of X.

Optimality of the MLE

In parametric estimation theory, a fundamental role is played by the local asymptotic normality (LAN) property since the work of Le Cam [24]. This general concept developed by
Le Cam is extended later by Le Cam and Yang [25] and Jeganathan [22] to local asymptotic mixed normality (LAMN) and local asymptotic quadraticity (LAQ) properties. These
notions are mainly dedicated to study the asymptotic efficiency of estimators of a given parametric model. The aim of this section is to check the validity of either LAN, LAMN or
LAQ properties for the global model in order to get the asymptotic efficiency of our maximum likelihood estimators studied in the previous section. Here we prove these properties
only for the global model θ = (b, α) when b is known to be symmetric. The same technique
applies for all the other cases considered in this paper where we have been able to obtain the
corresponding local asymptotic property.

4 OPTIMALITY OF THE MLE

19

Let us consider the Wishart process (Xt )t>0 ∈ Sd+ with parameters θ := (b, α), with
α ≥ d + 1 and b ∈ Sd .
√
√

dXt = [αI + bXt + Xt b] dt + X t dWt + dWt> X t , t > 0
(47)
X0 ∈ Sd+ .
We recall that Pθ denotes the distributions induced by the solutions of (47) on canonical
space C(R+ , Sd+ ) with the natural filtration FtX := σ(Xs , s ≤ t) and Pθ,t = Pθ |FtX denotes
the restriction of Pθ on the filtration Ft .
For α̃ ≥ d + 1 and b̃ ∈ Sd+ , we set θ̃ = (b̃, α̃),
Ht =

p
α̃ − α p
( Xt )−1 + (b̃ − b) Xt ,
2

and we introduce the log-likelihood function
 Z T

Z
dPθ̃,T
1 T
θ
Tr[Hs dWs ] −
=
Tr[Hs Hs> ]ds
`T (θ̃) = log
dPθ,T
2
0
0

(48)

Rt
The process (W̃t = Wt − 0 Hs> ds, t ≤ T ) is a d × d-Brownian motion under Pθ̃,T . In
the sequel, let us introduce the quantity δT := (δ1,T , δ2,T ) ∈ R2 where for i ∈ {1, 2} the
localizing rates satisfy |δi,T | → 0 when T → ∞. For all u := (u1 , u2 ) ∈ R × Sd , we define
δT · u := (δ1,T u1 , δ2,T u2 ) ∈ R × Sd . Now, we rewrite (48) with θ̃ = θ + δT · u
`θT (θ

√

p
δ1,T u1 ( Xs )−1
Tr
+ δT · u) =
dWs + δ2,T u2 Xs dWs
2
0
Z
1 T
T
−
Tr [δ2,T u2 Xs δ2,T u2 ] ds − δ1,T u1 Tr [δ2,T u2 ]
2 0
2
Z T


1
(δ1,T u1 )2 Tr (Xs )−1 ds.
−
8 0
Z

T



Hence, by using the definitions (13), (23) and (24) of the martingales processes (Nt )t≥0 and
(Mt )t≥0 and the processes (Rt )t≥0 and (Qt )t≥0 , it is easy to check that

 1
1
δ1,T u1 NT + Tr δ2,T u2 MT − Tr [δ2,T u2 RT δ2,T u2 ]
2
2
1
T
2
− δ1,T u1 Tr [δ2,T u2 ] − (δ1,T u1 ) (QT )−1
2
8
1
= ΛT (u) − ΓT (u),
(49)
2



where ΛT (u) = 21 δ1,T u1 NT + Tr δ2,T u2 MT
is a linear random function with respect to
`θT (θ + δT · u) =

u ∈ R × Sd+ with quadratic variation
1 2 2 −1
2
ΓT (u) = δ2,T
Tr[u22 RT ] + T δ1,T δ2,T u1 Tr [u2 ] + δ1,T
u1 QT .
4

4 OPTIMALITY OF THE MLE

20

Case −b ∈ Sd+ and α > d + 1

4.1

We first consider α > d + 1. In this ergodic case, we set δi,T = T −1/2 for i ∈ {1, 2}, and we
get from (19) and (21)




1
a.s.
ΓT (u) −→ Γ∞ (u) := Tr u22 R∞ + u1 Tr [u2 ] + u21 Tr (Q∞ )−1 ,
4

as T → +∞.

(50)

This yields the validity of the so called Raykov type condition. Hence, according to Theorem 1
in [28], relations (49) and (50) ensure the validity of the local asymptotic normality (LAN)
property, that is under Pθ we have


1/2
(ΛT (u), ΓT (u)) ⇒ Γ∞ (u)Z, Γ∞ (u) , as T → ∞,
(51)
with Z a standard normal real random variable. It is worth noting that the above convergence
can also be obtained using the proof of Theorem 2.1. In fact, we have already proven that
under Pθ


NT MT
√ ,√
⇒ (G̃, H̃)
(52)
T
T
where (G̃, H̃) is a centered Gaussian vector taking values in R × Sd+ such that


Eθ (G̃i,j G̃k,l ) = δjl (R∞ )i,k + δjk (R∞ )i,l + δil (R∞ )j,k + δik (R∞ )j,l ,
−1

Eθ (G̃i,j H̃) = 2δi,j and Eθ (H̃ 2 ) = Q∞ .
Therefore, LAN property (51) follows from relations (50) and (52).
We now consider the case α = d + 1 and set δ1,T = T −1 and δ2,T = T −1/2 . By using (33),
we get that under Pθ ,

 q


u1
1 2
2
(ΛT (u), ΓT (u)) ⇒ − Tr[b] + Tr u2 R∞ G̃ , Tr[u2 R∞ ] + u1 τ− Tr[b] , as T → ∞,
2
4
where τ− Tr[b] is defined as in Theorem 2.2 and G̃ is an independent matrix, whose elements
G̃i,j , 1 6 i, j 6 d, are independent standard normal variables. Hence, according to Le Cam
and Yang [25] and Jeganathan [22] this last convergence yields the LAQ property for this
ergodic case.

Case b = 0 and α > d + 1

4.2

We first assume α > d + 1. From (47) with b = 0 and (23), we have MT = XT − x − αId T .
From (43), it follows that as T → ∞
s
!
!
MT RT Q−1
d
d
NT
T
p
,
,
,
⇒
G, X10 − αId , R10 ,
,
α − (d + 1)
α − (d + 1)
log(T ) T T 2 log(T )
where X10 and R10 are defined as in Theorem 3.1. Thus, in the same way as in the previous
case if we set δ1,T = √ 1
and δ2,T = T −1 , then (ΛT (u), ΓT (u)) converges in law under Pθ
log(T )

to
1
2

s

h
i

d
1 
u21 d
0
u1 G + Tr u2 (X10 − αId ) , Tr u22 R1 +
α − (d + 1)
2
4(α − (d + 1))

!
, as T → ∞.

5 THE LAPLACE TRANSFORM AND ITS USE TO STUDY THE MLE

21

This ensures the validity of the LAQ property in this non-ergodic case.
When α = d + 1, we use the notation of Theorem 3.2 and get from (44)
!
 2 !
Q−1
NT
d
MT RT
d
T
⇒
τ1 .
,
, 2,
, X10 − αId , R10 ,
2
log(T ) T T log(T )
2
2
With δ1,T =

1
log(T )



and δ2,T = T −1 , we get that (ΛT (u), ΓT (u)) converges in law under Pθ to

h
i d2 u2 
d
1
0
2 0
1
u1 + Tr[u2 (X1 − αId )], Tr u2 R1 +
τ1 , as T → ∞.
4
2
8

This gives again the LAQ property.

5

The Laplace transform and its use to study the MLE

5.1

The Laplace transform of (XT , RT )

We present our main result on the joint Laplace transform of (XT , RT ), that can be of
independent interest. This Laplace transform is given by Bru [7], eq. (4.7) when b = 0 and
has been recently studied and obtained explicitly by Gnoatto and Grasselli [17]. Here, we
present another proof that enables us to get the Laplace transform for any α > d − 1, as well
as a more precise result concerning its set of convergence, see Remarks 5.1 and 5.2 below for
a further discussion.
Proposition 5.1. Let α > d − 1, x ∈ Sd+ , b ∈ Sd and X ∼ W ISd (x, α, b, Id ). Let v, w ∈ Sd
be such that
v
w
∃m ∈ Sd , − mb − bm − 2m2 ∈ Sd+ and
(53)
+ m ∈ Sd+ .
2
2
Then, we have for t > 0
h
 1
i
1
E exp − Tr[wXt ] − Tr[vRt
2
 2

exp − α2 Tr[b]t
1  0
=
exp − Tr (Vv,w
(t)Vv,w (t)−1 + b)x ,
α
2
det[Vv,w (t)] 2

(54)

with
Vv,w (t) =

∞
X
k=0

ṽ k
t2k+1
(2k + 1)!

!
w̃ +

∞
X
k=0

t2k

ṽ k
,
(2k)!

ṽ = v + b2 ,

and

w̃ = w − b.

√ −1
√
√
If besides ṽ = v√+ b2 ∈ Sd+,∗ ,√we have
V
(t)
=
(
ṽ)
sinh(
ṽt)
w̃
+
cosh(
ṽt) and then
v,w
√
0 (t) = cosh( ṽt)w̃ + sinh( ṽt) ṽ.
Vv,w
Before proving this result, we recall the following fact:
∀x, y ∈ Sd+ , Tr[xy] > 0,

(55)

√ √
√ √
which is clear once we have observed that Tr[xy] = Tr[ xy x] and xy x ∈ Sd+ . We also
recall a result on matrix Riccati equations, see Dieci and Eirola [11] Proposition 1.1.

5 THE LAPLACE TRANSFORM AND ITS USE TO STUDY THE MLE

22

Lemma 5.1. Let b̃ ∈ Sd and δ̃ ∈ Sd+ . Let ξ denote the solution of the following matrix
Riccati differential equation
ξ 0 + 2ξ 2 = b̃ξ + ξ b̃ + δ̃, ξ(0) ∈ Sd .

(56)

If ξ(0) ∈ Sd+ , the solution ξ is well-defined for any t > 0 and satisfies ξ(t) ∈ Sd+ .
Proof
5.1. Let
T > 0 be given. We first assume w, v ∈ Sd+ , which ensures that
 of Proposition

E exp − 21 Tr[wXT + vRT ] < ∞. We consider the martingale




1
Mt = E exp − Tr[wXT + vRT ] Ft , t ∈ [0, T ].
2
Due to the affine structure, we are looking for smooth functions β : R+ → R, γ, δ : R+ → Sd
such that
Mt = exp (β(T − t) + Tr[γ(T − t)Xt ] + Tr[δ(T − t)Rt ]) .
We necessarily have β(0) = 0, γ(0) = −w/2 and δ(0) = −v/2. Itô’s formula gives
dMt n
= − β 0 (T − t) − Tr[γ 0 (T − t)Xt ] − Tr[δ 0 (T − t)Rt ] + Tr[γ(T − t)(αId + bXt + Xt b)]
Mt
o
p
p
+ Tr[δ(T − t)Xt ] + 2 Tr[γ(T − t)2 Xt ] dt + Tr[γ(T − t)( Xt dWt + dWt> Xt )].

Since M is a martingale, the drift term should vanish almost surely. The drift term being
a (deterministic) affine function of (Xt , Rt ), we obtain the following system of differential
equations:
δ 0 = 0,
0

(57)
2

− γ + γb + bγ + 2γ + δ = 0,
0

− β + α Tr[γ] = 0.

(58)
(59)

The first equation gives δ(t) = −v/2. The second equation is a matrix Riccati differential
equation. We now consider ξ = m − γ with m satisfying (53). It solves (56) with b̃ = b + 2m,
δ̃ = −bm − mb − 2m2 + v/2 and ξ(0) = m + w/2. We know then by Lemma 5.1 that ξ is
well defined for any t > 0 and stays in Sd+ . In particular, γ is well defined for any t > 0.
We set γ̃ = γ + 12 b. We have γ 2 = γ̃ 2 − 12 (bγ̃ + γ̃b) + 14 b2 and thus γ̃ solves the following
matrix Riccati differential equation:
1
1
γ̃ 0 = 2γ̃ 2 − ṽ, γ̃(0) = − w̃, with ṽ = v + b2 and w̃ = w − b.
2
2


 

M1 (t) M2 (t)
0
−ṽ/2
We set M (t) =
= exp t
∈ M2d and get by Levin [26]
M3 (t) M4 (t)
−2Id
0
that


−1
1
1
γ̃(t) = M2 (t) − M1 (t)w̃ M4 (t) − M3 (t)w̃
.
2
2
We check that the matrix M4 (t) − 21 M3 (t)w̃ is indeed invertible. In fact, let




1
τ = inf t > 0, det M4 (t) − M3 (t)w̃ = 0 .
2

5 THE LAPLACE TRANSFORM AND ITS USE TO STUDY THE MLE

23



d
We have τ > 0 and for t ∈ [0, τ ), dt
[M4 (t) − 21 M3 (t)w̃] = −2 M2 (t) − 21 M1 (t)w̃ and thus




d
1
1
det M4 (t) − M3 (t)w̃ = −2 det M4 (t) − M3 (t)w̃ Tr[γ̃(t)].
dt
2
2


Rt
This gives det M4 (t) − 21 M3 (t)w̃ = exp(−2 0 Tr[γ̃(s)]ds) > 0, and we necessary get τ = +∞
since γ and thus γ̃ is well defined for t > 0.
Since

2k  k


2k+1 

0
−ṽ/2
ṽ
0
0
−ṽ/2
0
−ṽ k+1 /2
=
and
=
,
−2Id
0
−2Id
0
0 ṽ k
−2ṽ k
0
we get
M1 (t) = M4 (t) =

∞ 2k k
X
t ṽ
k=0

(2k)!

, M2 (t) = −

∞

∞

k=0

k=0

X t2k+1 ṽ k
1 X t2k+1 ṽ k+1
, M3 (t) = −2
.
2
(2k + 1)!
(2k + 1)!

√
√
and we have
M1 (t) = M4 (t) = cosh(t ṽ),
, ṽ is well defined
If ṽ = v +√b2 ∈ Sd+,∗
√
√
√
M2 (t) = − 12 ṽ sinh(t ṽ) and M3 (t) = −2( ṽ)−1 sinh(t ṽ). Now, we define
!
∞
∞
X
X
ṽ k
1
ṽ k
2k+1
t
V (t) = M4 (t) − M3 (t)w̃ =
t2k
.
w̃ +
2
(2k + 1)!
(2k)!
k=0

k=0

Since V 0 (t) = −2M2 (t) + M1 (t)w̃, we obtain that

1
1
γ̃(t) = − V 0 (t)V (t)−1 and thus γ(t) = − V 0 (t)V (t)−1 + b .
2
2
Last, we have β 0 (t) = − 21 α Tr[V 0 (t)V (t)−1 ] − 12 α Tr[b] and we obtain that
1
1
β(t) = − α log(det[V (t)]) − α Tr[b]t,
2
2
(t)]
= det[V (t)] Tr[V 0 (t)V (t)−1 ].
since d det[V
dt
It remains to show that we indeed have (54) for v and w satisfying (53). We define
exp(β(T −t)+Tr[γ(T −t)Xt ]+Tr[− v2 Rt ])
Et =
. By Itô’s formula, we have
exp(β(T )+Tr[γ(T )x])

p
p
dEt
= Tr[γ(T − t)( Xt dWt + dWt> Xt )].
Et
This is a positive local martingale and thus a supermartingale which gives E[ET ] 6 1, and we
want to prove that this is a martingale. To do so, we use the argument presented by Rydberg
in [34]. For L > 0, we define
τL = inf{t > 0, Tr[Xt ] > L},
L
and πL (x) = x1Tr[x]6L + Tr[x]
x1Tr[x]>L for x ∈ Sd+ . We consider (EtL , t ∈ [0, T ]) the solution
of
p
p
dEtL = EtL Tr[γ(T − t)( πL (Xt )dWt + dWt> πL (Xt ))], E0L = 1.

5 THE LAPLACE TRANSFORM AND ITS USE TO STUDY THE MLE
We clearly have E[ETL ] = 1. Besides, under PL given by

dPL
dP FT

24

= ETL , the process

p
dWtL = dWt − 2 πL (Xt )γ(T − t)dt, t ∈ [0, T ],
is a matrix Brownian motion. Since Et = EtL for t 6 τL , we have E[ET ] = E[ETL 1τL >T ]+E[ET 1τL 6T ].
By Lebesgue’s theorem, we get E[ET 1τL 6T ] → 0. On the other hand, E[ETL 1τL >T ] = PL (τL > T ).
L→+∞

Let us consider the Wishart process X̃ starting from x such that
p
p
h
i
dX̃t = αId + (b + 2γ(T − t))X̃t + X̃t (b + 2γ(T − t)) dt + X̃ t dWt + dWt> X̃ t .
We also define τ̃L = inf{t ∈ [0, T ], Tr[X̃t ] > L} with convention inf ∅ = +∞. The process X̃
solves the same SDE on [0, τ̃L ∧ T ] under P as X on [0, τL ∧ T ] under PL . We therefore have
PL (τL > T ) = P(τ˜L > T )

→

L→+∞

1,

which finally gives E[ET ] = 1.
Corollary 5.1. Let Y ∼ W ISd (y, α, b, a) be a Wishart process with parameters α > d − 1,
y ∈ Sd+ , a, b ∈ Md satisfying
ba> a = a> ab> and a invertible.

(60)

Let v, w ∈ Sd be such that
∃m ∈ Sd ,

1
ava>
awa> + m ∈ Sd+ and
− ab> a−1 m − m(a> )−1 ba> − 2m2 ∈ Sd+ .
2
2

(61)

Then, we have




Z T
1
E exp − Tr wYT + v
Ys ds
2
0

 1 

exp − α2 Tr[b]t
0
−1
> −1 >
> −1
−1
=
,
exp
−
Tr
(V
(t)V
(t)
+
(a
)
ba
)(a
)
ya
α
v,w
v,w
2
det[Vv,w (t)] 2

P
P∞ 2k ṽk
∞
2k+1 ṽ k
with Vv,w (t) =
k=0 t (2k)! and
k=0 t
(2k)! w̃ +
ṽ = ava> + (a> )−1 b2 a> ,

and

w̃ = awa> − (a> )−1 ba> .

Proof. We know that Y = a> Xa with x = (a> )−1 ya−1 and X ∼ W ISd (x, α, (a> )−1 ba> , Id ),
law

see e.g. equation (13) in [1]. We notice that (a> )−1 ba> = ab> a−1 ⇐⇒ ba> a = a> ab> and
thus (a> )−1 ba> ∈ Sd . We have








Z T
Z T
1
1
E exp − Tr wYT + v
Ys ds
= E exp − Tr awa> XT + ava>
Xs ds
,
2
2
0
0
which gives the result by applying Proposition 5.1.

5 THE LAPLACE TRANSFORM AND ITS USE TO STUDY THE MLE

25

By setting m̃ = a−1 m(a> )−1 , the condition (61) is equivalent to the existence of m̃ ∈ Sd ,
such that
1
v
w + m̃ ∈ Sd+ and − b> m̃ − m̃b − 2m̃a> am̃ ∈ Sd+ .
(62)
2
2
The case m = 0 gives back the finiteness of the Laplace transform when v, w ∈ Sd+ . If we
take m̃ = −w/2, we get also the finiteness when
v + b> w + wb − wa> aw ∈ Sd+ .

(63)

Another interesting choice is m = − 12 (a> )−1 ba> . We have m ∈ Sd from (60). This choice
gives the finiteness of the Laplace transform when v+b> (a> a)−1 b ∈ Sd+ and w−(a> a)−1 b ∈ Sd+ .
Let us note that ṽ = a(v + b> (a> a)−1 b)a> so that the first condition is the same as ṽ ∈ Sd+ .
Another interesting choice of m is given by the next remark.
Remark 5.1. Proposition 5.1 extends the result of Gnoatto and Grasselli [17] to α > d − 1,
and the sufficient condition (61) that ensures the finiteness of the Laplace transform is also
less restrictive, which is crucial in our study especially in the nonergodic case. In particular,
it does not assume a priori that v +b> (a> a)−1 b ∈ Sd+ . We can recover the result of [17] as folp
> −1 >
lows. Let us assume v+b> (a> a)−1 b ∈ Sd+ and take m = − (a ) 2 ba + 12 a(v + b> (a> a)−1 b)a> .
>
+
> −1
> −1 >
2
We have m ∈ Sd from (60) and it satisfies ava
2 − ab a m − m(a ) ba − 2m = 0 ∈ Sd .
Therefore, (61) holds if
q
w − (a> a)−1 b + a−1 a(v + b> (a> a)−1 b)a> (a> )−1 ∈ Sd+ .
This is precisely the condition stated in [17].


RT
Remark 5.2. It is possible to get similarly the Laplace transform of YT , 0 Ys ds when Y
solves
i
h
√
√
dYt = α + bYt + Yt b> dt + Y t dWt a + a> dWt> Y t , Y0 = y ∈ Sd+ ,
with a, b satisfying (60) and α − (d − 1)a> a ∈ Sd+ . Again, equation (13) in [1] gives
Y = a> Xa, where
law

h
i
√
√
dXt = α̂ + b̂Xt + Xt b̂> dt + X t dWt + dWt> X t , X0 = x,
with x = (a> )−1 ya−1 ∈ Sd , b̂ = (a> )−1 ba> ∈ Sd and α̂ = (a> )−1 αa−1 ∈ Sd . Repeating the
proof of Proposition 5.1, we observe that the Riccati equation (58) and equation (57) remain
unchanged while (59) is replaced by
1
1
β 0 = Tr[α̂γ] = − Tr[α̂V 0 (t)V (t)−1 ] − Tr[α̂b̂].
2
2
Therefore, we deduce that under the same condition (61), we have



Z T
1
E exp − Tr wYT + v
Ys ds
2
0
 1 

0
= exp(β(T )) exp − Tr (Vv,w
(t)Vv,w (t)−1 + (a> )−1 ba> )(a> )−1 ya−1 ,
2


5 THE LAPLACE TRANSFORM AND ITS USE TO STUDY THE MLE

26

Rt
t
0 (s)V
−1
> −1
with β(t) = − 12 0 Tr[(a> )−1 αa−1 Vv,w
v,w (s) ]ds − 2 Tr[α(a a) b] and Vv,w (t) defined
as in Corollary 5.1. Thus, the formula is no longer totally explicit. In Gnoatto and Grasselli [17], the result is stated with Tr[(a> )−1 αa−1 log(Vv,w (t))] instead of the first integral.
0 (s) and V
However, this replacement does not seem clear to us unless Vv,w
v,w (s) commute for
all s > 0 (this happens when the matrices ṽ and w̃ in Vv,w commute) or α = αa> a by using
the trace cyclic theorem.
Corollary 5.2. Let Y ∼ W ISd (y, α, b, a) be a Wishart process with parameters such that
ba> a = a> ab> and a invertible. Then,

Z T

Z
p
1 T
>
∀u ∈ Sd , E exp
Tr[u Ys dWs a]ds −
Tr[auYs ua ]ds
= 1.
2 0
0
h
i
RT
RT
√
Proof. We have 2 0 Tr[u Ys dWs a]ds = Tr[u(YT −y)]−αT Tr[ua> a]−Tr (ub + b> u) 0 Ys ds .
We apply Corollary 5.1 with w = −u and v = ub + b> u + ua> au. Therefore, (63) holds.
We then have w̃ = −(aua> + (a> )−1 ba> ) and ṽ = w̃2 and the result follows by simple
calculations.

5.2

Study of the MLE of b with the Laplace transform

We consider  : R+ → R∗+ a (deterministic) decreasing function such that limt→+∞ t = 0.
From the definition of the MLE of b (18), we get that
1
(b̂T − b) = L−1
( [XT − x − αT Id − bRT − RT b]).
2T RT T
T
Thus, we want to calculate the Laplace transform of (T [XT − x − αT Id − bRT − RT b], 2T RT )
in order to study the convergence of 1T (b̂T − b). For λ1 , λ2 ∈ Sd , we define
h

i
(64)
E(T, λ1 , λ2 ) := Eθ exp − T Tr[λ2 (XT − x − αT Id − bRT − RT b)] − 2T Tr[λ1 RT ]
h

i
= exp (T Tr[λ2 (x + αT Id )]) Eθ exp − Tr[T λ2 XT ] − Tr[(2T λ1 − T (λ2 b + bλ2 ))RT ] .
(65)
We now consider λ1 , λ2 ∈ Sd such that
λ1 − 2λ22 ∈ Sd+,∗ .

(66)

We define
vT = 2λ1 2T − 2(bλ2 + λ2 b)T ,

ṽT = vT + b2 ,

wT = 2λ2 T ,

w̃T = wT − b,

(67)

and have vT + bwT + wT b − wT2 = 2T (2λ1 − 4λ22 ) ∈ Sd+,∗ . Thus, by applying Proposition 5.1
with m = −T λ2 , we get that E(T, λ1 , λ2 ) is finite and given by

 1 

exp − α2 Tr[b]T
0
−1
E(T, λ1 , λ2 ) =
exp
−
Tr
(V
(T
)V
(T
)
+
b)x
α
vT ,wT
vT ,wT
2
det[VvT ,wT (T )] 2




× exp T Tr λ2 (x + αT Id )
(68)

5 THE LAPLACE TRANSFORM AND ITS USE TO STUDY THE MLE

27

with
p
p
p
VvT ,wT (T ) =( ṽT )−1 sinh( ṽT T )w̃T + cosh( v˜T T )
p
p
p
Vv0T ,wT (T ) = cosh( ṽT T )w̃T + sinh( ṽT T ) ṽT .
Besides, we have ṽT = (b − 2T λ2√
)2 + 2T (2λ1 − 4λ22 ) ∈ Sd+,∗ .
+,∗
When −b ∈ Sd and T = 1/ T , we can make explicit calculations and get
lim E(T, λ1 , λ2 ) = exp(− Tr[λ1 R∞ ] − Tr[2λ22 R∞ ]),

T →+∞

which gives another mean to prove Theorem 2.4. Here, we prove Theorem 3.4.
Proof of Theorem 3.4. Here, we focus on the case b = b0 Id with b0 > 0 and set T = e−b0 T .
Since the square root function is analytic on the set of positive definite matrices (see e.g.
[33], p. 134) we get that
p
2
ṽT = b0 Id − 2T λ2 + T (λ1 − 2λ22 ) + O(3T ),
b0
since the squares of each sides coincides up to a O(3T ) term. We observe that w̃T = 2T λ2 −b0 Id ,
√
2
and thus ṽT + w̃T = bT0 (λ1 − 2λ22 ) + O(3T ).
We now write

p −1  1
p
p
p
p
1
VvT ,wT (T ) =( ṽT )
exp( ṽT T )( ṽT + w̃T ) + exp(− ṽT T )( ṽT − w̃T )
2
2
p
p
p
p
1
1
Vv0T ,wT (T ) = exp( ṽT T )( ṽT + w̃T ) + exp(− ṽT T )(w̃T − ṽT ).
2
2
h
i
√
1
1
1
2
Since T exp( ṽT T ) → Id , we get T VvT ,wT (T ) →
b0 2b0 (λ1 − 2λ2 ) + b0 Id and
1
T

Vv0T ,wT (T )

T →+∞
1
(λ1 −
T →+∞ 2b0

→

T →+∞

2λ22 ) − b0 Id . This yields to

Vv0T ,wT (T )VvT ,wT (T )−1
We also have

exp − α
Tr[b0 Id ]T
2
α

det[VvT ,wT (T )] 2

+ b0 Id −→ (λ1 −


=

T →+∞

2λ22 )



1
(λ1 − 2λ22 ) + b0 Id
2b0

−1

1
1
h h
ii ,
→
α
1
2 T →+∞ det 1
det[−1
V
(T
)]
(λ
−2λ22 )+b0 Id
v
,w
1
T
T
T
b
2b
0

.

and there-

0

fore



−1 
exp − 2b10 Tr (λ1 − 2λ22 ) 2b12 (λ1 − 2λ22 ) + Id
x
0
h
i
lim E(T, λ1 , λ2 ) =
.
T →+∞
det 2b12 (λ1 − 2λ22 ) + Id

(69)

0



We now want to identify the limit. We know that X ∼ W ISd 2bx0 , α, 0, Id ; 4b12 has the
0
following Laplace transform

h
i
exp − Tr u(Id + 2b12 u)−1 2bx0
0
u ∈ Sd+ , E[exp(− Tr[uX])] =
.
det[Id + 2b12 u]
0

6 NUMERICAL STUDY

28

Let G̃ denote a d-square matrix independent from X, whose entries are independent and
follow a standard Normal distribution. By Lemma B.4, we have
√
√
E[exp(− Tr[λ1 X + λ2 ( X G̃ + G̃ X)])] = E[exp(− Tr[(λ1 − 2λ22 )X])].

Thus,
the
convergence in law of T (XT − x − αT Id − bRT − RT b), 2T RT to

 √ (69) shows
√
X, X G̃ + G̃ X under Pθ , which gives the claim of Theorem 3.4.

6

Numerical Study

In this section, we test the convergence of the MLE given by (15) and (18). To do so, we
consider a given large value of T and simulate the Wishart process exactly on the regular time
grid ti = iT
N , i = 0, · · · , N . This can be done by using the method presented in Ahdida and
Alfonsi [1], see also Alfonsi [3]. We take N sufficiently large and approximate the integrals RT
and Q−1
T applying the trapezoidal rule along this time grid. Thus, we will use the estimator
with the exact value of XT and these approximated values of RT and Q−1
T .
This section has three goals. First, we check numerically the convergence results that
we have obtained. Second, we investigate numerically the convergence of the MLE in some
nonergodic cases, where no theoretical result of convergence have been found. Last, we test
the estimation of the parameters of a full Wishart process (1). To do so, we estimate first a
with the quadratic variation and then the parameters α and b by using the MLE (15) on the
process (a> )−1 Xa−1 .

6.1

Numerical validation of the convergence results

Using the method mentioned above, we have checked the convergence results obtained in
this paper. Namely, we sample M = 10000 independent paths of X in order to draw an
histogram of the properly rescaled value of b̂i,j − bi,j or α̂ − α. We do not reproduce all these
graphics here, and present for example in Figure 1 an illustration of the convergence given
by Theorem 3.4.

5.0

7

4.5
6
4.0
5

3.5
3.0

4
2.5
3
2.0
1.5

2

1.0
1
0.5
0

0.0
-1.0

-0.8

-0.6

-0.4

-0.2

0.0

0.2

0.4

0.6

0.8

1.0

(a) Limit law of exp(0.05T )(b − b̂T )1,1 .

-1.0

-0.8

-0.6

-0.4

-0.2

0.0

0.2

0.4

0.6

0.8

1.0

(b) Limit law of exp(0.05T )(b − b̂T )1,2 .

0.5 0.1 ),
Figure 1: Asymptotic law of the error for the estimation of θ = b with for: x = ( 0.1
0.3
T = 100, N = 10000, α = 4.5 and b = 0.05Id .

6 NUMERICAL STUDY

6.2

29

Experimental convergence in a nonergodic case

In this paragraph, we try to guess the asymptotic behavior of the MLE in a nonergodic
case, where no theoretical convergence result is known. Namely, we observe in Figure 2
the asymptotic estimation error, when b = diag(0.1, 0.005) is diagonal with positive and
distinct terms on its diagonal and when we use the estimator (18). As one might have
guess, the convergence of the diagonal terms seems to be with an exponential rate, with the
exponential speed corresponding to its value. Namely, b̂11 seems to converge to b11 with a
speed of exp(0.1T ) while b̂22 seems to converge to b22 with a speed of exp(0.005T ). More
interesting is the antidiagonal term. One could have imagine that the convergence rate is
the slowest of these two rates. Instead, on our experiment, the convergence of b̂12 towards
b12 seems to happen with the rate exp(0.1T ). We have observed the same behaviour for
other parameter values. Of course, it would be hasty to draw a global conclusion from few
particular experiments. However, it is interesting to note that these numerical tests are a
way to guess or check the convergence rate of the MLE.

3.5

35

3.0

30

2.5

25

2.0

20

1.5

15

1.0

10

0.5

5

0.0
-1.5

-1.0

-0.5

0.0

0.5

1.0

0
-0.12

1.5

(a) Limit law of exp(0.1T )(b − b̂T )1,1 .

-0.10

-0.08

-0.06

-0.04

-0.02

0.00

0.02

0.04

0.06

(b) Limit law of exp(0.005T )(b − b̂T )2,2 .

3.5

3.0

2.5

2.0

1.5

1.0

0.5

0.0
-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

(c) Limit law of exp(0.1T )(b − b̂T )1,2 .
0.3 0.1 ), T = 100,
Figure 2: Asymptotic law of the error for the estimation of θ = b with x =( 0.1
0.2
N = 10000, α = 3.5 and b = diag(0.1, 0.005).

6.3

Estimation of the whole Wishart process

In this last part of the numerical study, we perform the estimation of all the parameters of the
Wishart process (1). We consider a case where a is upper triangular and (a> )−1 ba> is symmetric. We proceed as follows. First, we sample exactly a discrete path (XiT /N , 0 6 i 6 N ).

6 NUMERICAL STUDY

30

Then, we estimate the matrix a> a by using (3), where the quadratic variations are replaced
by their classical approximations and the integrals are replaced by the trapezoidal rule. By
a Cholesky decomposition we get then an estimator â of a. Then, we use the MLE (15) on
the path ((â> )−1 XiT /N â> , 0 6 i 6 N ). This gives an estimator of α and (a> )−1 ba> , and
therefore an estimator of b. As a comparison, we also calculate similarly the estimator of
α and b when a is known and has not to be estimated. To draw histograms or calculate
empirical expectations, we run M = 10000 independent paths of X.
We consider a sufficiently large value of T and are interested in looking at the convergence
with respect to N . First, we plot the the error on the estimator of a with respect to the
number of time step in Log-Log scale. We observe that the convergence to zero takes place
with experimental rate close to 1/2. This is in line with the general results on the estimation
of the diffusion coefficient, see Dohnal [12] and Genon-Catalot and Jacod [15]. Then, we
focus on the influence of the discretization and the unknown parameter a on the convergence
of the MLE of b and α. In Table 1, we give in function of N the Mean Squared Error
MSE(θbN |θ) = E[|θbN − θ|2 ] of the estimator θbN , with θ = (b, α). It is estimated with the
empirical expectation. First, we observe that the convergence of the estimator of α is roughly
the same whether we know a or not. This is expected since the estimation of α does not
depend on the estimation of a. Instead, the bias on b is much higher when a is estimated
than when a is known. However, it decreases also faster at an experimental order of 0.7 while
the bias when a is known decreases at an experimental order of 0.45. This latter rate is in
line with the rate of 1/2 obtained in dimension 1 by Ben Alaya and Kebaier [5]. In our case,
it seems that the influence of the estimation of a vanishes
around N = 5000. Last, we have
√
N
b
plotted in Figure 4 the limit law of the estimator T (θ − θ) with N = 10000.
This short numerical study shows that the estimator obtained by discretizing the continuous time estimator is efficient in practice. Of course, it would be nice to obtain general
convergence results in function of T and N , but we leave this for further research.

1

10

0

10

-1

10

-2

10

1

10

2

10

3

10

4

10

Figure 3: Log-Log representation of the empirical expectation
of E[Tr[(a − b
aN )2 ]]1/2 for

−1 0.2
0.8
0.5
1
1
x =( 0.5 1 ), T = 100, a =( 0 2 ), α = 4.5, b = 2 −2 , where the line is the simple linear
regression i.e. log(E[Tr[(a − b
aN )2 ]]1/2 ) ≈ 2.62 − 0.58 log(N ).

6 NUMERICAL STUDY

31

Number of time steps
E[Tr[(a

−b
aN )2 ]]1/2

b
a=a
b
a=b
aN
b
a=a
b
a=b
aN
b
a=a
b
a=b
aN
b
a=a
b
a=b
aN

MSE(bbN
1,1 |b1,1 )
MSE(bbN
2,2 |b2,2 )
MSE(bbN
1,2 |b1,2 )
MSE(b
αN |α)

20

50

100

200

500

1000

2000

1.7671
0.0745
0.7636
0.2554
3.4085
0.0075
0.0442
0.8448
0.8267

1.4311
0.0338
0.5266
0.1310
2.8722
0.0033
0.0568
0.3579
0.3496

1.1487
0.0181
0.3489
0.0664
2.1159
0.0017
0.0596
0.1993
0.1895

0.7913
0.0115
0.1891
0.0372
1.1995
0.0011
0.0352
0.1151
0.1095

0.4107
0.0082
0.0624
0.0231
0.3600
0.0008
0.0148
0.0614
0.0617

0.2472
0.0069
0.0273
0.0176
0.1264
0.0008
0.0075
0.0416
0.0410

0.1514
0.0061
0.0142
0.0153
0.0480
0.0007
0.0039
0.0308
0.0311

5000
0.0846
0.0058
0.0085
0.0139
0.0201
0.0007
0.0019
0.0230
0.0234

Table 1: Mean Squared Error for the estimation of θ = (b, α) with respect to N . Same
parameters as Figure 3.

0.35

0.50

0.30
0.45

0.25

0.40

0.35

0.20
0.30

0.25

0.15

0.20

0.10
0.15
0.10

0.05

0.05

0.00

0.00
-5

-4

-3

-2

-1

0

(a) Limit law of

1

√

2

3

-6

4

T (b − b̂N
T )1,1 .

-4

-2

0

2

(b) Limit law of

√

4

6

8

T (b − b̂N
T )2,2 .

0.30

1.0

0.25
0.9
0.8

0.20
0.7
0.6

0.15
0.5
0.4

0.10

0.3

0.05

0.2

0.1

0.00

0.0
-2.0

-1.5

-1.0

-0.5

0.0

(c) Limit law of

0.5

√

1.0

1.5

T (b − b̂N
T )1,2 .

2.0

-8

-6

-4

-2

0

(d) Limit law of

2

√

4

6

8

T (α − α̂TN ).

Figure 4: Asymptotic laws of the error for the estimation of θ = (b, α) for b
a=b
aN , N = 10000,
same parameters as Figure 3.

A PROOF OF PROPOSITION 1.1

A

32

Proof of Proposition 1.1

We denote bs = (b + b> )/2 (resp. ba = (b − b> )/2) the symmetric (resp. antisymmetric) part
of b. We have


Z T
Z
p
1
det[XT ]
1 T
−1
Tr[( Xs ) dWs ] = log
− Tr[b]T −
(α − 1 − d) Tr[Xs−1 ]ds,
2
det[x]
2 0
0
Z T
Z
p
p
p
1 T
s
Tr[b Xs dWs ] =
Tr[bs ( Xs dWs + dWs> Xs )]
2 0
0
Z
Tr[bs XT ] − Tr[bs x] αT
1 T
s
Tr[bs (bXt + Xt b> )]dt.
=
−
Tr[b ] −
2
2
2 0


RT
Rt
√
√
a
X
Thus, the only part to calculate is E exp( 0 Tr[−b Xs dWs ]) FT , and we set Mta = 0 Tr[ba Xs dWs ].
We now observe that hTr[As dWs ], Tr[Bs dWs ]i = Tr[As Bs> ]ds and are looking for the process
Γ that takes values in Sd and minimizes

hdMta −Tr[Γt (dXt −(αId +bXt +Xt b> )dt)]i = − Tr[ba Xt ba ] + 2 Tr[Γt (Xt ba − ba Xt )] + 4 Tr[Γ2t Xt ] dt.
We obtain that 2(Xt ba − ba Xt ) + 4(Xt Γt + Γt Xt ) = 0 and thus


1 a
a
Γt = L−1
(b
X
−
X
b
)
.
t
t
Xt
2
It satisfies Tr[Γt (Xt ba − ba Xt )] = −2 Tr[Γt (Γt Xt + Xt Γt )] = −4 Tr[Γ2t Xt ]. By construction,
we have hdMta − Tr[Γt (dXt − (αId + bXt + Xt b> )dt)], Tr[Γ̃(dXt − (αId + bXt + Xt b> )dt)]i = 0
for any Γ̃ ∈ Sd . Thus, there exists a Brownian
p motion β independent of X such that
dMta − Tr[Γt (dXt − (αId + bXt + Xt b> )dt)] = − Tr[ba Xt ba ] − Tr[Γt (ba Xt − Xt ba )]dβt . In
Rtp
fact, both processes (Xt , 0 − Tr[ba Xs ba ] − Tr[Γt (ba Xs − Xs ba )]dβs ) and
(Xt , Mta −

Z

t

p
p
Tr[Γs ( Xs dWs + dWs> Xs )])

0

solve the same martingale problem for which uniqueness holds. Therefore, we have


 Z T

p
E exp −
Tr[ba Xt dWt ] FTX
0

Z T
h
i 1Z T
a
a
a
a
>
= exp
Tr Γt ((αId + bXt + Xt b )dt − dXt ) −
Tr[b Xt b ] + Tr[Γt (b Xt − Xt b )]dt
2 0
0
 Z T



 Z T
1
−1 1 a
a
= exp −
Tr LXt
(b Xt − Xt b ) dXt −
Tr[ba Xt ba ]dt
2
2
0
0





Z
Z T
1 T
1
−1
a
a
a
a
a
s
+
Tr LXt
(b Xt − Xt b ) (b Xt − Xt b ) dt +
Tr[b Xt b ]dt ,
2 0
2
0

 1
1 a
a
since Tr L−1
= 2 Tr[Xt−1 (ba Xt − Xt ba )] = 0 by Lemma B.1 and
Xt 2 (b Xt − Xt b )
Tr[Γt (bs Xt + Xt bs )] = Tr[bs (Γt Xt + Xt Γt )] =

1
Tr [bs (ba Xt − Xt ba )] = Tr[ba Xt bs ].
2

B TECHNICAL LEMMAS

33

Using (6) and the previous calculations, we obtain


Z
α − α
Tr[bs XT ] − Tr[bs x] 1 T
det[XT ]
0
= exp
+
Tr[(bs )2 Xs ]ds
log
−
4
det[x]
2
2 0
Z T
Z
 T
α − α0 α + α0
αT
a
s
Tr[b Xs b ]ds −
−
Tr[Xs−1 ]ds −
−1−d
Tr[b]
4
2
2
0
0
Z
Z

 1 T
 −1 a
 
1 T
a
a
a
a
a
Tr L−1
(b
X
−
X
b
)
dX
−
Tr
L
(b
X
−
X
b
)
(b
X
−
X
b
)
dt .
+
t
t
t
t
t
t
t
Xt
Xt
2 0
4 0

0
Lθ,θ
T

−1 a
s
s
s
a
s
s
a
s
Last, we use L−1
Xt (b Xt + Xt b ) = b and Tr[LXt (b Xt − Xt b ) (b Xt + Xt b )] = 2 Tr[b Xt b ]
to obtain (7).

B

Technical lemmas

Lemma B.1. For X ∈ Sd+,∗ and a > 0, let LX,a and LX = LX,0 be the linear applications defined by (12) on Sd . If a Tr[X −1 ] 6= 1, then LX,a is invertible and we have
Tr[X −1 Y ]
−1
Tr[L−1
X,a (Y )] = 2(1−a Tr[X −1 ]) . Besides, the map (X, Y, a) 7→ LX,a (Y ) is continuous on
{(X, Y, a) ∈ Sd+,∗ × Sd × R+ , a Tr[X −1 ] 6= 1}.
Proof. The invertibility of LX,a is equivalent to its one-to-one property. Since X ∈ Sd+,∗ ,
there exists an orthogonal matrix OX and a diagonal matrix DX with positive elements such
> . We get
that X = OX DX OX
>
>
Y ∈ ker(LX,a ) ⇐⇒ OX DX OX
Y + Y OX DX OX
= 2a Tr[Y ]Id
>
>
⇐⇒ DX (OX
Y OX ) = 2a Tr[Y ]Id − (OX
Y OX )DX .

(70)


> Y O )D
>
Since DX is diagonal, we obtain for 1 6 i, k 6 d, (OX
X
X i,k = (OX Y OX )i,k (DX )k,k

>Y O )
>
>
and DX (OX
X i,k = (DX )i,i (OX Y OX )i,k . For k 6= i, (70) gives (OX Y OX )i,k = 0. For
> Y O ) (D ) = a Tr[Y ] and therefore
k = i, we get (OX
X i,i
X i,i
>
Tr[Y ] = Tr[OX
Y OX ] = Tr[Y ]a

d
X

1

i=1

(DX )i,i

= Tr[Y ]a Tr[X −1 ].

> Y O ) = 0, which gives Y = 0 and
Since a Tr[X −1 ] 6= 1, we obtain Tr[Y ] = 0 and then (OX
X i,i
−1
the invertibility of LX,a . Let c = LX,a (Y ). We have c + X −1 cX − 2a Tr[c]X −1 = X −1 Y ,
which gives 2(1 − a Tr[X −1 ]) Tr[c] = Tr[X −1 Y ]. Last, the continuity property is obvious since
(X, a) 7→ LX,a is continuous and L 7→ L−1 is continuous on {L : Sd → Sd linear and invertible}.

Lemma B.2. For X ∈ Sd+,∗ , LX is self-adjoint and positive definite:
Tr[LX (Y )Y ] ≥ 2λ(X) Tr[Y 2 ],
where λ(X) > 0 is the lowest eigenvalue of X. Besides, for a < 1/ Tr[X −1 ], LX,a is selfadjoint and positive definite.

B TECHNICAL LEMMAS

34

Proof. For Y, Z ∈ Sd , we have Tr[LX (Y )Z] = Tr[(XY +Y X)Z] = Tr[Y (XZ+ZX)] = Tr[Y LX (Z)]
and Tr[LX (Y )Y ] = 2 Tr[XY 2 ] ≥ 2λ(X) Tr[Y 2 ] since X − λ(X)Id ∈ Sd+ . The self-adjoint
property is then clear for LX,a , and the positive definiteness comes from Lemma B.1 and the
continuity of the eigenvalues of LX,a with respect to a.
>
Lemma B.3. For X ∈ Sd+,∗ , Y ∈ Md , L̄X (Y ) = L−1
X (Y X + XY )X is self-adjoint and
−1
>
positive. The linear application L̄X,a (Y ) = LX (Y X + XY )X − a Tr[Y ]Id is also positive
for a < 1/ Tr[X −1 ], and there is a positive cX,a > 0 such that
> 2
Tr[L̄X (Y )> Y ] ≥ cX,a Tr[(L−1
X (Y X + XY )) ].

Proof. Since L−1
X is self-adjoint, we have for Z ∈ Md
>
Tr[L̄X (Y )> Z] = Tr[L−1
X (Y X + XY )ZX] =

=

1
>
>
Tr[L−1
X (Y X + XY )(ZX + XZ )]
2

1
>
>
Tr[(Y X + XY > )L−1
X (ZX + XZ )] = Tr[Y L̄X (Z)].
2

Similarly, Tr[L̄X,a (Y )> Z] = Tr[L̄X (Y )> Z] − a Tr[Y ] Tr[Z] = Tr[Y > L̄X,a (Z)]. Besides, we
notice that

1
>
2
Tr[L̄X,a (Y )> Y ] =
(Y
X
+
XY
)]
−
2a
Tr[Y
]
Tr[(Y X + XY > )L−1
X
2


1
−1
>
>
=
Tr[LX,a (L−1
X (Y X + XY ))LX (Y X + XY )]
2
by Lemma B.1. This gives the claim since LX,a is positive definite by Lemma B.2.
The following lemma gives the Laplace transform of the matrix Normal distribution.
Lemma B.4. Let C ∈ Sd+,∗ and C[C] ∈ (Rd )⊗4 defined by
C[C]i,j,k,l = δik Cj,l + δil Cj,k + δjk Ci,l + δjl Ci,k .

(71)

We introduce the Md -valued random variables G̃ and G ∼ N (0, C[C]) of which components are Normal random variables with mean 0 such that

∀i, j, k, l ∈ {1, . . . , d},

E[G̃i,j G̃k,l ] = δik δjl ,

E[Gi,j Gk,l ] = C[C]i,j,k,l .

(72)

We have the following results.


1. For all c ∈ Sd , E exp(− Tr[cG]) = exp(2 Tr[c2 C]).
2. For C̃ ∈ Md such that C̃ C̃ > = C, C̃ G̃ + G̃> C̃ > and G have the same law.
√
√
> X)])] = E[exp(Tr[cL−1 (c)])
3. Let X ∈ Sd+,∗ . For c ∈ Sd , E[exp(− Tr[cL−1
X ( X G̃ + G̃
X

C SOME ASYMPTOTIC BEHAVIOUR OF WISHART PROCESSES

35

Proof. We focus on the first point. For all c ∈ Sd , we have
X




E exp(− Tr[cG]) = E exp(−
ci,j Gi,j ) .
16i,j6d

Moreover,

P

16i,j6d ci,j Gi,j


E

X

is a Normal random variable and its variance is given by

ci,j Gi,j

2 

=

16i,j6d

X

ci,j ck,l C[C]i,j,k,l = 4 Tr[c2 C].

16i,j,k,l6d

It follows from the moment generating function of the Normal distribution that


E exp(− Tr[cG]) = exp(2 Tr[c2 C]).
To prove the second point it is sufficient to notice that Tr[c(C̃ G̃ + G̃> C̃ > )] = Tr[2cC̃ G̃]
and
X
X
 X
2 
E
(cC̃)i,j G̃i,j
=
(cC̃)i,j (cC̃)k,l δik δjl =
(cC̃)2i,j = Tr[cC̃ C̃ > c].
16i,j6d

16i,j,k,l6d

16i,j6d

√
√
√
√
> X) and have XZ+ZX =
For the third point, we set Z = L−1
X G̃+G̃> X.
X ( X G̃+G̃
We also introduce c̃ = L−1
X (c) and have c̃X + X c̃ = c. Thus, we obtain
√
√
Tr[cZ] = Tr[(c̃X + X c̃)Z] = Tr[c̃( X G̃ + G̃> X)]
and therefore E[exp(− Tr[cZ])] = exp(2 Tr[c̃2 X]) = exp(Tr[c̃(c̃X + X c̃)]) = exp(Tr[c̃c]).

C

Some asymptotic behaviour of Wishart processes

Lemma C.1. Let X ∼ W ISd (x, α, b, Id ) with b ∈ Sd , x ∈ Sd+ and α > d − 1. Then XT
converges in law √
when T → +∞ if and only if −b ∈ Sd+,∗ . In this case, XT converges in law
to W ISd (0, α, 0, −b−1 ; 1/2).
Let X ∼ W ISd (x, α, b, Id ) with b ∈ Md , x ∈ Sd+ and α > d − 1. If −(b + b> ) ∈ Sd+,∗ ,
R∞
√
>
q∞ := 0 esb esb ds is well defined and XT converges in law to W ISd (0, α, 0, 2q∞ ; 1/2).
Proof. Let us first consider the case −b ∈ Sd+,∗ . From Proposition 4 in [2], we have for
v ∈ Sd+ ,
  

R
 −1
T
exp Tr −v Id + 2 0 e2bs ds v
eT b xeT b
E[exp(− Tr[vXT ])] =
h
R
 iα/2
T
det Id + 2 0 e2bs ds v
→

1

,
det [Id − b−1 v]α/2
√
which is the Laplace transform of W ISd (0, α, 0, −b−1 ; 1/2). Now, let us consider −b 6∈ Sd+,∗ .
Then, there exists an eigenvector v ∈ Rd \ {0} such that bv = λv with λ > 0. Then, we have
d
>
>
>
>
→ +∞.
dt E[v Xt v] = αv v + 2λE[v Xt v], and therefore E[v XT v]
T →+∞

T →+∞

C SOME ASYMPTOTIC BEHAVIOUR OF WISHART PROCESSES

36

In the case b ∈ Md with −(b + b> ) ∈ Sd+,∗ , we know that the norm of ebs decays
exponentially to 0 as s → +∞, see e.g. Problem 11.3.6 in Golub and Van Loan [18]. Using
1
again Proposition 4 in [2], we get that E[exp(− Tr[vXT ])] →
α/2 .
T →+∞ det[Id +2q∞ v]

Lemma C.2. • Assume α > d + 1 and b = 0. Then,
ZT
log(T )

Q−1
1
T
→ α−(d+1)
d log(T ) T →+∞

a.s. Besides,

converges almost surely to d, and we have
"
∀µ > 0, sup E exp
T >2

!#

µ

p
NT
log(T )

• Assume α = d + 1 and b = 0. Then, as T → +∞,



< ∞.

(73)

2

Q−1
T converges in law to

2
d log(T )

τ1 = inf{t > 0, Bt = 1}, where B is a Brownian motion. Besides,
in probability to d, and we have



µ
∀µ > 0, sup E exp
NT
< ∞.
log(T )
T >2

ZT
log(T )

=

2NT
log(T )

converges

(74)

We mention that the results on the convergence for QT are given in Donati-Martin et
al. [13]. However, their proofs is in a working paper by the same authors that we have not
been able to find. For this reason, we present here an autonomous proof.
√
√
Proof. We first consider the case α > d + 1. We have dXt = αId dt + Xt dWt + dWt> Xt
and thus
q
q
d(e−t Xet −1 ) = [αId − e−t Xet −1 ]dt + e−t Xet −1 dW̃t + dW̃t> e−t Xet −1 ,
with dW̃t = e−t/2 d(Wet −1 ). We observe that W̃ is a matrix Brownian motion, which gives
Y ∼ W ISd (x, α, −Id /2, Id ), where Yt = e−t Xet −1 for t > 0. Using equation (25) to the
process Y , we get


Z
Z
q
1
det[Yt ]
1 t
2 t
−1
log
= (α − 1 − d)
Tr[Ys ]ds − d +
Tr[ Ys−1 dW̃s ].
(75)
t
det[Y0 ]
t 0
t 0
p
Rt
Rt
Since Y is ergodic and h 0 Tr[ Ys−1 dW̃s ]i = 0 Tr[Ys−1 ]ds, we get that the left hand side con−1 ]]−d,
verges in probability to zero and the right hand side converges a.s. to (α−1−d)E[Tr[Y∞

√
det[Yt ]
1
where Y∞ ∼ W ISd (0, α, 0, 2Id ; 1/2) is the stationary law of Y . Therefore, t log det[Y0 ]




 det[e−t X
 det[X
det[Yt ]
1
1
et −1 ]
et −1 ]
converges a.s. to zero. Since 1t log det[Y
=
log
=
log
− d,
t
t
det[x]
det[x]

 0]
det[XT ]
ZT
1
we get that log(T
converges a.s. to d when T → +∞.
) = log(T ) log
det[x]
Now, we use (25) taken at time T = et − 1 and Dubins-Schwarz theorem: there is a
Brownian motion β such that for all t > 0,
α − (1 + d)
+
Qet −1 t
This gives that

α−(1+d)
→
Qet −1 t t→+∞

2βQ−1

et −1

t

1
= log
t

d a.s., and therefore



det[Xet −1 ]
det[x]


.

Q−1
1
T
→ α−(d+1)
,
d log(T ) T →+∞

a.s.

C SOME ASYMPTOTIC BEHAVIOUR OF WISHART PROCESSES

37

ZT
to prove
(25), wehave NT = Z2T − α−1−d
Q−1
T ≤ 2 and thus
2
 (73).
 From 
 It remains

µ
√
T ] 2 log(T )
6 E det[X
< ∞, since the moments of X are bounded.
E exp √ µ NT
det[x]
log(T )

Again we set t = log(T + 1), and for Λ ∈ [0, 1], we have from (75)
Z T
Z t
q
q
−1
NT =
Tr[ Xs dWs ] =
Tr[ Ys−1 dW̃s ]
0
0




Z
Z t
q
d
det[Yt ]
α−1−d t
1
−1
−1
+ t−
Tr[Ys ]ds .
Tr[ Ys dW̃s ] + (1 − Λ)
log
=Λ
2
det[x]
2
2
0
0
By Cauchy-Schwarz inequality, we get
"
!#
µ
E exp p
NT
log(T + 1)
"
µ #

√
µd(1−Λ) √
1
det[Yt ] (1−Λ) t
t
E2
≤e 2
det[x]



Z
Z
q
1
2µΛ t
α−1−d t
−1
−1
2
√
Tr[ Ys dW̃s ] − µ(1 − Λ)
Tr[Ys ]ds .
× E exp √
t 0
t
0

√
2µ √
We now take Λ = Λt = 21t −1 + 1 + 4t with t = (α−1−d)
in order to obtain
t

2
1 2µΛ
√ t
√
= µ(1 − Λt ) α−1−d
. We note that for t large enough, Λt ∈ [0, 1]. Besides, we
2
t
t
√
2µ
. From Theorem 4.1
have Λt = 1 − t + o(1/t), so that t(1 − Λt ) converges to α−1−d
t→+∞

in [29], the second expectation is then equal to 1, while the first one is bounded since Y is
ergodic. This yields to (73).
We now consider the case α = d + 1. We set again t = log(1 + T ) and have T = et − 1.
Thus,






det[XT ]
det[et Yt ]
det[Yt ]
ZT = log
= log
= log
+ dt.
det[x]
det[x]
det[x]
√
Again, Yt converges
in
law
towards
W
IS
(0,
α,
0,
2Id ; 1/2). Therefore, the ergodic theorem
d


det[Yt ]
1
gives that t log det[x] converges in probability to 0, which yields to the convergence in

2
ZT
2
probability of log(T
to
d.
We
now
turn
to
the
convergence
of
Q−1
T . We know from
)
d log(T )
Theorem 4.1 in Mayerhofer [29] that for T > 0 and λ > 0,



2λ
(2λ)2
−1
E exp
NT − 2
Q
= 1.
d log(1 + T )
2d log(1 + T )2 T
From (25), we have NT = ZT /2 and we write



(2λ)2
−1
1 =E exp λ − 2
Q
2d log(1 + T )2 T






(2λ)2
2λ
−1
Q
exp
NT − exp(λ)
+ E exp − 2
2d log(1 + T )2 T
d log(1 + T )


(2λ)2
−1
We now observe that exp − 2d2 log(1+T
Q
6 1 and that
)2 T
"
#
"
λ




λ#
dt
2λ
det[XT ] d log(1+T )
det[Y
]
t
E exp
NT
=E
= eλ E
.
d log(1 + T )
det[x]
det[x]

REFERENCES

Since Y has bounded moments and is stationary, supt>1 E

38


det[Yt ]
det[x]

λ
dt

< ∞. This gives the

uniform integrability (74) and that






2λ
(2λ)2
−1
Q
exp
→ 0.
E exp − 2
NT − exp(λ)
T →+∞
2d log(1 + T )2 T
d log(1 + T )
h

i
(2λ)2
−1
Therefore, limT →+∞ E exp λ − 2d2 log(1+T
Q
= 1, which gives the desired conver2
T
)
gence in law.
Acknowledgements. The authors would like to thank Arnaud Gloter (University of Evry)
and Marina Kleptsyna (University of Le Mans) for helpful discussions, and the two anonymous referees for their fruitful comments.

References
[1] A. Ahdida and A. Alfonsi. Exact and high-order discretization schemes for Wishart
processes and their affine extensions. Ann. Appl. Probab., 23(3):1025–1073, 2013.
[2] A. Ahdida, A. Alfonsi, and E. Palidda. Smile with the Gaussian term structure model.
Working paper, 2014.
[3] A. Alfonsi. Affine diffusions and related processes: simulation, theory and applications,
volume 6 of Bocconi & Springer Series. Springer, Cham; Bocconi University Press,
Milan, 2015.
[4] M. Ben Alaya and A. Kebaier. Parameter estimation for the square-root diffusions:
ergodic and nonergodic cases. Stoch. Models, 28(4):609–634, 2012.
[5] M. Ben Alaya and A. Kebaier. Asymptotic behavior of the maximum likelihood estimator
for ergodic and nonergodic square-root diffusions. Stoch. Anal. Appl., 31(4):552–573,
2013.
[6] M. Bru. Thèse 3ème cycle. Résistance d’Escherichia coli aux antibiotiques: Sensibilités
des analyses en composantes principales aux perturbations Browniennes et simulation.
PhD thesis, Université Paris Nord, 1987.
[7] M. Bru. Wishart processes. J. Theoret. Probab., 4(4):725–751, 1991.
[8] C. Cuchiero, D. Filipović, E. Mayerhofer, and J. Teichmann. Affine processes on positive
semidefinite matrices. Ann. Appl. Probab., 21(2):397–463, 2011.
[9] J. Da Fonseca, M. Grasselli, and F. Ielpo. Estimating the Wishart affine stochastic
correlation model using the empirical characteristic function. Stud. Nonlinear Dyn.
Econom., 18(3):253–289, 2014.
[10] J. Da Fonseca, M. Grasselli, and C. Tebaldi. Option pricing when correlations are
stochastic: an analytical framework. Review of Derivatives Research, 10:151–180, 2008.
[11] L. Dieci and T. Eirola. Positive definiteness in the numerical solution of Riccati differential equations. Numer. Math., 67(3):303–313, 1994.

REFERENCES

39

[12] G. Dohnal. On estimating the diffusion coefficient. J. Appl. Probab., 24(1):105–114,
1987.
[13] C. Donati-Martin, Y. Doumerc, H. Matsumoto, and M. Yor. Some properties of the
Wishart processes and a matrix extension of the Hartman-Watson laws. Publ. Res.
Inst. Math. Sci., 40(4):1385–1412, 2004.
[14] E. Fournié and D. Talay. Application de la statistique des diffusions à un modèle de
taux d’interêt. Finance, 12(2):79–111, 1991.
[15] V. Genon-Catalot and J. Jacod. On the estimation of the diffusion coefficient for multidimensional diffusion processes. Ann. Inst. H. Poincaré Probab. Statist., 29(1):119–151,
1993.
[16] A. Gnoatto. The Wishart short rate model. Int. J. Theor. Appl. Finance, 15(8):1250056,
24, 2012.
[17] A. Gnoatto and M. Grasselli. The explicit Laplace transform for the Wishart process.
J. Appl. Probab., 51(3):640–656, 2014.
[18] G. H. Golub and C. F. Van Loan. Matrix computations. Johns Hopkins Studies in the
Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD, third edition,
1996.
[19] C. Gourieroux and R. Sufana. Derivative pricing with Wishart multivariate stochastic
volatility. J. Bus. Econom. Statist., 28(3):438–451, 2010.
[20] C. Gourieroux and R. Sufana. Discrete time Wishart term structure models. J. Econom.
Dynam. Control, 35(6):815–824, 2011.
[21] S. Heston. A closed-form solution for options with stochastic volatility with applications
to bond and currency options. The Review of Financial Studies, 6(2):327–343, 1993.
[22] P. Jeganathan. Some aspects of asymptotic theory with applications to time series
models. Econometric Theory, 11(5):818–887, 1995. Trending multiple time series (New
Haven, CT, 1993).
[23] Y. A. Kutoyants. Statistical inference for ergodic diffusion processes. Springer Series in
Statistics. Springer-Verlag London, Ltd., London, 2004.
[24] L. Le Cam. Locally asymptotically normal families of distributions. Certain approximations to families of distributions and their use in the theory of estimation and testing
hypotheses. Univ. california Publ. Statist., 3:37–98, 1960.
[25] L. Le Cam and G. L. Yang. Asymptotics in statistics. Springer Series in Statistics.
Springer-Verlag, New York, second edition, 2000. Some basic concepts.
[26] J. J. Levin. On the matrix Riccati equation. Proc. Amer. Math. Soc., 10:519–524, 1959.
[27] R. S. Liptser and A. N. Shiryaev. Statistics of random processes. I,II. Applications of
Mathematics (New York). Springer-Verlag, Berlin, expanded edition, 2001. Applications,
Translated from the 1974 Russian original by A. B. Aries, Stochastic Modelling and
Applied Probability.

REFERENCES

40

[28] H. Luschgy. Local asymptotic mixed normality for semimartingale experiments. Probab.
Theory Related Fields, 92(2):151–176, 1992.
[29] E. Mayerhofer. Wishart Processes and Wishart Distributions: An Affine Processes Point
of View. ArXiv e-prints, Jan. 2012.
[30] B. Mond and J. E. Pečarić. On matrix convexity of the Moore-Penrose inverse. Internat.
J. Math. Math. Sci., 19(4):707–710, 1996.
[31] L. Overbeck. Estimation for continuous branching processes.
25(1):111–126, 1998.

Scand. J. Statist.,

[32] G. Pagès. Sur quelques algorithmes récursifs pour les probabilités numériques. ESAIM
Probab. Statist., 5:141–170 (electronic), 2001.
[33] L. C. G. Rogers and D. Williams. Diffusions, Markov processes, and martingales. Vol.
2. Cambridge Mathematical Library. Cambridge University Press, Cambridge, 2000. Itô
calculus, Reprint of the second (1994) edition.
[34] T. H. Rydberg. A note on the existence of unique equivalent martingale measures in a
markovian setting. Finance and Stochastics, 1(3):251–257, 1997.
[35] D. W. Stroock. Probability theory, an analytic view. Cambridge University Press, Cambridge, 1993.

