Oliver Huth

Felix Rüting

Svend-Age Biehs

Martin Holthaus

Shape-dependence of near-field heat transfer between a spheroidal nanoparticle and a flat surface

Keywords: . 44.40.+a Thermal radiation -78.67.Bf Nanocrystals, nanoparticles, and nanoclusters -41.20.Jb Electromagnetic wave propagation; radiowave propagation

We study the radiative heat transfer between a spheroidal metallic nanoparticle and a planar metallic sample for near-and far-field distances. In particular, we investigate the shape dependence of the heat transfer in the near-field regime. In comparison with spherical particles, the heat transfer typically varies by factors between 1/2 and 2 when the particle is deformed such that its volume is kept constant. These estimates help to quantify the deviation of the actual heat transfer recorded by a near-field scanning thermal microscope from the value provided by a dipole model which assumes a perfectly spherical sensor.

Introduction

The optical properties of metallic nanoparticles depend significantly on their shapes, as has been demonstrated, e.g., for elliptical gold particles attached to the apex of fiber-based probes for near-field optical microscopy [1]. Similarly, the thermal near-field radiation exchanged between a nanometer-sized particle at temperature T P and a closely spaced surface at different temperature T S is influenced by the particle's shape, although its linear extension may be significantly smaller than the dominant thermal wavelength. This shape-dependence occurs even when the particle's volume, and thus the total amount of radiating matter, is kept constant. In this paper we quantify this effect for the case of spheroidal metallic nanoparticles.

Thermally induced near fields have attracted much experimental and theoretical attention in last years [2,3]. Several experiments and models have been designed in order to measure and describe the radiative heat transfer in the near-field regime, i.e., for distances smaller than the thermal wavelength. In this regime tunneling of thermal photons leads to a magnitude of heat transfer which can exceed that achieved by black-body far-field radiation by several orders of magnitude. Possible applications of this phenomenon include, in particular, thermophotovoltaic devices [4,5,6,7,8].

The near-field radiative heat transfer between two dielectric bodies has been calculated within the framework of Rytov's fluctuational electrodynamics [START_REF] Rytov | Principles of Statistical Radiation[END_REF] for various geometries, including two semi-infinite planar bodies and a Present address: Laboratoire Charles Fabry, Institut d'Optique, CNRS, Université Paris-Sud, Campus Polytechnique, RD128, 91127 Palaiseau cedex, France layered structures (see, e.g., [START_REF] Polder | [END_REF]11,12,13,14,15,16,17]), two spheres (e.g., [18,19,20]), a sphere above a semi-infinite body (e.g. [13,21,22,23,24,25]), and certain other two-dimensional geometries [26]. Early experiments to detect the radiative heat transfer between two effectively semiinfinite bodies with flat surfaces have been performed by Hargreaves [27], and by Domoto et al. [28] Also a pioneering, but unconclusive experiment [29] by Xu et al. needs to be mentioned. An accurate measurement using glass plates separated by a micron-sized gap has been reported only recently by Hu et al. [30] Moreover, there now exist experimental setups measuring the near-field radiative heat transfer between a sphere and a flat sample [31,32,33], involving spheres with radii of some 10 µm.

Relying on the same basic principle, a near-field scanning thermal microscope (NSThM) has been developed by Kittel et al. [34,35,36] for recording the radiative heat transfer at probe-sample distances even down to some nanometers. The sensor of this device consists of the tip of a scanning tunneling microscope, functionalized to act as a thermocouple, so that here the sensor-sample geometry differs significiantly from those geometries for which the radiative heat transfer has been calculated exactly. The foremost part of such a tip typically has a radius of less than 50 nm.

If one regards this active part of the NSThM sensor as a sphere, one may employ the familiar dipole model (see, e.g., [13,21,22,23,24]) for describing the results obtained with such an instrument quantitatively. However, the actual shape of the sensor is somewhat prolonged, resembling more an ellipsoid than a sphere, and varies from specimen to specimen, due to the difficult production process. Such variations of the geometry will have consequences for the signal measured with the NSThM. This is what motivates our present investigation of the shape-dependence of the near-field radiative heat transfer: In this work, we study the change of the near-field radiative heat transfer in response to shape variations of ellipsoidal nanoparticles above a flat sample. Besides helping one to estimate errors implied by the use of the dipole model, there may also be other, more general nanotechnological applications.

The paper is organized as follows: In Sec. 2 we briefly review the dipole model for calculating the near-field heat transfer between a dielectric sphere and a planar sample. This step mainly serves to collect the required material in a form which can be easily generalized. In Sec. 3 we extend this dipole model to the geometry of an ellipsoidal particle above a flat sample. The near-field heat transfer for this geometry is then calculated explicitly in Sec. 4. In Sec. 5 we discuss the shape dependence of the heat transfer as the axes of the ellipsoid are varied, and compare our results to the heat transfer between a flat sample and a sphere.

Dipole model

The near-field heat transfer between a spherical particle with temperature T P and a sample with a planar surface and temperature T S can be estimated by means of a simplified dipole model developed previously by several authors (e.g., [2,3,13,21,22,23,24]). To this end, one first considers the thermally fluctuating electric and magnetic fields E f and H f outside the sample, being generated by thermally fluctuating charges in its interior. These fluctuating fields consist of a propagating and of an evanescent part. The heat flux radiated by the sample is given by the mean Poynting vector E f ×H f and yields the well known Kirchhoff-Planck radiation law, to which only the propagating modes contribute, whereas the evanescent modes do not figure here, since they are bound to the surface of the sample. Now, if an additional spherical particle with radius R is placed at a distance a above the sample, as depicted in Fig. 1, the fluctuating fields E f and H f induce an electric dipole moment p and a Foucault current j ed within the particle. As emphasized by Chapuis et al. [24], this eddy current j ed is particularly important for metallic particles. It causes losses inside the particle which can be described by means of an effective magnetic dipole moment m eff , in analogy to the electric dipole moment p. We assume a R, so that higher multipoles may be neglected; in principle, this restriction to sufficiently large distances a could be removed by including higher multipoles. Within the dipole approximation, the energy absorbed by the particle can then be written as

P S→P = ṗ • E f + ṁeff • H f ; (1)
the angular brackets denote an ensemble average. For a homogenous and isotropic particle the relations between the induced dipole moments and the fields are given by and

p = ε 0 α E p E f (2)
m eff = µ 0 α H p H f , (3)
where α E p and α H p symbolize the electric and magnetic polarizabilities of the particle; ε 0 and µ 0 are the permittivity and the permeability of the vacuum, respectively. In general, the polarizabilities have a directional dependence and therefore are described by a tensor. For the highly symmetric case of a sphere, this tensor reduces to a scalar multiple of the unit tensor, so that the polarizabilities are represented by scalar values in this case. It is essential to keep in mind that only non-magnetic materials are taken into account here, so that the absorption ṁeff • H f is solely ascribed to the loss due to eddy currents.

With the help of Eqs. ( 2) and (3) one obtains the expression

P S→P (T S ) = ∞ 0 dω ω Im(α E P )ε 0 |E f | 2 ω +Im(α H P )µ 0 |H f | 2 ω (4)
for the energy transferred from the sample to the particle [24,37]. Here |E f | 2 ω denotes the frequency-dependent correlation function of the electric field, i.e., the expectation value of the product of the Fourier transform E f (ω) of the fluctuating electric field and its complex conjugate; the correlation function |H f | 2 ω of the magnetic field is defined analogously.

Reversely, fluctuating charges inside the sphere of temperature T P give rise to energy transfer from the sphere to the flat surface. The power absorbed within the sample, denoted P P→S (T p ), can be calculated in the same manner as P S→P (T S ) above. The resulting expression differs from the expression (4) only through the temperature, so that P P→S (T p ) = P S→P (T P ). As the energy flux P P→S is directed oppositely to P S→P , the resulting overall heat transfer between the particle and the sample has the form [24,37] P S↔P = P S→P (T S ) -P S→P (T P ) .

(5)

The correlation functions 4) are known from fluctuational electrodynamics [START_REF] Rytov | Principles of Statistical Radiation[END_REF]. If one assumes that the sample occupies the infinite half space z ≤ 0, as in Ref. [2], one gets

|E f | 2 ω and |H f | 2 ω in Eq. (
E f,i E * f,i ω = Θ(ω, T S ) 2ω πε 0 c 2 κ<k0 d 2 κ (2π) 2 1 4k z0 e ⊥ ⊗ e ⊥ ii 1 -|r ⊥ | 2 + e pr ⊗ e pr ii 1 -|r | 2 + κ>k 0 d 2 κ (2π) 2 e -2γ a 2γ Im(r ⊥ ) e ⊥ ⊗ e ⊥ ii + Im(r ) e ev ⊗ e ev ii ( 6
)
for the electric correlation function, where

Θ(ω, T ) = ω e ωβ -1 (7)
is the Bose function with the inverse temperature β = 1/(k B T ). Moreover, the unit vectors

e ⊥ := 1 κ (-k y , k x , 0) t e pr := 1 k 0 κ k x k z0 , k y k z0 , κ 2 t e ev := 1 k 0 κ k x γ, k y γ, κ 2 t (8)
have been introduced, writing

κ 2 = k 2 x + k 2 y together with k z0 = k 2 0 -κ 2 and γ = κ 2 -k 2 0 . As usual, k 0 = ω/c
, where c is the velocity of light in vacuum. In addition, we write

r ⊥ = k z0 -k z k z0 + k z r = ε S k z0 -k z ε S k z0 + k z (9)
for the Fresnel reflection coefficients, with k z = εk 2 0 -κ 2 . The first integral in Eq. ( 6), with κ < k 0 , describes the propagating part of the fluctuating field, whereas the second integral with κ > k 0 describes the evanescent part. The expression for H f,i H * f,i ω is obtained from Eq. ( 6) by interchanging r ⊥ ↔ r and 1 ε0 ↔ 1 µ0 , due to a corresponding symmetry of the electric and magnetic Green's functions.

Knowing the correlation functions above a half space, the mean power absorbed by the particle above a planar surface can be calculated from Eq. ( 4), provided the polarizabilities α E P and α H P are given. In the case of a spherical particle of radius R these quantities can be derived from Mie scattering theory [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF]. Denoting the particle's relative permittivity as ε P , and introducing the dimensionless variables x = k 0 R and y = √ ε P k 0 R, one finds [19,[START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF] α E P = 2πR 3 (2ε P + 1) (sin(y) -y cos(y)) -y 2 sin(y) (ε P -1) (sin(y) -y cos(y)) + y 2 sin(y)

α H P = πR 3 3 x 2 -6 y 2 y 2 + 3y cot(y) -3 - 2x 2 5 ( 10
)
for x 1, implying that the particle's radius should be small compared to the dominant thermal wavelength, R λ th . If we demand even |y| 1, meaning that the radius be smaller than the skin depth at thermal frequencies, the above expressions reduce to

α E P = 4π R 3 ε P -1 ε P + 2 , ( 11
)
α H P = 2π 15 R 3 (k 0 R) 2 (ε P -1) . ( 12
)
Evidently, the expression for α E P equals the well known Clausius-Mossotti formula.

Before proceeding, we specify the relevant orders of magnitude. For bulk metals the relative permittivity is well described by the Drude model [START_REF] Ashcroft | Solid State Physics[END_REF]

ε = 1 - ω 2 p ω (ω + iω τ ) . ( 13
)
For gold at room temperature the plasma frequency is ω p = 1.4 × 10 16 s -1 , while the relaxation rate figures as ω τ = 3.3 × 10 13 s -1 . Taking the thermal frequency ω th ≈ 10 14 s -1 one then finds ( | |k 0 ) -1 ≈ 21 nm, setting a rough upper limit to the radius of a gold sphere for which the expressions ( 11) and ( 12) may still hold with reasonable accuracy. When treating such minuscule particles we have to modify the Drude permittivity (13) in order to account for surface scattering, since the bulk mean free path of the electrons, which is about 42 nm for gold, reaches the same order of magnitude as the spheres' diameters. For spherical particles of radius R the required correction is achieved by the replacement [START_REF] Tomchuk | [END_REF]

ω τ → ω τ = 3 4 v F R , ( 14
)
where v F is the Fermi velocity. Taking v F = 1.4 × 10 6 m/s for gold, and setting R = 20 nm, we obtain ω τ = 5.3 × 10 13 s -1 , less than two times the bulk value. On the other hand, size quantization effects become essential only for radii below 2 nm, and may therefore be neglected here.

Dipole model for ellipsoidally shaped particles

The near-field radiative heat transfer between a flat sample with temperature T S and an ellipsoidally shaped particle with temperature T P , as sketched in Fig. 2, can be expressed through equations similar to Eqs. ( 4) and ( 5), but instead of employing the expressions α E/H P for the polarizabilities of a sphere, the proper polarizability tensors for an ellipsoidal particle have to be taken into account.

ε S T S ε P T P R R s p a Figure 2
. Sketch of a spheroidally shaped particle with semiaxes Rs and Rp. The semi-axis with length Rp lies along the spheroid's axis of rotation.

Denoting the three semi-axes as R x , R y , R z , we require a max{R x , R y , R z } ≡ R max to substantiate the dipole approximation, while keeping R max smaller than the skin depth, thus allowing for approximately homogeneous fields inside the ellipsoids. It is textbook knowledge that the induced dipole moment p of an ellipsoidal particle is connected with the incident field E f through the relation [START_REF] Landau | Electrodynamics of Continuous Media[END_REF]

p i = ε 0 V E ε P -1 1 + (ε P -1)n i E f,i (15)
with V E = 4π 3 R x R y R z specifying the volume of the particle, and ε P its permittivity. It is assumed here that the components of p and E f are given in the principal-axis system of the ellipsoid. The polarizability tensor then is diagonal, with diagonal elements given by

α E E,i = V E ε P -1 1 + (ε P -1)n i . ( 16
)
The quantities n i are the so-called depolarization coefficients. They are written as [START_REF] Landau | Electrodynamics of Continuous Media[END_REF]

n i = R x R y R z 2 ∞ 0 ds 1 (s + R 2 i )R(s) (17) with R(s) = (s + R 2 x )(s + R 2 y )(s + R 2 z
). These coefficients depend only on the shape of the ellipsoid, not on its volume. They obey the relations

3 i=1 n i = 1 and n i ≥ 0 . ( 18
)
For the special case of a sphere with R x = R y = R z ≡ R, Eq. ( 17) yields n 1 = n 2 = n 3 = 1/3, so that Eq. ( 16) correctly reduces to the Clausius-Mossotti formula (11).

In the case of spheroid, i.e., of a rotational ellipsoid with two equal semi-axes, R x = R y ≡ R s and R z ≡ R p , the coefficients n i take the form [START_REF] Landau | Electrodynamics of Continuous Media[END_REF] n

1 (e) = n 2 (e) = 1 2 1 -n 3 (e) , n 3 (e) =    1-e 2 2e 3 ln 1+e 1-e -2e , R s < R p 1+e 2 e 3 e -arctan(e) , R s > R p , (19)
with

e 2 =    1 - R 2 s R 2 p , R s < R p R 2 s R 2 p -1 , R s > R p . ( 20
)
For calculating the effective magnetic polarizability of an ellipsoid, we follow a strategy outlined by Tomchuk and Grigorchuk [START_REF] Tomchuk | [END_REF], and use the identity

ω Im(α H E,i )µ 0 |H f,i | 2 ω = Re V E d 3 r j ed • E f,ed * ω , (21)
thus shifting the emphasis from the effective magnetic moment m eff to a fluctuating eddy field E f,ed . This requires that the material is non-magnetic, so that m eff is entirely due to eddy currents j ed . The eddy field obeys the equations

∇ × E f,ed (ω) = iωµ 0 H f (ω) , ( 22
) ∇ • E f,ed (ω) = 0 ; ( 23
)
in principle, the second equation constitutes a boundary condition for E f,ed . Because we assume R max to be smaller than the skin depth, H f can be considered as constant within the volume of the particle, so that the eddy field E f,ed depends linearly on the spatial coordinates. The resulting relation between E f,ed and

H f reads [40] E f,ed,x (ω) = iωµ 0 R 2 x zH f,y (ω) R 2 z + R 2 x - yH f,z (ω) R 2 x + R 2 y . ( 24
)
The other components of the electric eddy field are obtained by cyclic permutation of the indices. We emphasize that the underlying assumption of constant fields H f within the particle is not well fulfilled in cases where the particle's characteristic linear dimensions are on the order of the skin depth, which is about 21 nm for gold at room temperature. In such cases, the formalism outlined here may still yield the correct orders of magnitude, but not exact numbers.

Next, one has the identity [START_REF] Tomchuk | [END_REF]

j ed (ω) = ωε 0 ε P E f,ed (ω) , ( 25
)
which relates the induced eddy currents to E f,ed , with ε P denoting the imaginary part of the particle's permittivity. Thus, for a spheroid the absorbed energy ( 21) is found to be Re

V E d 3 r j ed • E f,ed * ω = ω 3 c 2 µ 0 ε P V E 10 × |H f,z | 2 ω R 2 s + |H f,x | 2 + |H f,y | 2 ω 2R 2 s R 2 p R 2 s + R 2 p . ( 26
)
The effective polarizability is finally read off by comparing this expression with the second term in Eq. ( 4). Again the polarizability tensor has non-zero entries only on its diagonal, with imaginary parts given by

Im(α H E,x ) = Im(α H E,y ) = (k 0 R s ) 2 R 2 p R 2 s + R 2 p V E ε P 5 , Im(α H E,z ) = (k 0 R s ) 2 10 V E ε P . ( 27
)
For a sphere with R s = R p = R, both expressions give the correct imaginary part of the Mie formula (12).

Im(α E E ) / Im(α P E ) R s / R p Im(α E,x E ) + Im(α E,z E ) Im(α E,x E ) Im(α E,z E ) Figure 3. Imaginary parts Im(α E E,x 5
) and Im(α E E,z ) of the electric polarizabilities (16) at ω th = 10 14 s -1 for a spheroidal gold particle, and their sum, vs. the ratio R s /R p . The spheroid's volume V E is kept constant at that of a sphere with radius R = 20 nm; surface scattering is taken into account. Data are normalized by the imaginary part of the electric polarizability Im(α E P ) of the sphere.

In Figs. 3 and4 we plot the imaginary parts of the electric and of the magnetic polarizability ( 16) and ( 27) at the thermal frequency ω th = 10 14 s -1 against the ratio R s /R p , keeping the volume V E constant at that of a sphere with radius R = 20 nm. The Drude expression (13) has been taken for the permittivity, with the plasma frequency ω p = 1.4 × 10 16 s -1 for gold at room temperature, while surface scattering has been taken into account through the replacement (14), inserting the geometric mean (R 2 s R p ) 1/3 for R. This replacement is not fully correct for nonspherical particles [START_REF] Tomchuk | [END_REF]: In principle, one then has different scattering times for different directions, leading to anisotropy ) and Im(α H E,z ) of the magnetic polarizabilities (27) at ω th = 10 14 s -1 for a spheroidal gold particle, and their sum, vs. the ratio R s /R p . The spheroid's volume V E is kept constant at that of a sphere with radius R = 20 nm; surface scattering is taken into account. Data are normalized by the imaginary part of the magnetic polarizability Im(α H P ) of the sphere.

of the permittivity. In order to estimate the resulting error, we compare in Fig. 5 the electric absorptivities provided by the Drude permittivity with the simple surface correction ( 14), and without any correction at all. While the effect of the correction is clearly visible, the general trends are not altered. Hence, a more refined correction taking anisotropy into account should not give drastically different results, at least not in the interval 1/5 ≤ R s /R p ≤ 5 considered. The curve progressions observed here can be understood intuitively. For the electric polarizability depicted in Fig. 3 the absorptivity Im(α E E,z ) decreases monotoni-cally with increasing ratio R s /R p , whereas the absorptivity Im(α E E,x ) increases. For R s /R p 1 the spheroidal particle becomes needle-like along the z-direction; hence, the polarizability in z-direction is much greater than that in x-or y-direction. On the other hand, when R s /R p 1 the spheroid is pancake-like and lies parallel to the x-y-plane; therefore the polarizability in the directions perpendicular to the z-direction dominates.

In the case of the magnetic polarizability shown in Fig. 4 one observes that both Im(α H E,x ) and Im(α H E,z ) go to zero in the formal limit R s /R p → 0. This is clear because these two quantities represent absorptivities caused by eddy currents perpendicular to the x-and to the zdirection, respectively. Accordingly, for thin needle-like spheroids with R s /R p 1 both absorptivities disappear, because no eddy current can arise then. In the opposite case, for R s /R p 1, eddy currents can be easily induced perpendicular to the z-direction in the pancake-like particle, whereas perpendicular to the x-direction no eddy currents occur, as before. Hence, for

R s /R p → ∞ one ex- pects Im(α H E,z ) → ∞ and Im(α H E,x ) → 0. Since Im(α H E,x
) is positive for all values of R s /R p and goes to zero in both limiting cases there has to be a maximum, which occurs near R s /R p = 1, as witnessed by Fig. 4.

Radiative heat transfer between an ellipsoidal nanoparticle and a planar surface

In accordance with Sec. 2, the radiative heat transfer from a planar sample to an ellipsoidal nanoparticle is expressed, in analogy to Eq. ( 4), as

P S→E (T S ) = ∞ 0 dω ω Im(α E E,i )ε 0 |E f,i | 2 ω +Im(α H E,i )µ 0 |H f,i | 2 ω , ( 28
)
where summation over repeated indices is implied. It is assumed here that the particle is oriented as in Fig. 2, so that the surface normal is parallel to a principal axis of the ellipsoid. As discussed before, this form ( 28) is quite general, requiring only the specification of the polarizabilities of the nanoparticle, and of the correlation functions of the fluctuating fields above the sample's surface. The expression P E→S for the energy flux from the nanoparticle back to the sample can again be obtained from Eq. ( 28) by substituting the temperature T P for T S .

When the specific expressions for a spheroid above a planar sample are inserted into Eq. ( 28), the propagating modes with κ < k 0 yield the contribution

P pr,S→E = 2 π ∞ 0 dω Θ(ω, T S )k 2 0 κ<k0 d 2 κ (2π) 2 1 4k z0 Im(α E E,x ) (1 -|r ⊥ | 2 ) + k 2 z0 k 2 0 (1 -|r | 2 ) +Im(α E E,z ) κ 2 k 2 0 (1 -|r | 2 ) +Im(α H E,z ) (1 -|r | 2 ) + k 2 z0 k 2 0 (1 -|r ⊥ | 2 ) +Im(α H E,x ) κ 2 k 2 0 (1 -|r ⊥ | 2 ) . ( 29
)
This reduces correctly to the known result for a sphere by setting R s = R p = R:

P pr,S→P = 2 π ∞ 0 dω Θ(ω, T S )k 2 0 κ<k 0 d 2 κ (2π) 2 1 4k z0 × Im(α E P + α H P ) (1 -|r | 2 ) + (1 -|r ⊥ | 2 ) . ( 30
)
Observe that sum of the polarizabilities figures here, because

ε 0 |E f | 2 ω = µ 0 |H f | 2 ω
for the propagating modes. Accordingly, the relative size of Im(α E P ) and Im(α H P ) determines whether the radiative heat transfer is dominated by the electric or by the magnetic part. From Eqs. ( 29) and ( 30) one can calculate the heat radiated by an ellipsoidal or a spherical nanoparticle in the absence of the sample [START_REF] Yu | [END_REF] by simply setting |r ⊥ | = |r | = 0 and inserting the temperature T P instead of T S .

The heat transfer mediated by the evanescent modes with κ > k 0 from the sample to the spheroidal particle takes the form

P ev,S→E = 2 π ∞ 0 dω Θ(ω, T S ) k 2 0 κ>k 0 d 2 κ (2π) 2 e -2γa 2γ Im(α E E,x ) Im(r ⊥ ) + γ 2 k 2 0 Im(r ) + Im(α E E,z ) κ 2 k 2 0 Im(r ) + Im(α H E,x ) Im(r ) + γ 2 k 2 0 Im(r ⊥ ) + Im(α H E,z ) κ 2 k 2 0 Im(r ⊥ ) . ( 31
)
Again this result leads directly to the corresponding expression for a sphere [24]:

P ev,S→P = ∞ 0 dω Θ(ω, T S ) k 2 0 κ>k 0 d 2 κ (2π) 2 e -2γa 2γ Im(α E P ) Im(r ⊥ ) + 2κ 2 -k 2 0 k 2 0 Im(r ) +Im(α H P ) Im(r ) + 2κ 2 -k 2 0 k 2 0 Im(r ⊥ ) . ( 32
)
Note that in general one has

ε 0 |E f | 2 ω = µ 0 |H f | 2 ω
for the evanescent modes. In contrast to the transfer by propagating modes, this implies that the magnetic contribution to the heat transfer, which is proportional to Im(α H P ), can dominate the electric one [24] even if Im(α E P ) > Im(α H P ). In the quasi-static regime, i.e., for distances a even smaller than the substrate skin depth, one has κ | |k 0 , so that Eq. ( 31) leads to the approximation

P ev,S→E ≈ 2 π ∞ 0 dω Θ(ω, T S ) k 2 0 ∞ 0 dκ 2π κ 2 e -2κa k 2 0 × Im(α E E,x + α E E,z )Im(r ) + Im(α H E,x + α H E,z )Im(r ⊥ ) . ( 33
)
Furthermore, using the approximate expressions

Im(r ) ≈ 2ε S |ε S + 1| 2 , ( 34
)
Im(r ⊥ ) ≈ ε S k 2 0 4κ 2 (35)
for the imaginary parts of the reflection coefficients, the integration over the wave number κ can be carried out, giving

P ev,S→E ≈ 2 π 2 ∞ 0 dω Θ(ω, T S ) × Im(α E E,x + α E E,z ) ε S |ε S + 1| 2 1 (2a) 3 + Im(α H E,x + α H E,z ) k 2 0 ε S 32a . ( 36
)
As expected, in the quasi-static regime one has P E ev,S→E ∝ a -3 for the electric and P H ev,S→E ∝ a -1 for the magnetic contribution; these power laws are the same as those for a sphere. Since the sums Im(α E E,x +α E E,z ) and Im(α H E,x +α H E,z ) appear in the frequency integral, one may expect that in this regime P ev,S→E behaves like the imaginary part of the sum of the polarizabilities at the thermal frequency when R s and R p are varied.

It needs to be stressed, however, that this approximate result (36) has to be taken with some caution, since the required short particle-sample distances may conflict with those required by the dipole approximation.

Discussion

Now we study the shape dependence of the radiative heat transfer between a spheroidal nanoparticle and a planar sample numerically, on the basis of the full expressions (29) and (31). To this end, we assume the temperatures T P = 100 K for the particle and T S = 300 K for the sample. The permittivity of the sample is described by the bulk Drude model ( 13) with parameters for gold; that of the particle by the Drude expresson with the plasma frequency of gold, and with the relaxation rate adapted to surface scattering, as before. Since we are interested in the variation of the heat transfer caused by alterations of the shape of the nanoparticle, we vary the ratio R s /R p in such a way that the volume V E of the respective spheroidal nanoparticle equals the volume V P of a sphere with radius R = 20 nm. We estimate that the range 1/5 ≤ R s /R p ≤ 5 is about the largest that is still compatible with the dipole approximation, and restrict ourselves to this interval. E to the radiative heat transfer between a spheroidal gold particle with temperature T P = 100 K and a planar gold surface with temperature T S = 300 K, normalized to the corresponding value P E P for a sphere. When the ratio R s /R p is varied, the particle's volume is kept constant at that of a sphere with radius R = 20 nm. Line styles distinguish data for particle-sample distances ranging from 50 nm to 2000 nm.

In Fig. 6 we depict the electric contribution to the total heat transfer P E E vs. R s /R p for several distances a between 50 nm and 2000 nm, normalized to the corresponding energy transfer P E P between the isochoric sphere and the sample. Evidently, for small distances the curves of P E E /P E P actually resemble those of Im(α E E,x + α E E,z ) in Fig. 3, as may have been conjectured from the quasistatic approximation (36), notwithstanding its somewhat shaky justification. On the other hand, for relatively large distances the progression of P E E /P E P is similar to that of Im(α E E,z ) alone. This behavior stems from the fact that in the latter regime, i.e., at distances such that the propagating modes dominate the heat transfer, one has |E z,pr | 2 |E x,pr | 2 + |E y,pr | 2 , as documented in Fig. 7. In contrast, in the near-field regime where the evanescent modes dominate the transfer one finds |E z,ev | 2 = |E x,ev | 2 +|E y,ev | 2 , so that Im(α E E,x ) and Im(α E E,z ) contribute in equal measure there.

Analogously, we plot in Fig. 8 E to the radiative heat transfer between a spheroidal gold particle with temperature TP = 100 K and a planar gold surface with temperature TS = 300 K, normalized to the corresponding value P H P for a sphere. When the ratio R s /R p is varied, the particle's volume is kept constant at that of a sphere with radius R = 20 nm. Line styles distinguish data for particle-sample distances ranging from 50 nm to 5000 nm.

The ratio of the magnetic to the electric contribution P H E /P E E is drawn in Fig. 9, again for distances from 50 nm to 5000 nm. For all distances considered the magnetic contribution to the heat transfer is much greater than the electric one, unless R s is much smaller than R p . This fact has been discussed in detail by Chapuis et al. [24] for a spherical metallic nanoparticle above a planar metallic sample. Here we find that the curves of the ratio P H E /P E E exhibit a maximum near R s /R p = 1, so that the dominance of the magnetic part is somewhat less pronounced for nanoparticles with shapes differing from a sphere. For markedly needle-like, prolate particles with a very small ratio R s /R p , the graphs suggest that the electric contribution could even dominate the magnetic one, because the induction of Foucault currents would be suppressed by the needle-like shape. We remark that for non-metallic nanoparticles and non-metallic samples the contribution due to induced eddy currents can be neglected anyway, so that for such materials it is always the electric contribution which dominates the heat transfer. We conclude from Figs. 6 and 8 that the near-field heat transfer between a metallic nanoparticle and a sample depends to a sizeable extent on the nanoparticle's shape even when the total radiating volume is kept constant. In the example considered, the electric contribution to the heat transfer is enhanced for strongly prolate spheroids by a factor of about 10 as compared to a perfect sphere, and by a factor of about 4 for strongly oblate ones. On the other hand, the magnetic contribution is substantially reduced in the needle-like case, but exceeds the sphere value by a factor of about 2 for pancake-like particles, always assuming the orientation specified in Fig. 2.

Conclusions

We have investigated the radiative heat transfer between a metallic spheroidal nanoparticle and a planar metallic probe for small and large distances, on the basis of the analytical expressions (29) and (31) for the heat transfer mediated by propagating and evanescent modes, respectively. By numerical analysis for the example of a spheroidal gold particle above a gold surface we have investigated in detail the precise shape dependence of the radiative heat transfer in both the near-field and the far-field regime. Furthermore, we have derived the approximate expression (36) for the nonretarded evanescent regime, which provides a qualitative understanding of the numerical results. . Heat transfer between a gold spheroid with temperature T P = 100 K and a planar gold surface with temperature T S = 300 K separated by 100 nm. The volume of the spheroid is kept constant at that of a sphere with radius R = 20 nm when R s /R p is varied. The horizontal line marks the value of the heat current between the isochoric sphere and the sample.

Figure 10 shows the absolute total heat transfer between a gold spheroid with temperature T P = 100 K and a planar gold surface with temperature T S = 300 K for the separation a = 100 nm, again fixing the spheroid's volume at that of a sphere with radius R = 20 nm when changing the aspect ratio. For R s /R p = 1/5 the heat transfer between spheroid and sample is only about half of that between sphere and sample, whereas for R s /R p = 5 the spheroid-sample transfer is roughly two times as efficient as that occurring in the sphere-sample geometry. The very fact that there exists such a marked shape dependence might be of interest for nanoscale thermal engineering, insofar as it appears possible to control the amount of heat transported at nanoscale distances by carefully designing the shapes of both the emitting and the absorbing pieces.

Our discussion is subject to several restrictions: There is the dipole approximation (1) made right at the outset, the assumption of constant fields H f inside the particle entering into Eq. ( 24), and the simplified correction (14) for surface scattering. Taken together, these simplifications render our analysis approximate, rather than exact, although they should still capture the essential physics. Further issues that should be investigated in future works concern the possible effects of spatial dispersion [43], and of surface roughness [START_REF] Biehs | Near-field heat transfer between a nanoparticle and a rough surface[END_REF].

With respect to the problem of quantifying the actual heat transfer in a near-field scanning thermal microscope, our study helps to pin down the error margin. A typical NSThM sensor tip is larger than the nanoparticles considered in our examples, and does possess an internal structure, but again the precise shape of the sensor will influence the heat current it records. In view of our model calculations, we estimate that the values obtained on the grounds of the dipole model under the assumption of an ideal spherical sensor may deviate, in the appropriate distance regime, by up to an order of magnitude from the true values, but probably not by more.

Figure 1 .

 1 Figure 1. Sketch of a spherical nanoparticle with radius R, temperature TP, and permittivity εP placed at a distance a above a planar sample with temperature T S and permittivity ε S .

Figure 4 .

 4 Figure 4. Imaginary parts Im(α H E,x) and Im(α H E,z ) of the magnetic polarizabilities(27) at ω th = 10 14 s -1 for a spheroidal gold particle, and their sum, vs. the ratio R s /R p . The spheroid's volume V E is kept constant at that of a sphere with radius R = 20 nm; surface scattering is taken into account. Data are normalized by the imaginary part of the magnetic polarizability Im(α H P ) of the sphere.

Figure 5 .

 5 Figure 5. Electric absorptivities for spheroidal gold particles with a fixed volume corresponding to that of a sphere with radius R = 20 nm, both with (corrected) and without (bulk data) surface scattering taken into account.

Figure 6 .

 6 Figure 6. Electric contribution P E E to the radiative heat transfer between a spheroidal gold particle with temperature T P = 100 K and a planar gold surface with temperature T S = 300 K, normalized to the corresponding value P E P for a sphere. When the ratio R s /R p is varied, the particle's volume is kept constant at that of a sphere with radius R = 20 nm. Line styles distinguish data for particle-sample distances ranging from 50 nm to 2000 nm.

  Figure 10. Heat transfer between a gold spheroid with temperature T P = 100 K and a planar gold surface with temperature T S = 300 K separated by 100 nm. The volume of the spheroid is kept constant at that of a sphere with radius R = 20 nm when R s /R p is varied. The horizontal line marks the value of the heat current between the isochoric sphere and the sample.

  Ratios |Ex| 2 +|Ey| 2 / |Ez| 2 (full line) and |Hx| 2 + |Hy| 2 / |Hz| 2 (dashed line) above a semi-infinite planar gold sample at T S = 300 K, as function of the distance a.shown in Fig.4, since |H z,ev | 2 = |H x,ev | 2 + |H y,ev | 2for small distances (see Fig.7). On the other hand, for relatively large distances a the plots of P H E /P H P are quite similar to that of Im(α H E,x ) alone, since then |H z,pr | 2 |H x,pr | 2 + |H y,pr | 2 .

		2
			electric
	+ |E y | 2 >/<|E z | 2 >)	1	magnetic				+ |H y | 2 >/<|H z | 2 >)
	log(<|E x | 2	0					log(<|H x | 2
		-1
			-8	-7	-6	-5
				log(a / m)
	1 2 3 Figure 7. 0 P E H 1000 nm 50 nm 2000 nm 3000 nm / P P H 5000 nm
		0	1	2	3	4	5
					R s / R p
	Figure 8. Magnetic contribution P H
							the normalized mag-P to the total heat transfer, for E /P H netic contribution P H distances a extending from 50 nm up to 5000 nm. Similar
							to the electric case, in the near-field regime the graphs of P H E /P H P resemble that of Im(α H E,x + α H E,z ) previously

7 Acknowledgments This work was supported by the Deutsche Forschungsgemeinschaft through Grant No. KI 438/8-1. S.-A. B. gratefully acknowledges a fellowship awarded by the Deutsche Akademie der Naturforscher Leopoldina under Grant No. LPDS 2009-7.


B Abdel

M F Blanc Mignon

M Roumie

A Siblini

J P Chatelon

M Korek

Physico-Chemical Characterization of Multilayer YIG Thin Film Deposited by RF Sputtering

Yttrium iron garnet (YIG) film made as a magneto-optical medium suffers from the problem of crack formation, caused by the heating process. YIG thin film is deposited by radio frequency RF magnetron sputtering; the obtained layer is amorphous and it needs annealing to be crystallized. After heat-treatment at 740 o C of the sample realized on quartz substrate, we observe cracks on the entire film surface. This is due to the large difference between the thermal expansion coefficient (5.5 x 10 -7 K -1 for quartz and 10 x 10 - 6 K -1 for YIG). In this paper we present a new fabrication method to reduce this problem, we make a multilayer to obtain at the end a uniformly unique layer with excellent crystalline structure. Such films have the possibility to reach a thickness of 500 nm. YIG films have been studied by Rutherford Backscattering Spectrometry (RBS), optic Ellipsometry and the Scan Electron Microscope. The RBS spectra were collected in channelling geometry with incident particles energy 2 MeV and 3.5 MeV. The thickness and the stoichiometric value of the thin films have been evaluated. Simulation of all spectra indicates a constant composition. Ellipsometry method is well adapted to model the thin film structure layers, and to measure the thickness of the film and the complex index of refraction. The theoretical ellipsometric value of the index of refraction is (2.22) while the experimental value is ranging from 2.2 to 2.3 for a wavelength of 1550nm.

1 Introduction Some recent works on preparation of magnetic thin films are useful for magnetooptical applications [1,2]. The YIG film made the objective of a great number of studies for reasons related to its magnetic and magneto-optical properties in the fields of optics and microwave frequencies. The YIG is a ferrimagnetic garnet crystal with the composition of Y 3 Fe 5 O 12 , with a substantial Faraday rotation in large parts of the optical and microwave spectrum [3]. YIG thin films are grown using various techniques: the conventional technology which is liquid phase epitaxy (LPE) [5], pulsed laser deposition (PLD) [6] and radio frequency (RF) magnetron sputtering [2][3][4]. In this work, we have used the RF-magnetron sputtering for its several advantages: dry process during deposition, possibility of high purity starting material, ability to sputter dielectric materials, commonly held industrial process and high compatibility with semiconductors technology. For this purpose, the radio frequency sputtering grown films can be used as high quality optical waveguides exhibiting the Faraday effect. Since the quartz substrate is cheaper then the GGG substrate, so we are interested to make a thin film of YIG on quartz substrate. The difficulty after heat-treatment at 740 o C to crystallize the film, is the cracks on the entire surface that are observed [7,8]. The stress during heat treatment is due to the large difference between thermal expansion coefficient of quartz and the film. For this reason we have decided to make a multilayer YIG thin film on a quartz substrate for optical and magneto-optical applications. We studied the deposition conditions on the physicochemical, magnetic and magneto-optical properties of the thin film of YIG, in order to determine the optimal parameters of deposition for obtaining a good quality film. The film must have some characteristics like smooth surface, correct stoechiometric values and suitable refractive index as necessary for magneto-optical applications. The essential aim of the new method of the present work is to carry out films of YIG on quartz substrate having a good surface quality with the minimum of cracks after a heat treatment.

Experimental

Film preparation

In order to optimize the surface of the thin films deposited by RF magnetron sputtering, the following conditions should be set up such as argon pressure and flow rate, RF power and the sputtering time. YIG thin film was deposited by using RF magnetron non reactive sputtering system onto a quartz substrate [3]. The working pressure is 2x10 -2 mbar and an argon gas of 99.999% is introduced into the sputtering chamber to create the plasma. The magnetron input power is 100 W. The rate of YIG film deposition is 0.3 μm/h. A thermal treatment promoted crystallization of the film; the samples are heated in vacuum. It is best not to exceed 300 degrees per hour in order to avoid thermal stress cracks in the thin film starting from room temperature. The film is heated at a rate of 4 o C/min till reaching 740 o C, and is maintained at this temperature for 2 hours.

We repeat this process many times (deposition and thermal treatment) until we obtained the required film thickness.

RBS Measurements A 5SDH pelletron tandem accelerator of 1.7 MV was used to perform

Rutherford Backscattering Spectrometry measurements on the samples under normal incident beam and in a random direction to avoid channeling. The RBS is based on collisions between atomic nuclei. By this technique, it is possible to determine the atomic mass and elemental concentrations versus depth below the surface. YIG thin film is bombarded with a beam of high energy particles ( 4 He + of energy 2 MeV) [9]. The scattering angle θ is 165° referring to the beam direction, where the solid angle of the detector is 5.45 × 10 -3 sr. The energy measured for the particle backscattering depends on the mass of the target atoms. RBS measurements was done and showed the influence of the multilayer film growth on the surface state of the film.

Ellipsometry Measurements

The characterization of the optical constants and thickness of YIG thin films is a major part of our research, and ellipsometry is the primary method to determine these quantities. The instrument used is a variable angle spectroscopic and it is sensitive to several material characteristics, such as: layer thickness, optical constants (refractive index). The spectroscopic ellipsometer features high accuracy determination of the ellipsometric angles across their wide spectral range allowing superior precision and sensitivity for the characterization of transparent substrates thin films with low index contrast. For the purposes of our research we are most interested in determining film thickness and optical constants. The concentration of the effects of the Fresnel interference is essential in thin films design [START_REF] French | Proceedings Of SPIE[END_REF] since it can lead to large changes in the reflection and transmission of both thin film and substrate. The obtained data are fitted according to optical model permit us to determine the complex index of refraction and the film thickness.

SEM Measurements

Scanning electron microscope (SEM) was performed with a JEM 100 C instrument equipped with a high resolution scanning device 9ASID-4D0 and a side-entry goniometer. An accelerating voltage of 40 KV was used images were recorded at a Kodak PXP-120 film. The resolution of the taken images is magnified from 100 times to 3000 times.

3 Results YIG thin films are deposited on quartz substrate by using a RF sputtering system. The control of the surface condition and the thickness of the films were measured using the (SEM) and a mechanical profilometer (table1), where the final thickness is 500 nm. The film crystallizes at a temperature 740 o C, correct X-ray diffraction patterns are observed as shown in figure 1; the diffraction peaks observed being indexed to the JCPDS card 43-0507 and they corresponds to YIG.

Table 1: the layer with its corresponding thickness

Number of layers Final thickness ( in nm )

First layer 150

Second layer 365

Third layer 500 Only at high magnification resolution (3000x), a cracks appear on the surface in figure 2.c. The thickness seems to be a limit, because crakes were observed for thickness higher than 500 nm, so the stress on the surface of the film can be reduced by the film thickness. The chemical stoichiometry of the film was analyzed by (RBS). The spectra recorded with the incident beam of 2 MeV energy of 4 He + and backscattering ions were analyzed for an angle θ = 165 o . Figure 3 represents the fitting of the RBS spectra, the four sharp discontinuities observed correspond to 4 He + ions backscattered from O, Si, Fe and Y atoms. exceeds the prediction by about 0.178. For the O atoms, the excess is somewhat higher and it is 0.14.

Such differences between experimental results and predictions for Fe and O are not large. In the case of Y atoms the difference is larger. They may be due to the dependence of the mass of the atoms and the thickness of the film.

Figure 4 represents the ellipsometry fitting data to determine the complex refractive index. A parameterized model was used to describe the YIG optical constants over the wide spectral range:

from the far UV (190 nm) to near-IR (2100 nm). The ellipsometry determine the refractive index (n) and extinction coefficient (k) with a high degree of certainty for YIG thin film. The ellipsometry results show that the complex index of refraction is between 2.2 and 2.3 at a wavelength 1550 nm (the theoretical value is 2.2) and it is independent on the thickness of the film.

Fig 1 :Fig 2 :

 12 Fig 1: X-ray diffraction pattern of YIG multilayer and monolayer films.

Figure 2 .

 2 Figure 2.a presents the SEM photograph of YIG film on quartz substrate shows a formation of cracks on the surface due to the difference of thermal expansion between the YIG film and the substrate. In order to reduce the crake's formation a thin film is made by a multilayer annealed in each layer at 740 o C shown in figure 2.b. No crakes appear on the surface of the film for the same SEM resolution. The comparison between the SEM photos of the global method with those of the new method "multilayer" for the same magnification (500x) does not reveal the presence of cracks.

  present at the surface and in the film; also it is possible to simulate the RBS spectra using the SIMNRA program to obtain the stoechiometic value. Simulations and fits of the spectra are done where the simulated spectrum overlies the experimental data. The obtained stoichiometry of YIG thin film annealed at 740 o C is Y 2.68 Fe 5.178 O 12.14 . The difference for Y atoms of the experimental value with respect to the theory is only 0.32. The experimental value for Fe atoms

Fig 3 :Fig 4 :

 34 Fig 3: Fitting of RBS spectra for YIG thin film.


Jan Siem

Koopman
email: s.j.koopman@feweb.vu.nl

Neil Shephard
email: neil.shephard@economics.ox.ac.uk

Drew Creal
email: dcreal@feweb.vu.nl

Testing the assumptions behind importance sampling

Keywords: Extreme value theory, Importance sampling, Simulation, Stochastic volatility. JEL classification codes: C10, C22, C32, C51

Importance sampling is used in many areas of modern econometrics to approximate unsolvable integrals. Its reliable use requires the sampler to possess a variance, for this guarantees a square root speed of convergence and asymptotic normality of the estimator of the integral. However, this assumption is seldom checked. In this paper we use extreme value theory to empirically assess the appropriateness of this assumption. Our main application is the stochastic volatility model, where importance sampling is commonly used for maximum likelihood estimation of the parameters of the model.

Introduction

This paper develops tools for evaluating the importance sampling techniques that are popular in econometrics. There are two situations where importance sampling is commonly applied. It is often used to numerically evaluate the likelihood function in missing data models, when an analytic expression for the likelihood function does not exist. This is of primary importance in the literature on maximum likelihood estimation of stochastic volatility (SV) models. Importance sampling is also popular in Bayesian applications when computing expectations of the posterior distribution.

In a fundamental contribution, [START_REF] Geweke | Bayesian inference in econometric models using Monte Carlo integration[END_REF] argued that importance sampling should only be used in settings where the variance of the importance weights is known to exist. Failure of this condition can lead to slow and unstable convergence of the estimator as the central limit theorem governing convergence fails to hold. Robert and Casella (2004, p.90) provide detailed examples of importance samplers that fail this condition and note that ignoring the problem can result in strongly biased estimators. While the variance conditions can be checked analytically in low dimensional problems, proving that they are met in high dimensions can be challenging.

Many econometrics and statistics papers have recently been written which, in effect, a priori assume that this condition holds. This is the case in the stochastic volatility literature; see, e.g., [START_REF] Hendry | Likelihood evaluation for dynamic latent variable models[END_REF], [START_REF] Danielsson | Accelerated Gaussian importance sampler with application to dynamic latent variable models[END_REF], [START_REF] Sandmann | Estimation of stochastic volatility models via Monte Carlo maximum likelihood[END_REF], [START_REF] Elerian | Likelihood inference for discretely observed non-linear diffusions[END_REF], [START_REF] Durham | Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes (with discussion)[END_REF], [START_REF] Durham | Monte Carlo methods for estimating, smoothing, and filtering one-and two-factor stochastic volatility models[END_REF], and [START_REF] Asai | Asymmetric multivariate stochastic volatility[END_REF]. It is also assumed to hold in the macroeconomics literature on dynamic stochastic general equilibrium models; see, e.g. [START_REF] Creal | Bayesian analysis of DSGE models[END_REF].

In this paper, we develop methods to check for the existence of the variance of the importance weights. This will be based on the application of extreme value theory. 1 Limit theorems from extreme value theory imply that we can learn about the variance of the weights by studying the behavior of their distribution in the right hand tail. Earlier work on testing whether the variance of the importance weights exists also using extreme value theory was conducted by [START_REF] Monahan | Testing the behaviour of importance sampling weights[END_REF][START_REF] Monahan | Numerical Methods of Statistics[END_REF]. Monahan examined the tail-behavior of the weights by building a test statistic from the [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF] estimator. We take an alternative approach and develop three different test statistics using the well-known principles of the Wald, score and likelihood ratio tests. We also propose graphical diagnostics for a complete insight.

The new test procedures are carried out for two Monte Carlo studies and for the maximum likelihood analysis of two SV models. The latter models have been our main motivation to investigate the problem in detail. It will be shown that sometimes the variance of the importance weights does not exist when estimating parameters of SV models. Although our main focus will be on likelihood estimation, our methods apply generally to all applications of importance sampling including Bayesian approaches.

A review of importance sampling is provided in Section 2. In Section 3, we discuss the statistical estimation of tail indexes for extremes and introduce our test statistics. Section 4 includes results from some Monte Carlo experiments that are conducted to evaluate the power and size of the proposed tests. We apply the test statistics and graphical diagnostics to maximum likelihood estimation of the parameters of two SV models in Section 5. Section 6 concludes.

1 We should note that extreme value theory has been applied in the context of financial economics and insurance in order to determine the thickness of the tail of financial returns. Furthermore, it is used to compute various measures of risk. References to this literature include Embrechts, Kluppelberg, and Mikosch (1997, Section 6.5.2) and [START_REF] Danielsson | Tail index and quantile estimation with very high frequency data[END_REF]. The application of the theory is much harder in that context for the i.i.d. assumption does not hold.

Importance sampling was discussed as early as [START_REF] Kahn | Methods of reducing sample size in Monte Carlo computations[END_REF] and [START_REF] Marshall | The use of multi-stage sampling schemes in Monte Carlo computations[END_REF], while it was popularised in the influential monograph by Hammersley and Handscomb (1964, Section 5.4). It was first used in econometrics by [START_REF] Kloek | Bayesian estimates of equation system parameters: an application of integration by Monte Carlo[END_REF] in their work on computing posterior densities. Consider an integral that cannot be solved analytically

I = X h(x)f (x)dx, (1)
where h(x) is a function and f (x) is a probability density. An importance sampling density g(x)

is selected from which it is easy to simulate and whose support is also X . The integral I can be approximated by

I = 1 N N i=1 w i , where w i = w(x i ),
for a large number of N draws and with

w(x) = h(x)f (x) g(x)
and x i ∼g(x).

(2)

The strong law of large numbers implies that I will be a consistent estimator of I without any additional assumptions. By construction we know that the {w i } N i=1 are i.i.d. and that E[w(x)] = I. As a result, a simple application of Kolmogorov's strong law of large numbers (e.g. [START_REF] Geweke | Bayesian inference in econometric models using Monte Carlo integration[END_REF] and Geweke (2005, p. 114)) shows that I a.s.

→ I,

as N → ∞, whatever importance sampler we design. In practice, the normalizing constant of f (x) or g (x)

in (2) may not be known requiring a self-normalized importance sampling estimator.2

In order to easily measure the precision of I and to guarantee that the rate of convergence to

I is N 1/2
, it is helpful to have a Gaussian central limit theorem for I. We know that a necessary condition for the application of the Lindeberg-Lévy central limit theorem is the existence of a variance for the importance weights. This would allow us to conclude that

√ N I -I d → N (0, γ 2 ).
Expressions for the asymptotic variance γ 2 of the importance sampling estimator can be found in Geweke (2005, p. 114) and Cappé, Moulines, and Rydén (2005, p. 287). However, the existence of the variance of the importance weights is by no means guaranteed. If it does not exist, the importance sampling estimator may convergence erratically making estimates unreliable. An approach for testing for the existence of the variance is considered in the next section.

3 Extreme value theory and inference

Extreme value theory

We begin this section briefly reviewing the results from extreme value theory that are relevant for this paper; see [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF] for a detailed discussion. Assume we have a sequence w 1 , w 2 , ..., w N of i.i.d random variables from a distribution function F . The {w i } N i=1 need not represent importance weights, although that will be the case in our work below. Denote the endpoint of the distribution function F as w E = sup {w : F (w) < 1} ≤ +∞. [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF] defined a new sequence of n random variables Z 1 , ..., Z n created from those w i that lie above a threshold u, i.e. Z i = w i -u where n ≤ N . [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF] proved that as N and u get larger the distribution of the excesses {Z i } n i=1 is well approximated by the generalized Pareto distribution. More formally, he showed that lim

u→w E sup 0<z<w E -u |F u (z) -GP D (z; ξ, β)| = 0,
where F u (z) is the distribution function of the excesses conditional on threshold u and w > u.

GP D denotes the generalized Pareto distribution with density

f (z; ξ, β) = 1 β 1 + ξ z β -1 ξ -1 , z ∈ D(ξ, β) > u, β > 0. (3)
Here

D(ξ, β) =    [0, ∞), ξ ≥ 0, [0, -β/ξ], ξ < 0.
The threshold u plays a crucial role in this result and its selection will be discussed in Section 4. We continue by assuming that u is a given threshold.

The shape parameter ξ of the generalized Pareto distribution characterizes the thickness of the distribution's tails. In the context of extreme value theory, the value of ξ is determined by the thickness of the tails of F . This means that we can extract information about the thickness of the tails of the importance weights by making inference on the shape parameter using the excesses {Z i } n i=1 . In fact, it can be shown that when ξ > 0 then E[w j ] = ∞ for j ≥ 1/ξ. Consequently, the number of moments that the distribution of the importance weights has can be determined by focusing on ξ. In particular, we are interested in establishing that ξ ≤ 1/2. We know by construction of the importance sampler that ξ ≤ 1. The cases of ξ < 0 deal with situations where the {w i } N i=1 have some upper bound. This is of some relevance in importance sampling in the case of the sampler being bounded, that is

f (x) ≤ M g(x), (4)

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

for some finite choice of M > 0. Then if (4) is true for x ∈ X , the support of the importance weights {w i } N i=1 will be bounded and the existence of all moments is guaranteed. This is, in practice, quite unusual. In most cases of real interest, the sampler is not bounded, so ξ ∈ [0, 1].

Previous work: Monahan's test

Monahan (1993,2001) proposed testing if an importance sampler's weights have a finite variance using the [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF] estimator. The Hill estimator is an alternative to the generalized Pareto distribution for estimating the tail parameter ξ. It is defined as

ξ Hill = 1 k k j=1 log w (N -j+1) -log w (N -k) , (5)
where w (1) ≤ w (2) ≤ ... ≤ w (N ) are the order statistics from a sample of importance weights.

The tuning parameter k determines the number of large order statistics used in (5) meaning that it plays a similar role as the threshold parameter u above.

The Hill estimator is largely regarded as a nonparametric estimator although it can be interpreted as a maximum likelihood estimator under appropriate conditions on the distribution function F . Asymptotic normality of the Hill estimator has been studied by [START_REF] Haeusler | On asymptotic normality of Hill's estimator for the exponent of regular variation[END_REF] who have shown that it depends on the rate of growth of k relative to N . Monahan used their result

k/ξ 2 ξ Hill -ξ d → N (0, 1),
and chose the values k = 2N 1/3 and k = 4N 1/3 to build a test statistic. Monahan's test has as its null hypothesis that the variance of the importance weights is infinite making it a test of

H 0 : ξ ≥ 1/2 versus H 1 : ξ < 1/2.
In Section 4, we compare Monahan's test with the tests we propose below. To make the comparison easier, we have reversed the null hypothesis of his test to make it consistent with ours.

Proposed methods

We propose testing for an infinite variance by constructing test statistics from the maximum likelihood estimates of the generalized Pareto distribution. Our paper builds on the work of [START_REF] Smith | Maximum likelihood estimation in a class of non-regular cases[END_REF][START_REF] Smith | Estimating tails of probability distributions[END_REF] who studied the properties of the maximum likelihood estimator of θ = (ξ, β) ′ in this setting. We construct three test statistics using his results. Formally, the hypotheses we will be interested in deciding between is

H 0 : ξ ≤ 1 2 , and
H 1 : ξ > 1 2 .
The null hypothesis implies the existence of the variance while the alternative denies it. In practice it is helpful to make the null a point hypothesis (see, for example, Cox and Hinkley

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT (1974, pp. 331-334)), so making the comparison between

H 0 : ξ = 1 2 , and
H 1 : ξ > 1 2 . ( 6
)
As we have a parametric model for the weights over a threshold, we will use a likelihood function to carry out the testing. This approach to inference is often called the peaks over threshold method. We will study the behavior of the score and likelihood ratio tests of the hypothesis.

Estimation and Wald test

Starting with N importance weights, we first find the n importance weights that lie above the threshold and we construct the excesses z 1 , ..., z n from these large importance weights. The log-likelihood for a sample z 1 , ..., z n over the threshold u equals

log f (z; θ) = -n log β -1 + 1 ξ n i=1 log 1 + ξβ -1 z i . (7)
Maximum likelihood estimation of the parameter vector θ = (ξ, β) ′ is discussed in [START_REF] Smith | Maximum likelihood estimation in a class of non-regular cases[END_REF][START_REF] Smith | Estimating tails of probability distributions[END_REF]. The asymptotic distribution of the maximum likelihood estimator θ is given by

√ n θ -θ d → N 0, I (θ) -1 where I (θ) -1 = (1 + ξ) 2β 2 β β 1 + ξ , (8)
where nI (θ) denotes the expected information matrix. [START_REF] Smith | Maximum likelihood estimation in a class of non-regular cases[END_REF] has shown that likelihood inference is standard for this problem as long as ξ > -1/2. This covers the null hypothesis value H 0 : ξ = 1 2 and all the values under the alternative. [START_REF] Davison | Models for exceedances over high thresholds (with discussion)[END_REF] describe how to reduce maximum likelihood estimation of (7) to a one-dimensional parameter search.

Substituting β = ξ/τ into (7) and differentiating with respect to ξ and τ produces

ξ = ξ (τ ) = 1 n n i=1 log (1 + τ z i ) , (9)
n τ = 1 + 1 ξ (τ ) n i=1 z i 1 + τ z i . (10)
Equation ( 9) can then be substituted into (10) and solved numerically for an estimate τ . Given τ , the maximum likelihood estimates ξ and β are easily obtained. A sensible starting value for τ = ξ/β can be found by imposing the null hypothesis ξ = 0.5 and noting that the expected value of Z is

E[Z] = β ξ 2 Γ(ξ -1 -1) Γ(1 + ξ -1 ) = 2β, for ξ = 0.5. ( 11
)
Iterations can therefore start at τ = 1/z where z = n -1 n i=1 z i .

Once the maximum likelihood estimator of θ is found the Wald test can be computed to test the null hypothesis (6). As an alternative to the Wald test formulation, we can consider the

t-statistic t = n 3 β2 ξ - 1 2 , ( 12
)
where β is the maximum likelihood estimate of β. The null hypothesis is rejected when the t-statistic takes a large positive value compared to a standard normal random variable.

Estimation under null hypothesis and score test

The maximum likelihood estimator of β under the null hypothesis of ξ = 0.5 can be found by univariate Fisher scoring. The score value of β is given by

s β r = - n β + 3 2 β -2 n i=1 z i 1 + 0.5β -1 z i , = nβ -1 3 n n i=1 z i 2β + z i -1 ,
Using the results in ( 11) once again, we can take 0.5z as the initial value for β. After convergence the restricted estimate of β is obtained which we denote by βr .

The one-sided score statistic will be used for testing the null hypothesis as this will be computationally simple. It is based on the score value of ξ under the null hypothesis and is given by

s ξ r = 4 n i=1 log 1 + 0.5β -1 z i -3β -1 n i=1 z i 1 + 0.5β -1 z i = 4 n i=1 log 1 + z i 2β - 6 n n i=1 z i 2β + z i , for ξ = 0.5.
This score value for the null hypothesis ( 6) is a function of β and it can be evaluated when β is replaced by its (restricted) maximum likelihood estimator βr . We know from the results of [START_REF] Smith | Maximum likelihood estimation in a class of non-regular cases[END_REF][START_REF] Smith | Estimating tails of probability distributions[END_REF] that

s ξ * = 1 √ 2n s ξ r d → N (0, 1), (13)
which provides a simple test. We reject the null H 0 : ξ = 1 2 when ŝξ * is significantly positive and where ŝξ * is s ξ * with β replaced by βr .

Likelihood ratio test

When estimation of θ has taken place under the constraint that ξ ≥ 0.5, to deliver θ, together with estimation under the null hypothesis, the likelihood ratio statistic can also be used to test

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT the null LR = 2 log f (z; θ) -log f (z; βr , ξ = 0.5) .
Of course, the null is on the boundary next to the alternative and so the limiting distribution of the statistic is

LR d → 0.5 χ 2 0 + χ 2 1 , under H 0 .
The χ 2 0 distribution is a degenerate distribution with all its probability mass at zero. There exists a 0.5 probability that the likelihood ratio statistic will be zero (as ξ will be negative). The rest of the time the statistic will have a χ 2 1 distribution (e.g. [START_REF] Chernoff | On the distribution of the likelihood ratio[END_REF] and [START_REF] Gourieroux | Kuhn-Tucker, likelihood ratio and Wald tests for nonlinear models with inequality constraints on the parameters[END_REF]).

Monte Carlo evidence

To check the effectiveness of extreme value theory in practice for testing the assumptions behind importance sampling we consider two experiments. The first is based on sampling from a normal density to approximate a normal density with a different variance. Our second experiment is the classic example of estimating a Student's t distribution using a normal candidate density.

Normal density experiment

We take f (x) to be the normal density N (0, 1) and the importance sampling density g(x) is N {0, (1 + ǫ) -1 } where ǫ is a positive number. The function h(x) is constant and equal to one.

Thus for this problem I = 1 in (1). The question is whether the variance of the importance weights exists? The weight function is given by

w(x) = h(x)f (x) g(x) = 1 √ 1 + ǫ exp ǫ 2 x 2 , which obviously has E[w(x)] = 1, while E[w(x) 2 ] = 1 √ 2π √ 1 + ǫ 1 + ǫ ∞ -∞ exp ǫx 2 exp[- 1 2 x 2 (1 + ǫ)]dx = 1 √ 2π 1 √ 1 + ǫ ∞ -∞ exp[- 1 2 x 2 (1 -ǫ)]dx,
which is bounded only if ǫ < 1. The variance of the importance weights will not exist if ǫ ≥ 1, meaning that the central limit theorem will not hold for importance samplers built in that range of values.

<INSERT TABLE 1 HERE>

In Table 1, we present the diagnostic tests for a large number of different thresholds. They

show that for small values of ǫ importance sampling is valid while for values of ǫ ≥ 1 the variance

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

of the importance weights does not exist. For our test statistics to be effective, we require a simulation sample size N that is sufficiently large (say, 100, 000). The Wald test appears to have the most power. It can potentially be oversized if the threshold is selected too high. It can also perform poorly as the threshold gets extremely small. Tests that are slightly oversized are not necessarily problematic. Rejecting importance samplers whose weights are barely finite may not be a bad practice as these algorithms can still be unstable. The LR statistic has nearly the same power and it behaves more robustly for different values of u than the Wald statistic.

This implies that computing more than one statistic is a good idea. Both of these statistics are considerably more powerful than the Monahan statistic computed using either k = 4N 1/3

or k = 2N 1/3 . We note that the Monahan statistic performs better than our tests when the number of draws N is small (e.g. N ≤ 5, 000).

We have chosen the threshold u in an arguably ad hoc manner to both facilitate the Monte Carlo study and to understand how different values of u impact the tests. In practice, the threshold can be chosen in alternative ways and users may want to compare different selection methods. Methods that are common in the extreme value literature are a mean residual life plot or by minimizing the distance between the empirical distribution function and the generalized Pareto distribution evaluated at the maximum likelihood estimates; see, e.g. [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF] and [START_REF] Gonzalo | Which extreme values are really extreme[END_REF]. We recommend graphing the maximum likelihood estimates of ξ along with confidence bands for a large number of thresholds (say 50 values from 0.01N to 0.50N ). We illustrate this strategy in our application to the stochastic volatility model in Section 5.

Our proposed test statistics do not provide strict guarantees for a successful importance sampling procedure. We believe that graphical diagnostics, although informal, can play a com- 4.2 Normal approximation to the Student's t distribution

<INSERT TABLE 2 HERE>

Our second Monte Carlo experiment is the classic example of approximating a Student's t distribution with a normal distribution. This example can be used to evaluate the power of the tests. A similar experiment was also considered by [START_REF] Monahan | Testing the behaviour of importance sampling weights[END_REF]. We let f (x) be a Student's t distribution with ν degrees of freedom and where g(x) is a normal density with mean zero and variance ν/(ν + 1).3 Once again the function h(x) is constant and equal to one. The weight function is simply

w(x) = Γ ν+1 2 √ 2 Γ ν 2 √ ν + 1 1 + x 2 /ν -(ν+1)/2
exp -0.5( ν+1 ν )x 2 . In this example, the quantity E[w(x) 2 ] does not exist.

Table 2 contains the diagnostic tests for different values of the degrees of freedom and the threshold. The performance of the tests are mostly in agreement with the results from the normal experiment above. The Wald test is the most powerful and performs the best when N is reasonably large. The difference in power between this test and the others is significant when ν ≥ 20. Figure 2 plots the series of graphical diagnostics for two values of the degrees of freedom. The recursive variance estimates easily detect problems when ν = 5 while ν = 20 also appears doubtful. Problems created by large values of ν can be hard to recognize with plots, making the combination of test statistics and graphical diagnostics valuable.

<INSERT FIGURE 2 HERE>

An illustration: stochastic volatility model

Stochastic volatility (SV) models have gained considerable interest in econometrics as a class of models which allows volatility to change through time. These models are often stated in continuous time and appear frequently in the option pricing literature. [START_REF] Ghysels | Stochastic volatility[END_REF] and Shephard (2005) review the literature on this topic from an econometric viewpoint.

The likelihood function for SV models is an integral than cannot be solved analytically but it can

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

be evaluated by Monte Carlo simulation. Papers that carry out maximum likelihood analysis of SV models using importance sampling include [START_REF] Danielsson | Stochastic volatility in asset prices: estimation with simulated maximum likelihood[END_REF], [START_REF] Shephard | Likelihood analysis of non-Gaussian measurement time series[END_REF], [START_REF] Durbin | Monte Carlo maximum likelihood estimation of non-Gaussian state space models[END_REF], [START_REF] Sandmann | Estimation of stochastic volatility models via Monte Carlo maximum likelihood[END_REF], [START_REF] Durham | Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes (with discussion)[END_REF], [START_REF] Durham | Likelihood-based specification analysis of continuous-time models of the short-term interest rate[END_REF], [START_REF] Liesenfeld | Univariate and multivariate stochastic volatility models: estimation and diagnostics[END_REF], [START_REF] Brandt | On the relationship between the conditional mean and volatility of stock returns: a latent VAR approach[END_REF], [START_REF] Durham | Monte Carlo methods for estimating, smoothing, and filtering one-and two-factor stochastic volatility models[END_REF], [START_REF] Asai | Asymmetric multivariate stochastic volatility[END_REF], [START_REF] Richard | Efficient high-dimensional importance sampling[END_REF], and [START_REF] Jungbacker | Monte Carlo estimation for nonlinear non-Gaussian state space models[END_REF].

The SV model we use to illustrate our discussion will be a discrete time, log-normal, univariate SV model for a sequence of returns y t for t = 1, ..., T . It can be thought of as an Euler discretisation of a continuous time SV where the spot volatility follows a log-normal Ornstein-Uhlenbeck process (e.g. [START_REF] Hull | The pricing of options on assets with stochastic volatilities[END_REF]). It is given by

y t = µ y + exp(α t /2)ε t , α t+1 = µ α + φ (α t -µ α ) + σ η η t , α 1 ∼ N {0, σ 2 η /(1 -φ 2 )}, (14)
where the disturbances ε t and η t are independent of α t for t = 1, . . . , T . The unconditional mean of the returns is given by µ y while the mean of the log-volatility is µ α . The degree of volatility persistence is measured by the autoregressive parameter φ and the parameter σ η is the standard deviation of the disturbances to the log-volatility process. We will consider two cases for ε t . The first model, denoted by SV-N , assumes that ε t ∼ N (0, 1) while the second model, denoted SV-t, has ε t ∼ t ν where ε t is scaled to have unit variance.

Both SV models require the estimation of the parameter vector Θ. Inference for the unknown parameters will be based on f (y; Θ) where y = (y 1 , ..., y T ) ′ . Clearly

f (y; Θ) = X T f (y|α; Θ)f (α; Θ)dα, (15)
where α = (α 1 , ..., α T ) ′ . The analytic form of f (y|α; Θ) and f (α; Θ) can be deduced from ( 14), but the integral in (15) can only be solved numerically due to the non-Gaussian nature of the model. Our approximation of f (y; Θ) by importance sampling is given by

f (y; Θ) = 1 N N i=1 w i , where w i = w(α i ), (16)
with

w(α) = f (y|α; Θ)f (α; Θ) g(α|y; Θ) and α i ∼g(α|y; Θ).
It is assumed that g(α|y; Θ) is strictly positive for all α ∈ X T . We take the importance sampler g(α|y; Θ) as the Laplace approximation to the posterior of α 1 , . . . , α T given the data y 1 , . . . , y T ; see, e.g., [START_REF] Shephard | Likelihood analysis of non-Gaussian measurement time series[END_REF] and Robert and Casella (2004, p. 107). The initial value of α 1 is a draw from the unconditional distribution of α t . The maximum likelihood estimator of Θ is referred to as Θ. When maximising the log-likelihood function, we use N = 3, 000 draws where the same underlying uniform random variables are used for each likelihood evaluation.

In We start by considering the SV-N model whose estimated parameters are reported in Table 3 together with their 95% confidence intervals. We estimated both models on two different sample sizes, the full sample and a smaller sample corresponding to 2nd January 2002 through 31st

July 2007. To illustrate our diagnostic procedure, we take the simulated maximum likelihood estimate Θ = Θ as fixed. Conditional on this estimate, we increase the simulation size from N = 3, 000 to N + = 100, 000 and produce a set of diagnostic graphs. The question we ask is if the variance of the weights exists for Θ = Θ, basing the diagnostic graphics on these i.i.d.

100, 000 draws. The largest 100 importance weights are presented in Figure 3 together with a histogram of the remaining 99, 900 smaller weights. The third graph presents recursive estimates of the variance of the importance weights. The final graph presents estimates of ξ along with an asymmetric 95% confidence interval for different thresholds u ranging from the largest 1 percent to 50 percent of the N + importance weights.

<INSERT FIGURE 3 HERE>

These diagnostic graphs provide evidence that the importance sampler is unreliable as they indicate a maximum likelihood estimate of ξ which is larger than 0.5. The variance of the importance weights does not exist for the SV-N model. This result is in line with the more informal assessment given in Figure 3(iii). This records the recursive variance estimator, which contains several large jumps. Standard asymptotics to measure the uncertainty of the importance sampler's estimate of the log-likelihood function is problematic in this case.

We now turn our attention to the empirically more interesting SV-t model and look at the diagnostic graphs of the model with estimated parameters. The estimated parameter values are provided in Table 3 and show a large improvement in the fit of the model over the full sample.

Figure 4 presents the same diagnostic graphs as in Figure 3. It is clear that these graphs indicate a better behaved importance sampler. In particular, the estimated shape parameters provide

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

evidence that two moments exist for the importance weights. When these diagnostic graphs are compared with the same graphs for the estimated SV-N model, we conclude that the importance sampler for the SV-t model is more reliable.

<INSERT FIGURE 4 HERE>

In Table 4, we present the test statistics for both SV models. The test statistics are computed for two different thresholds where we have chosen u so that we have the 10,000 and 40,000 largest weights, respectively. The score and t-value statistics are signed tests for the null hypothesis

H 0 : ξ = 1 2 .
We suspect that no variance exists for the importance weights when these tests have significant positive values. Both tests are asymptotically standard normally distributed, so having a 95% critical value of 1.64. The LR test follows an asymptotic 0.5 χ 2 0 + χ 2 1 distribution with the 95% critical value 2.69.

<INSERT TABLE 4 HERE>

The results given in Table 4 are consistent across the choice of statistic. For the SV-N model, the results are poor with all the statistics strongly rejecting the existence of a variance for the full sample. As we reduce the number of observations and the dimension of the integral decreases, the importance function for the SV-N begins to perform better. The point estimates in Table 3 appear reasonable for the SV-N model, despite the fact that the weight function has an infinite variance. We believe users should still proceed with caution. Robert and Casella (2004, p. 95) emphasize that when the variance does not exist the range of Monte Carlo estimates computed using different sequences of random numbers will be enormous. For an SV model, this means that the estimated log-volatilities may be far from the true, unobserved values. [START_REF] Robert | Stochastic volatility: an experimental approach[END_REF] demonstrates how an importance sampler with a poorly chosen importance density can entirely miss-estimate the true log-volatilities for many sequences of random numbers.

Conclusion

Importance sampling is increasingly used in econometrics to estimate likelihoods and to compute moments of the posterior distribution in Bayesian analyses. Behind its use lies the assumption that the variance of the importance sampling weights exists. In this paper, we have proposed test (iv) estimated shape parameters with 95% confidence intervals (dotted lines) for 50 different thresholds going from the largest 1% to 50% of the N + weights.

  plementary role. They are a simple way to detect severely deficient importance sampling algorithms. A graph of the largest 100 weights indicates the seriousness of the outliers. The remaining weights can be represented via a histogram. Plotting estimates of the variance of the weights recursively can show the impact that each weight has on the estimate. As the sample increases, the variance should converge to a constant if it exists. Plots like these do not always detect problems and we recommend viewing them for many different sequences of random numbers.<INSERT FIGURE 1 HERE>Graphs for the normal density experiment with ǫ = 0.5, 1.0, and 1.2 are presented in the three rows of Figure1. The values of ǫ were chosen to compare the output from importance samplers representing both the null and alternative hypotheses. When ǫ = 0.5, the variance is finite and the recursive estimate of the variance convergences smoothly. The graphs for the cases ǫ = 1.0 and 1.2 are clearly informative in this example. The large jumps in the recursive variance estimates correspond to the outliers in the left hand graphs of the middle and bottom rows.

Figure 1 :Figure 2 :Figure 3 :Figure 4 :

 1234 Figure1: Graphical diagnostics for the normal experiment. Top, middle, and bottom rows represent ǫ = 0.5, 1.0, 1.2 respectively. The left hand side pictures are of the largest 100 weights, the middle pictures are a histogram of the remaining weights and the right hand side shows a recursive estimator of the variance of the weights. N = 100, 000 draws.

  our empirical study we use the return series y t of daily Standard & Poor's 500 stock index closures. The historical return series is for the period 3rd January 1985 to 31st July 2007 and was obtained from Yahoo Finance. Our sample consists of 5, 693 daily observations. Importance sampling is being carried out over 5, 693 dimensions in this case. The continuously compounded (raw) on the stock index are not adjusted for dividends and they are expressed in percentage terms and therefore given by y t = 100(ln P t -ln P t-1 ) where P t denotes the closing price of the Standard & Poor's 500 index on day t.

	<INSERT TABLE 3 HERE>
	5.1 Test statistics and graphical diagnostics

  Table 4 also includes test statistics computed from models estimated on roughly one-seventh the sample size (818 observations). At this point, the test statistics fail to reject the null hypothesis. Conversely, the test statistics indicate that the importance function for the SV-t model has a finite variance for all sample sizes.

This self-normalized estimator induces a small bias but it remains consistent.Liu (2001, p.

35) andRobert and Casella (2004, p.95) note that the self-normalized estimator may be preferred as it typically has smaller mean square error than the unbiased estimator.

We believe we have implemented the same experiment as Monahan. A potential missprint in his paper may imply an alternative normal importance density with variance (ν + 1)/ν rather than the importance density with variance ν/(ν + 1) given here.

Acknowledgements

We would like to thank the Editor and three referees for their insightful comments. Furthermore, we thank Shimrit Abraham for excellent computational support and Charles S. Bos for helpful All the computations reported in this paper were carried out using the OxMetrix 5.0 programming environment of Doornik (2007).

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

statistics and graphical diagnostics to check this assumption in high dimensional settings when it is often difficult to evaluate analytically. We have applied our new tests on empirical problems and find this assumption is often overturned. This suggests that this is an important problem of which researchers should be aware. We hope that our methods will provide researchers with a better understanding of when their importance samplers are reliable. u = 0.55N ǫ = 0.5 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 ǫ = 0.8 0.00 0.46 0.07 0.00 0.00 0.00 0.60 0.21 0.00 0.00 ǫ = 1.0 0.13 0.97 0.76 0.01 0.02 1.00 1.00 1.00 0.00 0.00 ǫ = 1.2 0.71 1.00 0.99 0.05 0.06 1.00 1.00 1.00 0.08 0.08 ǫ = 1.5 0.99 1.00 1.00 0.23 0.19 1.00 1.00 1.00 0.46 0.33 ǫ = 3.0 1.00 1.00 1.00 0.90 0.73 1.00 1.00 1.00 1.00 0.95 ǫ = 5.0 1.00 1.00 1.00 0.98 0.89 1.00 1.00 1.00 1.00 1.00 u = 0.50N ǫ = 0.5 0.00 0.00 0.00 --0.00 0.00 0.00 -ǫ = 0.8 0.00 0.19 0.01 --0.00 0.02 0.00 -ǫ = 1.0 0.03 0.87 0.47 --0.93 1.00 0.99 -ǫ = 1.2 0.41 1.00 0.94 --1.00 1.00 1.00 -ǫ = 1.5 0.95 1.00 1.00 --1.00 1.00 1.00 -ǫ = 3.0 1.00 1.00 1.00 --1.00 1.00 1.00 -ǫ = 5.0 1.00 1.00 1.00 --1.00 1.00 1.00 -u = 0.40N ǫ = 0.8 0.00 0.02 0.00 --0.00 0.00 0.00 -ǫ = 1.0 0.00 0.48 0.11 --0.02 0.83 0.42 -ǫ = 1.2 0.08 0.92 0.63 --0.99 1.00 1.00 -ǫ = 1.5 0.64 1.00 0.98 --1.00 1.00 1.00 -ǫ = 3.0 0.82 1.00 1.00 --1.00 1.00 1.00 -ǫ = 5.0 0.97 1.00 1.00 --1.00 1.00 1.00 -u = 0.30N ǫ = 0.8 0.00 0.01 0.00 --0.00 0.00 0.00 -ǫ = 1.0 0.00 0.20 0.03 --0.00 0.10 0.01 -ǫ = 1.2 0.01 0.68 0.30 --0.55 0.99 0.97 -ǫ = 1.5 0.28 0.97 0.85 --1.00 1.00 1.00 -ǫ = 3.0 0.91 1.00 0.98 --1.00 1.00 1.00 -ǫ = 5.0 1.00 1.00 1.00 --1.00 1.00 1.00 -u = 0.10N ǫ = 0.8 0.00 0.01 0.00 --0.00 0.00 0.00 -ǫ = 1.0 0.00 0.05 0.01 --0.00 0.00 0.00 -ǫ = 1.2 0.00 0.19 0.08 --0.01 0.37 0.22 -ǫ = 1.5 0.03 0.47 0.34 --0.72 0.98 0.99 -ǫ = 3.0 0.75 0.98 0.76 --1.00 1.00 1.00 -ǫ = 5.0 0.97 1.00 1.00 --1.00 1.00 1.00 -u = 0.01N ǫ = 0.8 0.00 0.01 0.01 --0.00 0.00 0.00 -ǫ = 1.0 0.00 0.01 0.02 --0.00 0.00 0.01 -ǫ = 1.2 0.00 0.00 0.05 --0.00 0.00 0.08 -ǫ = 1.5 0.01 0.00 0.10 --0.04 0.00 0.37 -ǫ = 3.0 0.08 0.00 0.34 --0.81 0.01 0.99 -ǫ = 5.0 0.17 0.00 0.51 --0.98 0.02 1.00 -- 1.00 1.00 1.00 0.00 0.00 1.00 1.00 1.00 0.00 0.00 ν = 100 1.00 1.00 1.00 0.00 0.00 1.00 1.00 1.00 0.00 0.00 u = 0.10N ν = 1 0.98 1.00 1.00 --1.00 1.00 1.00 -ν = 2 0.89 1.00 1.00 --1.00 1.00 1.00 -ν = 5 0.37 1.00 0.87 --1.00 1.00 1.00 -ν = 20 0.00 0.50 0.05 --0.00 0.46 0.06 -ν = 100 0.00 0.00 0.00 --0.00 0.00 0.00 -u = 0.01N ν = 1 0.21 0.02 0.56 --1.00 0.90 1.00 -ν = 2 0.15 0.05 0.47 --0.97 0.84 1.00 -ν = 5 0.04 0.15 0.24 --0.49 0.78 0.91 -ν = 20 0.00 0.14 0.01 --0.00 0.02 0.00 -ν = 100 0.00 0.01 0.00 --0.00 0.00 0.00 --

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT