Marie- Laurence Grespier

Poncet

Vo Van Chuong

Ngoc Linh

Hong Dan

Minh Phuong

Loic Piffet

Sébastient Dutercq

Thuy Nga

Xuan Lan

Hiep Thuan

Trang Dai

Thuy Lynh

Thanh Binh

Xuan Hieu

Van Ly

Plusieurs personnes m'ont aidé durant ce travail de thèse. La première personne que je tiens à vivement remercier est mon Directeur de thèse, le professeur Richard Émilion, pour le choix du sujet, pour sa confiance en moi, sa patience, et son apport considérable sans lequel ces travaux n'auraient pas pu être menés à terme. Je lui suis reconnaissant pour tout le temps qu'il a consacré à répondre à mes questions et à corriger ma rédaction. Ce fut pour moi une expérience extrêmement enrichissante.

J'adresse mes vifs remerciements à mon codirecteur de thèse, le professeur Romain Abraham, Directeur du laboratoire MAPMO, pour le choix du sujet, pour ses explications et ses précieux conseils qui m'ont éclairé, pour son accueil et son aide dans le laboratoire durant toute ces années.

Je tiens à vivement remercier les professeurs Jean-François Delmas et Sophie Dabo-Niang d'avoir accepté et d'accomplir la délicate tâche de rapporteurs de cette thèse.

Mes vifs remerciements aux membres du jury d'avoir accepté d'évaluer ce travail de recherche.

Je remercie très spécialement Dr. Ted Soubdhan, Maître de conférences en physique à l'université d'Antilles-Guyane qui nous a introduit à la problématique

Introduction

Context

The aim of the present thesis is to propose some probabilistic models for sequences of solar radiation which is defined as the energy given off by the sun (W/m 2 ) at the earth suface. Our main model concerns a Stochastic Differential Equations (SDE) in random environment, the latter being modelized by a hidden Markov chain. Statistical fitting of such models hinges on filtering equations that we establish in order to update the estimations in the steps of EM algorithm. Experiments are done using real large datasets recorded by some terrestrial captors that have measured solar radiation. Such a modelling problem is of greatest importance in the domain of renewable energy where short-term and very short-term time horizon prediction is a challenge, particularly in the domain of solar energy.

Random aspects

Probabilistic models turn out to be relevant as the measured solar radiation is actually a global radiation, or total radiation, which results from two components, a deterministic one and a random one, namely -the direct radiation which is the energy coming through a straight line from the sun to a specific geographical position of the earth surface. At a given time this deterministic radiation can be computed quite precisely and as it roughly corresponds to a measurement during a perfectly clear-sky weather, it is also known as the extra-terrestrial radiation -the diffuse radiation which is reflected by the environment and depends on meteorolgical conditions, and is therefore highly random. Both components can be measured by captors.

The total solar radiation can also be studied indirectly by considering its dimensionless form, the so-called clearness index (CI), which is defined as the ratio of the total radiation to the direct radiation and thus is a nice descriptor of the atmospheric transmittance. Our approach will therefore consist in considering a discrete (resp. continous) sequence of solar observations as a path of a discrete-time (resp. continuous-time) stochastic process.

The following figure (Figure 1) illustrates our arguments: the iobserved irregular falls are due to frequent cloud passages which depend on some random conditions such as wind speed, type of clouds and some other meteorological variables:

Two possible approaches

In the present understanding, the establishment of meteorological radiation models are usually based on physical processes as well as on statistical techniques [Gueymard 1993, Kambezidis 1989, Muneer 1997, Psiloglou 2000, Psiloglou 2007].

The physical modelling studies the physical processes occurring in the atmosphere and influencing solar radiation. Accordingly, the solar radiation is absorbed, reflected, or diffused by solid particles in any location of space and especially by the earth, which depends on its arrival for many activities such as weather, climate, agriculture, . . . . The physical calculation method is exclusively based on physical considerations including the geometry of the earth, its distance from the sun, geographical location of any point on the earth, astronomical coordinates, the composition of the atmosphere, . . . . The incoming irradiation at any given point takes different shapes. The second approach, "statistical solar climatology" branches into multiple aspects: modelling of the observed empirical frequency distributions, forecasting of solar radiation values at a given place based on historical data, looking for statistical interrelationships between the main solar irradiation components and other available meteorological parameters such as sunshine duration, cloudiness, temperature, and so on. Our work has clearly taken the second approach.

HMM and SDE

Stochastic characteristics of solar radiation highly depend on some unobserved meteorological events such as frequency, height and type of the clouds and their optical properties, atmospheric aerosols, ground albedo, water vapor, dust and atmospheric turbidity (Woyte et al. (2007)) while a sequence of solar radiation can be observed and measured at a given station. This has suggested us to model a random sequence of clearness index (resp. a stochastic process of solar radiation) by using a HMM which is a pair of correlated stochastic processes: the first (unobserved) one, called the state process, is a finite-state Markov chain in discrete-time (resp. in continoustime) representing meteorological regimes while the second (observed) one depends on the first one and describes the sequence of clearness index (resp. the process of solar radiation) as a discrete process (resp. a continuous one, solution of a SDE). The idea of using HMM and SDE in the study of solar radiation sequences was mentioned by T. Soubdhan and R. Emilion in [START_REF] Soubdhan | [END_REF], Soubdhan 2011]. After a classification of daily solar radiation distributions, the authors thought that the sequence of class labels can be governed by a HMM in discrete time with some underlying unobservable regimes. The same authors have also proposed a SDE to model a continuous-time clearness index sequence but their data-driven approach fails for prediction during high variability regimes. However our work has been developed starting from these ideas. Our results can be summarized as follows:

1. We propose a discrete time HMM to model a daily (resp. hourly, resp. monthly) clearness index sequence.

2. We propose a continuous time HMM to model the clearness index process over a time interval [0, T ] in a solar day.

3. We propose a continuous time HMM and a SDE to directly model the total solar radiation process over a time interval [0, T ].

Estimation procedures

To fit our models to observed data, the estimation procedures will combine Maximum Likelihood Estimators (MLE) and Expectation Maximization (EM) algorithm for partially observed systems [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF], Celeux 1989].

Filtering

A crucial notion in our estimation procedure is that of filter which is a time-indexed increasing family of σ-algebras, each one being generated by the events occured up to time t. The filtering process is defined as the family of conditional expectations w.r.t. these σ-algebras. A large part of our contribution deal with recursive equations of the filtering process needed in the estimation algorithms. They hinge on the work of [START_REF] Dembo | [END_REF][START_REF] Campillo | [END_REF][START_REF] Elliott | [END_REF][START_REF] Elliott | [END_REF]].

Continuous-time filtering equations will be approximated by robust versions, following an approach due to [Clark 1977] and using some results of [START_REF] James | [END_REF][START_REF] Krishnamurthy | [END_REF][START_REF] Clark | [END_REF].

Reference probability method. Girsanov theorem

A great part of our computations concerns the so-called reference probability method which refers to a procedure where a probability measure change is introduced to reformulate the original estimation and control task into a new probability space (fictitious world) in which well-known results for identically and independently distributed (i.i.d.) random variables can be applied. Then the results are reinterpreted back to the original probability space (real world) by applying [Elliott 1995, chp. 1]. The Radon-Nykodim derivative of the new probabily measure w.r.t. the original one is given by the famous Girsanov theorem in both its discrete and continuous time version.

Thesis organization

Our thesis is divided into five Chapters following this introductory part.

Chapter 1.

In the first chapter we present some backgroung notions concerning solar radiation: direct, diffuse and global radiations, clearness index. The computation of the direct solar radiation is detailed. The end of the chapter briefly presents some points concerning measurement devices and datasets the we have dealt with.

Chapter 2.

In the second chapter, we recall some mathematical results that will be needed in chapters 3 and 4: conditional Bayes formula, Ito product, Ito formula, Girsanov theorem, HMM, EM algorithm.

Chapter 3.

In this third chapter we introduce three models for clearness index sequences (CISs): -DTM-K, a Discrete-Time Model for discrete daily CISs, (K h ) h=1,2,... -CTM-k, a Continuous-Time Model for continuous processes of CI, (k t ), t ∈ [0, T ], and its Discrete-Time Approximate Model, DTAM-k, obtained from time dicretization by uniformly partitionning [0, T ] into intervals of width ∆. v For each model, we define the state process, the observation process and the parameter vector. The state process of these models are finite-state homogeneous Markov chains. For CTM-k, the transition matrix of the chain is a rate matrix. For DTAM-k, the ∆ width in the time partition is chosen to be small enough so that the transition matrix of the chain be a stochastic matrix. The observation process is a function of the chain which values are corrupted by a Gaussian noise (for DTM-K and DTAM-k) and by a standard Brownian motion (for CTM-k).

The filtering equations are established with complete proofs. Computations to obtain MLE updating formulas in the iterations of EM algorithm are detailed. Using DTAM-k, we first establish the computable approximation of the continuous time equations in CTM-k, and then we provide the estimates for the noisy variance.

Chapter 3 ends with some experiments with real data. Parameters of DTM-K are estimated from La Réunion island (France) data with daily CISs having similar characteristics while parameters of the CTM-k approximated by parameters of DTAM-k are estimated from Guadeloupe island (France) data which were sampled at 1Hz (i.e. at each second).

Chapter 4.

In this fourth chapter, we propose our main model, a continuous-time HMM for the total solar radiation sequence (y t ) t≥0 under the random effects of meteorological events, denoted CTM-y.

The state process is similar to the CTM-k case but the observation process (y t ) modelling total solar radiation process, is assumed to be of the SDE:

dy t = [g(X t )I t -y t ]dt + σ(X t )y t dW t ,
where I t is the extraterrestrial radiation received at time t, (W t ) is a standard Brownian motion and g(X t ), σ(X t ) are functions of the Markov chain (X t ).

Again, the change-of-measure technique and the steps of EM algorithm establishing the filtering equations for updating the parameter vector g, are fully detailed.

Here too, we propose an approximation of state filter equation and we build a Discrete-Time Approximate Model (DTAM-y) to provide discrete-time approximate equations. Our computations hinge on a robust discretization of continuous-time filters recently obtained by [Elliott 2010, chap. 1]. Estimation of the noisy variance is studied.

Experimentations with real data and parameter estimations are performed from various samples of data sampled at 1Hz. Using the model with estimated parameters, we generate some simulations of solar radiation process paths. 1.6.2 Data observed in Guadeloupe and La Réunion islands . . . . 9

Résumé

Dans ce chapitre, nous rappelons d'abord quelques notions de physique en énergie solaire : rayonnement solaire extraterrestre, calcul du rayonnement direct, rayonnement diffus, rayonnement total ou global, indice de clarté. Nous parlerons brièvement des instruments de mesure du rayonnement et nous décrirons enfin les données réelles que nous avons utilisées.

Figure 1.1: The types of solar radiation.

Introduction

Solar radiation emission from the sun into every corner of space appears in the form of electromagnetic waves that carry energy at the speed of light. The solar radiation is absorbed, reflected, or diffused by solid particles in any location of space and especially by the earth (Figure 1.1). This process depends on many environment conditions such as weather, climate, pollution, . . . . The incoming radiation at any given point takes different shapes depending on its geographical location, its astronomical coordinates, its distance from the sun, the composition of the local atmosphere and the local topgraphy. This section provides some basic concepts, definitions, and astronomical equations which are used in our thesis. These concepts, definitions and equations are referenced from [START_REF] Liu | [END_REF], Psiloglou 2000, Sen 2008, Tovar-Pescador 2008].

1.2 Extraterrestrial solar radiation

Extraterrestrial normal radiation

The extraterrestrial normal radiation, denoted I 0 , also called top of the atmosphere radiation, is the solar radiation arriving at the top of the atmosphere. It can simply be considered as the product of a solar constant denoted by I CS and a correction factor of the earth's orbit, namely its excentricity, denoted by ε:

I 0 = I CS • ε.
(1.1) Figure 1.2: Horizontal plane for the extraterrestrial horizontal radiation.

The World Radiation Center has adopted that I CS = 1367 W/m 2 with an uncertainty of 1% [START_REF][END_REF]] and introducing I CS is justified as follows. As already said, the sun radiation is subject to many absorbing, diffusing, and reflecting effects within the earth's atmosphere which is about 10 km average thick and, therefore, it is necessary to know the power density, i.e., watts per meter per minute on the earth's outer atmosphere and at right angles to the incident radiation. The density defined in this manner is referred to as the solar constant I CS . It is equivalent to the energy from the sun, per unit time, received on a unit area of surface perpendicular to the direction of propagation of the radiation, at mean earth-sun distance, outside of the atmosphere. The excentricity ε, as suggested by [Spencer 1972], is given by ε = 1.00011 + 0.034221 cos Γ + 0.00128 sin Γ + 0.000719 cos 2Γ + 0.000077 sin 2Γ, (1.2) where the day angle Γ (in radians) is equal to:

Γ = 2πn d -1 365 ,
n d denoting the number of the day in the year (1 for first of January, 365 for December 31).

A simple approximation for ε was suggested by [Duffie 1980[START_REF][END_REF]:

ε = 1 + 0.033 cos 2πn d 365 .
(1.3)

Extraterrestrial horizontal radiation

At time t of a day, the amount incident radiation per horizontal surface area unit along the zenith direction, called the extraterrestrial horizontal radiation and denoted I t , is related to the extraterrestrial normal radiation I 0 as follows:

I t = I 0 cos θ z , (1.4)
where θ z is the zenith angle at time t between the normal to the surface and the direction of the direct beam (Figure 1.2). The zenith angle calculation is described in Section 1.3 assuming, for sake of simplicity of modelling, that land is horizontal. The solar day is defined as the time that is needed by the Sun to achieve a complete tour of the Earth [Lanini 2010]. This does not necessarily correspond to 24 hours and varies from year to year. Figure 1.3 illustrates how the relative Sun position is moving across the sky: pictures of the Sun taken by an immobile photographer at the same time of the day have been superimposed. It can be seen that after one whole year of observations, the Sun is computing a eight-shape circuit. The principal causes of this phenomena are the elliptical shape of the terrestrial orbit around the Sun and the tilt of the Earth in relation to the plane of its orbit.

Zenith angle calculation

Equation of time

As a consequence, at 12 noon the Sun does not have the same position in the different months. The curve described by the so called the equation of time, (E t ), was proposed by [Spencer 1972] and then truncated by [Iqbal 1986]: E t = 229.18(0.000075+0.001868 cos Γ -0.032077 sin Γ -0.014615 cos 2Γ -0.04089 sin 2Γ),

(1.5)

where Γ = 2π(n d -1)

365

, n d = 1, 2, . . . , 365.

Apparent solar time

Most meteorological measurements are recorded in terms of local standard time.

In many solar energy calculations, it is necessary to obtain irradiation, wind, and temperature data for the same instant. It is, therefore, necessary to compute local apparent time, which is also called the true solar time. Solar time is the time to be used in all solar geometry calculations. It is necessary to apply the corrections due to the difference between the local longitude, L loc , and the longitude of the standard time meridian, L stm . The apparent time, L at , can be calculated by considering the standard time, L st according to [Iqbal 1986] as:

L at = L st ± (L stm -L loc ) + E t . (1.6)
where E t is caculated as in (1.5).

In this expression, "+" applies to west direction and "-" applies to east direction. All terms in the above equation are to be expressed in hours.

Hour angle

The hour angle, denoted ω, is the angular displacement of the sun east or west of the local meridian due to rotation of the earth on its axis at 15 0 per hour as morning negative and afternoon positive [Iqbal 1986]:

ω = 15(12 -L at ), (1.7)
where L at is the apparent time which is calculated as in (1.6).

Declination

The solar declination, denoted δ, is the angle between a line joining the centers of the Sun and the Earth to the equatorial plane, depends on the date and on the location (Figure 1.4): north direction has positive value, its maximum is equivalent to +23.45 0 at the summer solstice and its minimum to -23.45 0 at the winter solstice. We can consider the following expressions for the approximate calculation of δ ( [Iqbal 1986]) as:

δ = 23.45 sin 360(284 + n d ) 365 , (1.8)
where n d is the day number of year: n d = 1, 2, . . . , 365.

Zenith angle

The zenith angle, denoted θ z , is the angle between the vertical and the line to the sun i.e., the angle of incidence of beam radiation on a horizontal surface (see again Figure 1.2). At solar noon zenith angle is zero, in the sunrise and sunset this angle is 90 0 . The zenith angle θ z can be actually calculated by [Sen 2008, page 86]:

cos θ z = cos φ cos δ cos ω + sin φ sin δ, (1.9)
where φ is the latitude of the location and δ calculated by (1.8), ω calculated by (1.7).

Total solar radiation

Total (global) solar radiation, G t , is the sum of the direct beam, I b , and the diffuse solar radiation, I d , on a horizontal surface (Figure 1.5). Solar radiation from the sun after traveling in space enters the atmosphere at the space-atmosphere interface, where the ionization layer of the atmosphere ends. Afterwards, a certain amount of solar radiation is absorbed by the atmosphere, by the clouds, and by particles in the atmosphere. A certain amount is reflected back into the space, and a certain amount is absorbed by the earth's surface. The combination of reflection, absorption (filtering), refraction, and scattering result in highly dynamic radiation levels at any given location on the earth. As a result of the cloud cover and scattering sunlight, the radiation received at any point is both direct (or beam) and diffuse (or scattered), see again Figure 1.1.

Direct solar radiation

Direct solar radiation is defined as the radiation which travels in a straight line from the sun to the earth's surface. It is the solar radiation received from the sun without scatter by the atmosphere and without any disturbances. The quantity of direct solar radiation reaching any particular parts of the earth's surface is determined by the position of the point, time of year, shape of the surface, . . . . To model this would require knowledge of intensities and direction at different times of the day . . . . As examples, we can refer to the models in [Psiloglou 2000, Psiloglou 2007].

Diffuse solar radiation

After the solar radiation enters the earth's atmosphere, it is partially scattered and partially absorbed. The scattered radiation is called diffuse radiation. Again, a portion of this diffuse radiation goes back to space and a portion reaches the ground.

Diffuse radiation is first intercepted by the constituents of the air such as water vapor, CO 2 , dust, aerosols, clouds, etc., and then it is released as scattered radiation in many directions. This is the main reason why diffuse radiation scattering in all directions and closed to the earth's surface as a source does not give rise to sharp shadows. When the solar radiation in the form of an electromagnetic wave hits a particle, a part of the incident energy is scattered in all directions and is called diffuse radiation. Diffuse radiation is scattered out of the solar beam by gases (Rayleigh scattering) and by aerosols (which include dust particles, as well as sulfate particles, soot, sea salty particles, pollen, etc.). The Reflected radiation is mainly reflected from the terrain and is therefore more important in mountainous areas.

Diffuse radiation occurs when small particles and gas molecules diffuse part of the incoming solar radiation in random directions without any alteration in the wavelength of the electromagnetic energy. Diffuse cloud radiation would require modeling of clouds and this is considered as quite impossible because of a great daily variability.

Clearness index

The ratio of the total solar radiation G t to the extraterrestrial horizontal radiation I t is defined as the clearness index and is denoted k t :

k t = G t I t , (1.10)
Clearness index is the quantity needed to focus on the analysis of fluctuations in solar radiance. It gives the ratio of the actual energy on the ground to that initially available at the top of the atmosphere accounting, therefore evaluating at time t the transparency of the atmosphere. Alternatively, this index can be considered as an instantaneous class membership degree, the class being an ideal perfect clear-sky day, the more this index is closed to one, the more the day is clear at time t.

For long-term predictions, the clearness index is often considered over a given time interval ∆t. It is denoted by K ∆t and is defined as the relation between the horizontal total radiation on the ground and the extraterrestrial horizontal radiation over the same time interval ∆t:

K ∆t = ∆t G s ds ∆t I s ds .
(1.11)

The usually used integration periods are the day and the hour, termed daily clearness index and hourly clearness index, respectively.

Solar radiation measurement

This section is designed to be a concise introduction for the instrumentation used to measure the components of solar radiation as well as for the climatic characteristics and geographical location of the areas where the observed data were recorded.

Solar radiometers

Pyrheliometer

The Pyrheliometer is a solar radiometer which is used to measure the "direct normal radiation" I bn (note that I b = I bn cos θ z ). Pyrheliometers have a narrow aperture (generally between 5 0 and 6 0 total solid angle), admitting only beam radiation with some inadvertent circumsolar contribution from the Sun's aureole within the field of view of the instrument, but still excluding all diffuse radiation from the rest of the sky. Pyrheliometers must be pointed at, and track the Sun throughout the day. Their sensor is always normal to the direct beam, so that I bn is often called "direct normal radiation" (Figure 1

Data observed in Guadeloupe and La Réunion islands

The total solar radiation measurements used in our estimation procedures were performed in two French islands, namely Guadeloupe and La Réunion, located in the West Indies and the Indian Ocean, respectively. These areas are exposed to an important solar radiation and are characterized by a humid tropical climate. Guadeloupe island is located at 16 0 15N latitude and 60 0 30W longitude. The average solar load for a horizontal surface is between 4 kWh/m 2 and 7 kWh/m 2 per day. The air temperature varies between 17 0 C and 33 0 C. Relative humidity ranges from 70% to 80% and the trade winds are relatively constant throughout the year. The total solar radiation measurements were performed in this island in 2006 by a Pyranometer from KIPP&ZONEN, model SP-Lite, a sensor having a response time inferior to 1 s. The SP-Lite measures the solar energy received from the entire hemisphere (180 0 field of view). These data were recorded, pretreated, analyzed ans interpreted by Dr. T. Soubdhan, Assistant Professor in Physics at University of Antilles-Guyane, Guadeloupe.

La Réunion is a Southern Hemisphere volcanic island with an average temperature oscillating from 24 0 C to 32 0 C in the coastal regions and oscillating from 15 0 C to 22 0 C in the regions located above an altitude of 1500 m in the interior of the island. The combination of a very steep terrain, with large variations in altitude, and prevailing trade winds from south-southeast induce local contrasts in weather patterns at ground level. The radiations were captured by a SPN1 Pyranometer ([Delta-T-Devices 2012]), a sensor rated as a "good quality" one by World Meteorological Organization. This sensor is actually based on a set of seven thermopiles, symmetrically arranged below a shadow dome according to a specific geometry, ensuring by that way that, at any time of the day, wherever in the world the measurement is made, there is always one sensor fully exposed to the sun and one sensor fully shadowed. The recorded data used for our daily clearness index sequences modelling were sampled at 0.1 s at Moufia campus location (20 0 54S, 55 0 29E) from 2009 to 2011 in the setting of a La Réunion region project titled RCI-GS. They were managed by software engineer M. Delsaut and were averaged to give one collected point per minute for final storage purpose.

Both data collections and storages used a datalogger from CAMMPBELL SCI-ENTIFIC (the "burning" sunshine recorders were first developed by John Francis Campbell in 1853 and later modified in 1879 by Sir George Gabriel Stokes, . . . ).

Chapter 2

Mathematical recalls

Résumé

Nous supposons connu les ésprances conditionnelles, les chaînes de Markov à temps discret et à temps continu, les martingales et le calcul stochastique. Nous précisons dans ce chapitre formalisme, notations et méthodes utilisés dans les prochains chapitres en suivant la présentation de [START_REF] Elliott | [END_REF]] : modèles Markoviens cachés (HMM) à temps discret et méthode de changement de probabilité dite méthode de la probabilité référente basée sur le théoreme de Girsanov en version discrète, rappels de certains résultats de calcul stochastique, HMM à temps continu et méthode de la probabilité référente, algorithme EM en temps continu.

Abstract

Conditional expectations, discrete-time and continuous time Markov chains, martingales and stochastic calculus are assumed to be known. We precise in this chapter the formalism, the notations and the methods used in the next chapters, following the presentation in [START_REF] Elliott | [END_REF]] : Hidden Markovian Models (HMM) in discrete time and reference probability method based on a discrete version of Girsanov theorem, recalls of some results in stochastic calculus, continuous time HMM and reference probability method, EM algorithm in continuous time.

Conditional expectations

Radon-Nikodym derivative

Theorem 2.1. (Radon-Nikodym) If P and P are two probability measures on (Ω, B) such that for each B ∈ B, P (B) = 0 implies P (B) = 0, then there are exists a nonnegative random variable Λ, such that P (C) = C ΛdP for all C ∈ B. We write dP /dP | B = Λ.

For a proof, see [START_REF] Wong | [END_REF]].

Definition 2.1. Let X ∈ L 1 and A be a sub-σ-field of B. If X is non-negative and integrable we can use the Radon-Nikodym to deduce the existence of an A-measurable random variable, denoted by E(X|A), which is uniquely determined except on an even of probability zero, such that

A XdP = A E(X|A)dP, (2.1)
for all A ∈ A.

E(X|A) is called the conditional expectation of X give A. For a general integrable random variable we define E(X|A) as E(X + |A) -E(X -|A).

The following is a list of classical results. If A 1 and A 2 are two sub-σ-fields of B such that A 1 ⊂ A 2 , then

E(E(X|A 1 )|A 2 ) = E(E(X|A 2 )|A 1 ) = E(X|A 1 ). (2.2) If X, Y, XY ∈ L 1 , and Y is A-measurable, then E(XY |A) = Y E(X|A). (2.3) If X and Y are independent, then E(Y |σ(X)) = E(Y ).
(2.4)

Jensen inequality

Theorem 2.2. Let {Ω, F, P } be a probability space and G a subfield of F. Let φ : R → R be convex and let X be an integrable random variable such that φ(X) is integrable. Then, we have

φ(E(X|G)) ≤ E(φ(X)|G).
A proof can be found in [Elliott 1982].

Conditional Bayes formula

Theorem 2.3. Suppose (Ω, F, P ) is a probability space and G ⊂ F is a sub-σ-field. Suppose P is another probability measure absolutely continuous with respect to P and with Radon-Nikodym derivative dP /dP = Λ. Then if φ is any P integrable random variable

E[φ|G] = ψ where ψ = E[Λφ|G] E[Λ|G] if E[Λ|G] > 0
and ψ = 0 otherwise.

Theorem 2.3 was proved by [Elliott 2010, Theorem 3.2].

Martingale difference sequence

Definition 2.2. Let (X h ) = {X h , h = 1, 2, . . . } be an adapted discrete stochastic process on a filtered probability space (Ω, F, (F h ), P ). (X h ) is called a martingale difference sequence (MDS) if it satisfies the following two conditions:

i. E(|X h |) < ∞,
ii. E(X h |F h-1 ) = 0, a.s., for all h.

So, a stochastic serie is an MDS if its expectation with respect to the past is zero. MDS is an extremely useful concept in modern probability theory because it implies much milder restrictions on the memory of the sequence than independence. Most of limit theorems that hold for an independent sequence also hold for an MDS.

The definition implies that:

• if (Y h ) is a martingale, then (U h ) = (Y h -Y h-1
) will be an MDS,

• (V h ) = (X h -E[X h |F h-1 ]
) is an MDS.

The sequences (U h ), (V h ) above are also called sequences of martingale increments.

Binary vector representation of a Markov chain

Consider a (discrete or continuous time) Markov chain (X time ) with finite state space S X = {s 1 , s 2 , ..., s N }.

Let 1 (X time =s i ) denote the so-called indicator function defined as 1 (X time =s i ) = 1 if X time = s i and 1 (X time =s i ) = 0 if X time = s i . The prime symbol denoting transpose, we consider the following binary vector representation:

X time = (1 (X time =s 1 ) , 1 (X time =s 2 ) , ..., 1 (X time =s N ) ) ′ ,
so that at any time, just one component of X time is one while the others are zero.

Then, for sake of simplicity in computations, we will now consider the chain (X time ) which is derived from (X time ), the state space of (X time ) being the set of unit vectors e i with all components 0 but 1 at the i-th component:

S = {e 1 , e 2 , ..., e N }.
Let • , • denote the usual inner product in R N . Noticing that X time , e i = 1 (X time =s i ) , i = 1, 2, . . . , N , we will write the vector representation of the Markov chain as X time = ( X time , e 1 , X time , e 2 , ..., X time , e N ) ′ .

Throughout this thesis, we will assume without loss of generality, that the state space of the finite-state Markov chain (X time ) is a set of unit vectors defined as above and that X 0 is given or its distribution π 0 is known.

Hidden Markov models

A Hidden Markov Model (HMM) is a pair of stochastic processes called the state process and the observation process, respectively. The state process is a hidden, that is an unobserved, homogeneous Markov chain modelling the environment, each state of the chain representing a specific regime of the environment. The observation process is a real valued function of the chain corrupted by a Gaussian noise (in discrete time) or is assumed to satisfy a stochastic differential equation (in continuous time). Such processes will be defined on a complete filtered probability space (Ω, F, (F time ), P ).

Discrete-time HMM

In this section we present a discrete time HMM that will be used to model daily clearness index sequences.

Consider a system whose states are described by a discrete-time homogeneous Markov chain (X h ) h=0,1,2,... , called the state process, with state space S = {e 1 , e 2 , ..., e N }.

For h = 1, 2, . . . , we will write

X h = ( X h , e 1 , X h , e 2 , . . . , X h , e N ) ′ .
Recall that X h , e i = 1 (X h =e i ) , i = 1, 2, . . . , N .

Remark 2.1. For i = 1, 2, . . . , N , we have

E( X h , e i ) = N j=1
e j , e i P (X h = e j ) = P (X h = e i ).

(2.5)

Let X h denote the σ-algebra generated by {X 0 , X 1 , . . . , X h } and let A = (a ji ) ∈ R N ×N denote the probability transition matrix of (X h ) h=1,2,... defined as

a ji = P (X h = e j |X h-1 = e i ), i, j = 1, 2, . . . , N.
Note that a ii = 1j =i a ji (i = 1, 2, . . . , N ).

The equation of the process (X h ) h=1,2,... , the so-called state equation, will be obtained from the following lemma.

Lemma 2.1. Let (V h ) h=1,2,... be a sequence defined by

V h X h -AX h-1 .
(2.6)

Then (V h ) h=1,2,... is a sequence of martingale increments.

Recall that a sequence of martingale increments is a random discrete series whose expectation with respect to the past is 0 (see Section 2.2).

Proof. The Markov property implies that

P (X h = e j |X h-1 ) = P (X h = e j |X h-1 ).
(2.7)

From (2.7) and Remark 2.1, we have

E(X h |X h-1 ) = E(X h |X h-1 ) = AX h-1 . (2.8) Thus, E(V h |X h-1 ) = E(X h -AX h-1 |X h-1 ) = AX h-1 -AX h-1 = 0.
(2.9)

From Lemma 2.1, the state process (X h ) h=1,2,... can be represented by the state equation: (2.10) where (V h ) h=1,2,... is a sequence of martingale increments. We suppose that the state process (X h ) is not observed directly, but rather observed through a function of the Markov chain (X h ), say (K h ) h=1,2,... . Definition 2.3. A discrete-time HMM is a pair of processes (X h , K h ) h=1,2,... determined by the following equations:

X h = AX h-1 + V h ,
X h = AX h-1 + V h , (2.11) K h = b(X h ) + α(X h )w h , (2.12)
where

(V h ) is sequence of martingales, (w h ) is a sequence of i.i.d. N (0, 1) random variables, w h is independent of σ{X 1 , X 2 , . . . , X h , K 1 , K 2 , . . . , K h } and b(X h ) = X h , b , α(X h ) = X h , α , with b = (b 1 , b 2 , . . . , b N ) ′ , α = (α 1 , α 2 , . . . , α N ) ′ .
The parameter vector of the model is defined as the vector:

θ = (a ji , 1 ≤ j = i ≤ N ; b 1 , b 2 , . . . , b N ; α 1 , α 2 , . . . , α N ),
We shall assume α i = 0 and thus without loss of genrality that α i > 0, 1 ≤ i ≤ N . We will assume that θ belongs to a compact set Θ ⊂ R N ×N +N and will denote by P θ a probability measure on (Ω, F) for which the process (X h , K h ) h=1,2,... satisfies equations (2.12) with parameter θ.

In practice the number of states N will be suggested by the user.

Filtrations, number of jumps, occupation time and level sums

The following notions will be useful for estimating the model parameters.

Definition 2.4. For h = 1, 2, . . . , let

G K h σ{X 1 , X 2 , . . . , X h , K 1 , K 2 , . . . , K h }, Y K h σ{K 1 , K 2 , . . . , K h }
be the σ-algebras generated by (X h , K h ) and (K h ), respectively.

These σ-algebras containing all the available information up to time h form increasing sequences, and thus are filtrations. G K h is called the filtration of complete data and Y K h the filtration of incomplete observation data. Note that Y K h ⊂ G K h ⊂ F h for all h. Estimating the model parameters requires the computation of the conditional expectations of the following quantities given the observation history Y K h .

Definition 2.5. Let (X h , K h ) h=1,2,... be a discrete-time HMM and let f (•) be any bounded function. Let us define:

1. the number of jumps of the Markov chain from e i to e j until time h by (2.13) 2. the occupation time of the Markov chain in state e i until time h by (2.14) 3. and the level sums of the observation process in state e i up to time h by

J ij h = h l=1 X l-1 , e i X l , e j ,
O i h = h l=1 X l , e i ,
T i h (f ) = h l=1 f (K l ) X l , e i .
(2.15)

Reference Probability Method of measure change

In the so-called Reference Probability Method [START_REF] Elliott | [END_REF][START_REF] Elliott | [END_REF], the main technique for obtaining ML (Maximum Likelihood) estimates of parameters is the change of measure which is a discrete time version of Girsanov theorem (see Theorem 2.5). To achieve such a mathematical objective, we are going to work in a "fictitious world" (Ω, (F h ), P ), called the reference probability space, with a reference probability measure P which is determined by a fix parameter vector θ 0 ∈ Θ and then we will relate these results in the "real world" (Ω, (F h ), P θ ) by using a back change of measure and by applying Bayes's rule.

The reference probability measure P , a convenient measure to work with, is defined such that under P the observations are i.i.d. random variables. Working under P , we will reformulate the initial estimate of the parameter vector θ in a "fictitious world" (Ω, (F time ), P ) where the well-known results for i.i.d. random variables can be applied. Then the results will be reinterpreted back to the real world (Ω, (F h ), P θ ) with the initial probability measure P θ (see [Elliott 2010, chap. 1, pg. 3-11] for more details).

Introduce a reference probability measure P from P θ such that under P :

(C3.1) (K h ) h=1,2,... is a sequence of N (0, 1) i.i.d. random variables which are inde- pendent of the Markov chain (X h ), (C3.2) (X h ) h=1,2,... is a Markov chain with transition matrix A so that E(X h |G K h-1 ) = AX h-1
, where E is the expectation under P . The existence of P with the characteristics above is described by Lemma 2.2 below.

Lemma 2.2. Let P θ be the measure determined by the parameter vector θ (Definition 2.3) and let P be a new probability measure such that

P P θ G K h = Λ K,θ h = h l=1 λ l ,
(2.16)

where λ l = X l ,α φ(K l ) φ K l -X l ,b X l ,α
, l = 1, 2, . . . , h , φ(•) is the N (0, 1) density function.

Then P satisfy conditions (C3.1), (C3.2) stated above.

Proof. Our proof starts by observing that the existence of P such that

(P /P θ )| G K h-1 = Λ K,θ h
follows from Kolmogorov's Extension Theorem. We first prove that Λ K,θ h is integrable and a martingale under P θ by (a) and (b) below.

(a) For l = 2, 3, . . . , h, we have

K l = X l , b + X t , α w l , so that w l = K l -X l ,b Xt,α with w l ∼ N (0, 1)). Moreover as G K l-1 ⊂ G K l-1 ∨ X l , it follows that E θ λ l |G K l-1 =E θ    E θ   X l , α φ (K l ) φ K l -X l ,b X l ,α G K l-1 ∨ X l   G K l-1    =E θ    +∞ -∞ X l , α φ (K l ) φ(w l ) φ(w l )dw l G K l-1    = +∞ -∞ φ(K l )dK l =1, (2.17) because w l is independent of G K h-1 . From this, E θ (Λ K,θ h ) = E θ h l=1 λ l = 1 < ∞. (b) By (2.17), E θ (Λ K,θ h |G K h-1 ) = Λ K,θ h-1 E θ λ h |G K h-1 = Λ K,θ h-1 . We now note that P (K h ≤ x|G K h-1 ) = E(1 (K h ≤x) |G K h-1 )
. By a version of Bayes theorem,

E(1 (K h ≤x) |G K h-1 ) = E θ (Λ K,θ h 1 (K h ≤x) |G K h-1 ) E θ (Λ K,θ h |G K h-1 ) = Λ K,θ h-1 Λ K,θ h-1 E θ λ h 1 (K h ≤x) |G K h-1 E θ λ h |G K h-1
.

(2.18) From (2.18) and (2.17), we have

P (K h ≤ x|G K h-1 ) =E θ   X h , α φ (K h ) φ K h -X h ,b X h ,α 1 (K h ≤x) G K h-1   = +∞ -∞ X h , α φ (K h ) φ(w h ) φ(w h )1 (K h ≤x) dw h = x -∞ φ(K h )dK h . Then P (K h ≤ x|G K h-1 ) = P (K h ≤ x)
), K h ∼ N (0, 1) under P . Let X h = σ{X 0 , X 1 , . . . , X h }, by considering similarly,

P (K h ≤ x|X h ) =E θ    E θ   X h , α φ (K h ) φ K h -X h ,b X h ,α 1 (K h ≤x) G K h   X h    =E θ    +∞ -∞ X h , α φ (K h ) φ(w h ) φ(w h )1 (K h ≤x) dw l X h    = +∞ -∞ X h , α φ (K h ) φ(w h ) φ(w h )1 (K h ≤x) dw h = x -∞ φ(K h )dK h ,
a quantity which is independent of X h . Condition (C3.1) follows.

Consider now condition (C3.2). We have

E(X h |G Y h-1 ) = E θ (Λ K,θ h X h |G Y h-1 ) E θ (Λ K,θ h |G Y h-1 ) =E θ   X h , α φ (K h ) φ K h -X h ,b X h ,α X h G K h-1   =E θ (X h |G K h-1 ) =AX h-1 .
So, under P , (X h ) remains a Markov chain with transition matrix A. This completes our proof.

Remark 2.2. Starting with the probability measure P we can recover the measure P θ by observing that

dP θ dP G K h = Λ K,θ h = h l=1 φ K l -X l ,b X l ,α X l , α φ (K l ) , (2.19) Note that Λ K,θ h Λ K,θ h = 1.

Normalized and unnormalized filters

Definition 2.6. For any

G K h -adapted sequence (H h ) (for instance H h ≡ J ij h , O i h or T i h (f )), define γ(H h ) E(Λ K,θ h H h |Y K h ), (2.20) π(H h ) E θ (H h |Y K h ), (2.21)
where Λ K,θ h is determined as in (2.19). The processes γ(H h ) and π(H h ) are called the unnormalized filter and the normalized filter of the process (H h ), respectively.

In the problem of parameter estimation, we have to compute the normalized filters π(H h ), for

H h ≡ J ij h , O i h or T i h (f ).
We first work on the reference probability space (Ω, (F h ), P ) to establish the equation of unnormalized filters γ(H h ). These results will be used to obtain the normalized filters π(H h ) in the real world (Ω, (F h ), P θ ). By using Bayes rule:

E θ (H h |Y K h ) = E(Λ K,θ h H h |Y K h ) E(Λ K,θ h |Y K h ) , (2.22)
which is equivalent to

π(H h ) = γ(H h ) γ(1) . (2.23)
On the other hand, it is not possible to compute directly the filter γ(H h ), for

H h ≡ J ij h , O i h or T i h (f ), but we can compute its associated process γ(H h X h ) = E(Λ K,θ h H h X h |Y K h ) and note that γ(H h ) =γ(H h X h , 1 ), because 1 = N i=1 X h , e i = X h , 1 , = γ(H h X h ), 1 .
Consequently, we can rewrite (2.23) as

π(H h ) = γ(H h X h ), 1 γ(X h ), 1 . (2.24)
Thus, for obtaining the normalized filters

π(H h ) (H h ≡ J ij h , O i h or T i h (f )
) in the real world (Ω, (F h ), P θ ), we must determine the equation of unnormalized filter processes γ(H h X h ). This will be computed in Chapter 3.

Some recalls on stochastic calculus

Ito product rule

For two semimartingales X t and Y t , Ito product rule gives

X t Y t = t 0 X s dY s + t 0 Y s dX s + [X, Y ] t , (2.25)
where

[X, Y ] t = lim n→∞ (in prob.) X 0 Y 0 + 0≤k<2 n (X t(k+1)2 -n -X tk2 -n ) × (Y t(k+1)2 -n -Y tk2 -n )
is the quadratic variation of X t and Y t .

Remark 2.3. If the process [X, Y ] t has a compensator denoted by X, Y t , it will be called the predictabe quadratic variation of X t and Y t . If the martingale part of the semimartingale is discontinuous, then

[X, Y ] t = X 0 Y 0 + 0<s≤t ∆X s ∆Y s .
Remark 2.4. If X t and Y t are Ito processes, then:

X t = X 0 + t 0 µ X s ds + t 0 σ X s dW s , Y t = X 0 + t 0 µ Y s ds + t 0 σ Y s dW s ,
where W t is a standard Brownian motion, and

[X, Y ] t = X 0 H 0 + t 0 σ X s σ Y s ds.
(2.26)

Ito formula

Let (ξ t ) t∈[0,T ] be a random process satisfying the stochastic differential equation

(SDE) dξ = a t dt + b t dW t , (2.27)
where (W t ) is a Wiener process, and the nonanticipative functions a t , b t are such that

P T 0 |a t |dt < ∞ = 1, (2.28)
P T 0 b 2 t dt < ∞ = 1. (2.29) Theorem 2.4. Let f (t, x) : [0, +∞) × R → R be a continuous function satisfying f t ∈ C 1 ([0, +∞)) and f x ∈ C 2 (R). Then the process f (t, ξ t ) satisfies the SDE df (t, ξ t ) = f ′ t (t, ξ t ) + f ′ x (t, ξ t )a t + 1 2 f ′′ xx (t, ξ t )b 2 t dt + f ′ x (t, ξ t )b t dW t .
(2.30) Theorem 2.4 which is Theorem 4.4 of [START_REF] Lipster | Statistics of random processes, parts i general theory[END_REF]]. The formula given by (2.30) was obtained by K. Itô and is called the Itô formula.

Girsanov theorem

Theorem 2.5. (Girsanov, [Elliott 2010, page 355]) Suppose that

(Y t ), t ∈ [0, T ], is a standard Brownian motion on a filtered space {Ω, F, (F t ), Q}. Let f : Ω×[0, T ] → R be a predictable process such that 1 T 0 |f t | 2 dt < ∞ a.s. E exp 1 2 T 0 f 2 t dt < ∞ 1 Note that if ft is bounded (|ft| ≤ L < ∞) then E exp 1 2 T 0 f 2 t dt ≤ exp T 2 L 2 < ∞ Let Λ t = exp t 0 f s dY s - 1 2 t 0 |f s | 2 ds , (2.31) and suppose that E Q [Λ t ] = 1 for all t ∈ [0, T ].
If P is the probability measure on {Ω, F} defined by

dP dQ = Λ T ,
then the process defined by

W t = Y t - t 0 f s ds.
is a standard Brownian motion on {Ω, F, F t , P }.

For a proof, refer to [Elliott 1982].

Remark 2.5. By Ito formula, the process (Λ t ) defined as in (2.31) satisfies the following equation:

Λ t = 1 + t 0 Λ s f s dY s , (2.32
)

and E Q [Λ t ] = 1 [Krishnamurthy 2002].
Remark 2.6. Let (Y t ) be a process determined on {Ω, F, F t , P } by the SDE:

dY t = f t dt + dW t , (2.33)
where (W t ) is a standard Brownian motion and the function f is defined as in Theorem 2.5. If Q is the measure defined by

dQ dP = (Λ t ) -1 = exp - t 0 f s dW s - 1 2 t 0 |f s | 2 ds , then (Y t
) is a standard Brownian motion under Q and we have, by Itô's formula,

(Λ t ) -1 = 1 - t 0 (Λ s ) -1 f s dW s .
(2.34)

Continuous-time homogeneous Markov chain

Consider a continuous-time Markov chain (X) t∈(0,T ) with state space S = {e 1 , e 2 , ..., e N }. Let p ji (t, h) be the transition probability from state e i at time t to state e j at time t + h:

p ji (t, h) = P (X t+h = e j |X t = e i ).
The Markov chain is said homogeneous if p ji (t, h) does not depend on t and only depends on h :

P (X t+h = e j |X t = e i ) p ji (h).
The time interval between state transition from state e i to state e j is a random variable H ji and note that its distribution is exponential with a given parameter λ ji [Gross 2012].

Therefore

P (H ji ≤ h) = 1 -exp(-λ ji h),
and

p ji (h) = 1 -exp(-λ ji h).
Let h = ∆t be a small increment of time so that:

p ji (∆t) = 1 -exp(-λ ji ∆t) ≈ 1 -(1 -λ ji ∆t) = λ ji ∆t.
Assume that p ii (∆t) is approximately an affine function of ∆t:

p ii (∆t) = 1 -λ i ∆t.
We can define the transition rates using

λ ji lim ∆t→0 p ji (∆t) ∆t , λ i lim ∆t→0 1 -p ii (∆t) ∆t . Note that p ii (∆t) + j =i p ji (∆t) = 1.
From this, we have the following relationship

λ i = lim ∆t→0 1 -p ii (∆t) ∆t = lim ∆t→0 j =i p ji (∆t) ∆t = j =i λ ji , ∀i.
The column sums of the probability transition matrix P (∆t) = (p ji (∆t)) is 1.

Define the transition rate matrix for (X t ) as

A ∆      -λ 1 λ 12 . . . λ 1N λ 21 -λ 2 . . . λ 2N . . . . . . . . . . . . λ N 1 λ N 2 . . . -λ N      .
It follows from the relations above that:

A ∆ = lim ∆t→0 P (∆t) -I ∆t ,
where I is the N × N unit matrix.

From this, it is seen that

P (∆t) ≈ ∆tA ∆ + I.
This approximation will be used to establish the transition probability matrix in ours approximate models later (Section 3.3.1 and Section 4.4.1).

Continuous-time HMM

Following [START_REF] Dembo | [END_REF][START_REF] James | [END_REF][START_REF] Elliott | [END_REF][START_REF] Elliott | [END_REF], we review now some of the standard definitions and notions of a continuous-time HMM based on a SDE with a noise variance taken to be one.

The state process of the model is a continuous-time homogeneous Markov chain (X t ) t∈[0,T ] with state space S = {e 1 , e 2 , ..., e N } and transition probabilities determined by:

P ji (∆) = P (X t+∆ = e j |X t = e i ) = λ ji ∆ + o(∆), if j = i, ∆ → 0, 1 + λ ii ∆ + o(∆), if i = j, ∆ → 0. (2.35)
where

λ ii = - j =i λ ji , i = 1, 2, . . . , N .
In this case, the transition matrix A ∆ = (λ ji ) ∈ R N ×N is a transition rate matrix (or infinitesimal generator, see Section 2.7). It is well-known that (X t ) now has a semimartingale representation [Elliott 2010, page 198]:

X t = X 0 + t 0 A ∆ X s ds + V t , (2.36)
where

(V t ) is an F t -martingale.
The process (X t ) is not observed directly, it is observed through a scalar process (Y t ) which is assumed to satisfy a stochastic differential equation.

Definition 2.7. A continuous-time HMM is a pair of processes (X t , Y t ) t∈[0,T ] such that

X t = X 0 + t 0 A ∆ X s ds + V t ,
(2.37)

dY t = h(X t , Y t , t)dt + dW t , (2.38)
where

(V t ) is an F t -martingale, (W t ) is a standard Brownian motion which is inde- pendent of (X t ) and h t = h(X t , Y t , t) is a predictable process such that T 0 |h t | 2 dt < ∞ a.s. h t is bounded.
Equation 2.38 is called the observation equation of the model.

Filtrations, number of jumps, occupation time and level sums

We next precise some notations in the continuous-time standard HMM which will be used in the problem of parameter estimation.

Definition 2.8. Let (X t , Y t ) t∈[0,T ] be a continuous-time standard HMM, define:

1. the filtration of complete data

G Y t σ{X s , Y s ; 0 ≤ s ≤ t}, (2.39)
2. the filtration of incomplete data

Y Y t σ{Y s ; 0 ≤ s ≤ t}, (2.40)
3. the number of jumps of the Markov chain from e i to e j until time t

J ij t = t 0 X s , e i dX s , e j ,
(2.41)

4. the occupation time of the Markov chain in state e i until time t (2.42) 5. the level sums of the observation process in state e i up to time t

O i t = t 0 X s , e i ds,
T i t (f ) = t 0 f (Y s ) X s , e i dY s , (2.43)
where f (•) is any bounded function.

Remark 2.7. We have

Y Y t ⊂ G Y t , for all t ∈ [0, T ].

Change of measure

Let θ be the parameter set of the model (X t , Y t ) (Definition 2.8) and write P θ for the measure determined by the set θ.

Again, in the computations of parameter estimation, we will also use the reference probability method with the change-of-measure technique similar to the discrete-time case.

For this, we first introduce a reference probability measure P from the initial measure P θ by putting:

dP dP θ G Y t = Λ Y,θ t = exp - t 0 h s dW s - 1 2 t 0 |h s | 2 ds .
(2.44) By Girsanov's theorem (see Theorem 2.5, Remark 2.6), (Y t ) is a standard Brownian motion under P .

We will work with a standard Brownian motion (Y t ) defined on the reference probability space (Ω, F, P ).

Then, the obtained results will be reinterpreted in the initial probability space (Ω, F, P θ ) by a back change of measure form P to P θ :

dP θ dP G Y t = Λ Y,θ t = exp t 0 h s dY s - 1 2 t 0 |h s | 2 ds . (2.45)
Specifically, the results obtained on (Ω, F, P ) for the unnormalized filters

γ(H t ) = E(Λ Y,θ t H t |Y Y t ), H t ≡ J ij t , O i t , T i t (f ), will be related to (Ω, F, P θ ) for ob- taining the normalized filters π(H t ) = E θ (H t |Y Y t
) by Bayes's rule: 1) .

E θ (H t |Y Y t ) = E(Λ Y,θ t H t |Y Y t ) E(Λ Y,θ t |Y Y t ) , (2.46) or π(H t ) = γ(H t ) γ(
(2.47)

Remark 2.8. As in the discrete-time HMM, we have

π(H t ) = γ(H t X t ), 1 γ(X t ), 1 .
(2.48)

Remark 2.9. From (2.45) and by applying Ito formula (see Section 2.6.2, Theorem 2.4) with

ξ t = t 0 h s dY s -1 2 t 0 |h s | 2 ds, f (t, ξ t ) = Λ Y,θ t = e ξt , we obtain Λ Y,θ t = 1 + t 0 Λ Y,θ s h s dY s , (2.49) and Λ Y,θ t = 1 - t 0 Λ Y,θ s h s dW s . (2.50) Note that Λ Y,θ t Λ Y,θ t = 1.

Parameter estimation

Because the state process is unobserved, the parameter estimation can be only based on incomplete data

Y • time (Y • time ≡ Y K h or Y • time ≡ Y Y t ).
In this case, the Maximum Likelihood Estimator cannot be easily handled [START_REF] James | [END_REF], Charalambous 2000[START_REF] Elliott | [END_REF]] so that several algorithms were proposed by statisticians and engineers, among the most successful ones is the EM algorithm.

We first describe the likelihood function, the pseudo log-likelihood used in the EM algorithm and the two main steps of this algorithm. We then summarize without proofs some relevant theorems on the stationary and the converging properties of the sequence of estimates.

Likelihood function

Let {P θ , θ ∈ Θ} be a family of probability measures and let P be a reference probability mesure defined as in Section 2.4, the likelihood function for computing an estimate θ ′ of parameter vector θ based on the available information Y • time is:

L(θ) = E dP θ dP Y • time ,
and the ML estimate θ M LE is defined by

θ M LE = arg max θ∈Θ L(θ).

Pseudo log-likelihood function

The EM algorithm provides an iterative numerical method which can be used to generate a sequence {θ (p) , p ≥ 0} for updating the ML estimates of parameter vector θ. To this end, we first consider the following log-likelihood function:

L(θ) = log E dP θ dP Y • time .
Due to the monotonicity of the log function, the maximization of L(θ) is equivalent to the maximization of L(θ).

Let P y θ denote the restriction of P θ to Y • time . It can be proved [START_REF] Dembo | [END_REF][START_REF] Lipster | Statistics of random processes, parts i general theory[END_REF]], that

dP y θ ′ dP y θ = E θ dP θ ′ dP θ Y • time .
Then

L(θ ′ ) -L(θ) = log dP y θ ′ dP y -log dP y θ dP y = log dP y θ ′ dP y θ , and
L(θ ′ ) -L(θ) = log E θ dP θ ′ dP θ Y • time .
Jensen inequality then yields

L(θ ′ ) -L(θ) ≥E θ log dP θ ′ dP θ Y • time ,
Therefore, for all θ, θ ′ ∈ Θ, we have obtained

L(θ ′ ) -L(θ) ≥ Q(θ ′ , θ),
with equality if and only if θ = θ ′ , where

Q(θ ′ , θ) E θ log Λ •,θθ ′ time |Y • time with Λ •,θθ ′ time dP θ ′ dP θ G • time .
(2.51)

The function Q(θ ′ , θ) as defined in (2.51) is called the pseudo log-likelihood function.

EM Algorithm

Starting with an initial parameter vector θ (0) , each iteration of EM algorithm consists in the two following steps: p) and compute the pseudo log-likelihood function

E-Step (Expectation Step) : From (2.51), set θ = θ ′ (
Q(θ ′ , θ ′ (p) ) = E θ ′(p) log Λ •,θ ′(p) θ ′ time |Y • time , M-Step (Maximization Step) : Find θ ′ (p+1) ∈ arg max θ ′ ∈Θ Q(θ ′ , θ ′ (p)
). Repeat from E-Step with p = p + 1, unless a stopping test is satisfied. The stationary and converging properties of the EM algorithm were evaluated by [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF][START_REF] Dembo | [END_REF]] via Theorem 2.6, Theorem 2.7 and Theorem 2.8 below.

Theorem 2.6. For every sequence {θ (p) , p = 0, 1, 2, . . . } generated by EM algorithm,

L(θ (p+1) ) ≥ L(θ (p) ) (2.52)
with equality if and only if

θ (p+1) = θ (p) .
Theorem 2.6 is Theorem 1 in [START_REF] Dembo | [END_REF]] and the infinite dimensional analogue of Theorem 1 in [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF]].

Theorem 2.7. Suppose that {θ (p) , p = 0, 1, 2, . . . } is an instance generated by EM algorithm such that (1) the sequence L(θ (p) ) is bounded, and

(2) Q(θ (p+1) |(θ (p) )) -Q(θ (p) |(θ (p) )) ≥ λ(θ (p+1) -θ (p)
) 2 for some λ > 0 and all p.

Then the sequence {θ (p) , p = 0, 1, 2, . . . } converges to some θ * in the closure of Θ.

The proof can be found in [Dempster 1977, Theorem 2].

Theorem 2.8. Assume that

(1) Θ θ 0 = {θ ∈ Θ : L(θ) ≥ L(θ 0 )} is compact for any θ 0 ∈ Θ.
(2) L(θ) is continuous in Θ and differentiable in the interior of Θ.

(3) Q(θ ′ , θ) is continuous with respect to θ and θ ′ .

(4) All the EM instances {θ (p) , p = 0, 1, 2, . . . } are in the interior of Θ.

Then, all the limit points of any instance {θ (p) , p = 0, 1, 2, . . . } of the EM algorithm (and there exists at least one such limit point) are stationary points of L(•), having the same value L * (of L(•)). Furthermore, L(θ (p) ) converges monotonically to L * . Theorem 2.8 was proved by [Dembo 1986, Theorem 2].

Chapter 3

Stochastic models for clearness index processes

Résumé

Dans ce chapitre, nous proposons un modèle de type HMM à temps discret, noté DTM-K, pour une suite d'indices journaliers de clarté (K h ) h=1,2,... , ainsi qu'un modèle de type HMM à temps continu, noté CTM-k, pour un processus d'indice de clarté (k t ) t∈[0,T ] sur un intervalle de temps [0, T ]. Les estimations numériques des paramètres de CTM-k nécessite de construire un modèle approché à temps discret, noté DTAM-k, obtenu en discrétisant le temps par une partition de [0, T ].

Pour chacun de ces modèles, nous définissons un processus des états, un processus des observations, et un vecteur de paramètres. Nous détaillons les équations de filtrages pour obtenir les formules de mise à jour des estimations dans les itérations de l'algorithme EM, en particulier pour estimer la variance du bruit.

Les applications numériques ont été faites en utilisant des données provenant de mesures effectuées à La Réunion pour DTM-K et en Guadeloupe pour CTM-k et DTAM-k.

Les résultats de ce chapitre ont été présentés à la conférence [Tran 2013] (voir aussi [Tran a]).

Abstract

In this chapter, we propose a discrete-time HMM model, denoted DTM-K, for a sequence of daily clearness index (K h ) h=1,2,... , and also a continuous-time HMM model, denoted CTM-k, for a process (k t ) t∈[0,T ] of clearness index on a time interval [0, T ]. Numerical estimations of CTM-k parameters requires to build a discrete-time approximating model, denoted by DTAM-k, obtained by discretizing time, using a partition of [0, T ].

For each of these models, we define a state process, an observation process, and a parameter vector. We detail the filtering equations in order to get estimation update formulas used in the iterations of EM algorithm, in particular for estimating the noise variance.

Experiments for DTM-K (resp. for CTM-k and DTAM-k) are done using data coming from measurements performed in La Réunion island (resp. in Guadeloupe island).

The results of this chapter have been presented in [Tran 2013] conference (also see [Tran a]). 3.3).

Modelling a daily clearness index sequence

The empirical distribution of a daily Clearness Index Sequence (CISs) during a period (for instance during October 2010 Figure 3.1a) suggests that the daily CIS distribution could be a Gaussian mixture, each Gaussian component corresponding, may be, to some specific meteorological regime. This has lead us to modelize the dynamic of the sequence by a discrete-time HMM where 1. the unobserved state process is a Markov chain representing the dynamic of regimes, each daily index belonging to a regime, several daily indices belonging eventually to a same regime, 2. the observed process is such that, given (or within) regime i, the various observed daily clearness indices are outcomes of a Gaussian distribution whose mean b i and standard deviation α i depend on regime i, i = 1, 2, . . . , N .

Actually, each regime corresponds to a Gaussian component of the suggested Gaussian mixture, and in terms of probabilistic classification, each regime corresponds to a (Gaussian) class. The advantage of considering a HMM is that it provides a parametric description of the random dynamic of the regimes, which is not the case in a classification setting.

In this section we propose a discrete-time HMM (X h , K h ) h=1,2,... to model a daily clearness index sequence in random environment. We first describe the state process, the observation process and the parameters of the model denoted DTM-K. Then we will detail the pseudo log-likelihood function, the filtering equations and the computations used in the EM algorithm. Experiments with real data will be presented in Section 3.4.

State process

The random dynamic of meteorological regimes will be modellized by a state process, an unobserved discrete-time, finite-state homogeneous Markov chain (X h ) h=0,1,2,... (Figure 3.2). The state space of (X h ) is the set of unit vectors:

S = {e 1 , e 2 , . . . , e N }, e i = (0, . . . , 0, 1 i-th , 0, . . . , 0) ′ ∈ R N .
Recall that we can write X h = ( X h , e 1 , X h , e 2 , . . . , X h , e N ) ′ and the process (X h ) h=1,2,... can be represented by the following equation, called the state equation,

X h = AX h-1 + V h , (3.1)
where (V h ) h=1,2,... is a sequence of martingale increments.

Observation Process and model parameters

The random values of a daily CIS (K h ) are modelled by the so-called observation process as follows. In regime i, that is when the Markov chain is in state e i , (i = 1, 2, . . . , N ) the daily clearness index K h will be considered as an outcome of a Gaussian distribution N (b i , α 2 i ) depending on regime i. In other words:

1 (X h =e i ) K h = 1 (X h =e i ) (b i + α i w h ), h = 1, 2, . . . , (3.2)
where w h are N (0, 1) i.i.d. random variables and b i , α i are parameters to be estimated.

The DM-K model of daily CIS under the random effects of meteorological events will be the HMM (X h , K h ) h=1,2,... , where the state process (X h ) is determined by (3.1) and the observation process (K h ) is defined as

K h = N i=1 1 (X h =e i ) K h = N i=1 1 (X i =e i ) (b i + α i w h ) , h = 1, 2, . . . . (3.3) If b = (b 1 , b 2 , . . . , b N ) ′ and α = (α 1 , α 2 , .
. . , α N ) ′ , then we have equivalently

K h = X h , b + X h , α w h , h = 1, 2, . . . . (3.4)
The parameter vector for the proposed model DM-K is

θ = (a ji , 1 ≤ j = i ≤ N ; b 1 , b 2 , . . . , b N ; α 1 , α 2 , . . . , α N ),
where the number of states N will be suggested by the user.

Parameter estimation

As mentioned in Chapter 2, our aim is to determine a new parameter vector

θ ′ = {(a ′ ji ), 1 ≤ j = i ≤ N ; b ′ , α ′ } which maximizes the pseudo log-likelihood function at each iteration of EM algorithm, where b ′ = (b ′ 1 , b ′ 2 , . . . , b ′ N ) ′ , α ′ = (α ′ 1 , α ′ 2 , . . . , α ′ N ) ′ .
Updating the transition probabilities a ji is as follows ( [Elliott 1995, chap. 3, pg. 68-70] and [START_REF] James | [END_REF]]):

a ′ ji = π(J ji h ) π(O i h ) 1 ≤ j = i ≤ N, (3.5)
where π(J ji h ) and π(O i h ) are the normalized filters (see Section 2.5.3) of the number of jumps and the occupation time, respectively.

Note that in (3.5), π is built using θ (corresponding to θ ′(p) in p-th iteration of the EM algorithm), and not θ ′ (corresponding to θ ′(p+1) ).

We now consider the update from b and α to b ′ and α ′ , respectively.

Pseudo log-likelihood function

Given the filtration Y K h of observed data at day h, consider the conditional expectation of the log-likelihood function:

Q(θ ′ , θ) = E θ log Λ K,θθ ′ h |Y K h , (3.6) where Λ K,θθ ′ h = dP θ ′ dP θ G K h . We first have Λ K,θθ ′ h = h l=1 X l , α φ K l -X l ,b ′ X l ,α ′ X l , α ′ φ K l -X l ,b X l ,α , (3.7)
where φ(•) denotes the N (0, 1) density. From (3.6) and (3.7), we get:

Q(θ ′ , θ) = E θ h l=1 log 1 X l , α ′ - 1 2 K l -X l , b ′ X l , α ′ 2 Y K h + R(θ, Y K h ), (3.8)
where the function R(θ, Y K h ) does not depend on θ ′ . p) and rewrite (3.8) as

3.1.3.2 Computations in EM algorithm E-Step (Expectation Step) : Set θ = θ ′ (
Q(θ ′ , θ ′ (p) ) = N i=1 E θ ′ (p) h l=1 X l , e i log 1 α ′ i - 1 2 X l , e i α ′2 i (K l -b ′ i ) 2 Y K h + R(θ ′(p) , Y K h ) = N i=1 π(O i h ) log 1 α ′ i - 1 2α ′2 i π(T i h (K 2 h )) -2b ′ i π(T i h (K h )) + b ′2 i π(O i h ) + R(θ ′(p) , Y K h ), where
O i h = h l=1 X l , e i , (3.9)
T i h (K h ) = h l=1 K h X l , e i , (3.10) T i h (K 2 h ) = h l=1 K 2 h X l , e i , (3.11) and π(H h ) = E θ ′(p) (H h |Y K h ) for H h ≡ O i h , T i h (K h ) or T i h (K 2 h ). M-Step (Maximization Step) : Let us find now θ ′ (p+1) ∈ arg max θ ′ ∈Θ Q(θ ′ , θ ′ (p) ). 1. Taking derivative of Q(θ ′ , θ ′ (p) ) with respect to b i , i = 1, 2, . . . , N , we obtain ∂ ∂b ′ i Q(θ ′ , θ ′ (p) ) = - 1 2α ′2 i -2π(T i h (K h )) + 2b ′ i π(O i h ) . (3.12) Now ∂ ∂b ′ i Q(θ ′ , θ ′ (p) ) = 0 yields b ′ i = π(T i h (K h )) π(O i h ) . (3.13) 2. Similarly, for i = 1, 2, . . . , N , ∂ ∂α ′ i Q(θ ′ , θ ′ (p) ) = 0 yields α ′2 i = 1 π(O i h ) π(T i h (K 2 h )) -2b ′ i π(T i h (K h )) + b ′2 i π(O i h ) .
(3.14)

Updating parameter

For H h ≡ O i h , T i h (K h ) or T i h (K 2 h ), as π(H h ) = γ(H h X h ),1 γ(X h
),1 , (3.5), (3.13), (3.14) become:

a ′ ji = γ(X h J ij h ), 1 γ(X h O i h ), 1 , b ′ i = γ(X h T i h (K h )), 1 γ(X h O i h ), 1 , 1 ≤ j = i ≤ N, (3.15) α ′2 i = γ(X h O i h ), 1 -1 × γ(X h T i h (K 2 h )), 1 -2b ′ i γ(X h T i h (K h )), 1 + b ′2 i γ(X h O i h ), 1 . (3.16)
It remains to establish some recursive equations to compute γ(

X h J ji h ), γ(X h O i h ), γ(X h T i h (K h )) and γ(X h T i h (K 2 h ))
. This is done in the next section.

Filtering equations

Definition 3.1. Let Λ K,θ h be the process determined as in (2.19), h = 1, 2, . . . . Define

q h = γ(X h ) = E(Λ K,θ h X h |Y K h ).
(3.17)

The process (q h ) is called the state filter process of model (X h , K h ).

In oder to obtain a recursive equation for (q h ), h = 1, 2, . . . , introduce the following quantities:

B (i) h φ K h -b i α i α i φ (K h ) , i = 1, 2, . . . , N, B h diag(B (1) h , B (2)
h , . . . , B

(N ) h ).

Lemma 3.1. The recursive equation of state process q h is given by

q h = B h Aq h-1 , h = 1, 2, . . . (3.18)
with initial condition q 0 = π 0 .

Proof. By definition of (q h ) h=1,2,... we have

q h =E(Λ K,θ h X h |Y K h ) =E Λ K,θ h-1 φ K h -X h ,b X h α X h , α φ (K h ) X h Y K h = N i=1 E Λ K,θ h-1 X h , e i |Y K h-1 B (i) h e i = N i=1 E E Λ K,θ h-1 X h , e i |G K h-1 Y K h-1 B (i) h e i , Noticing that Y K h-1 ⊂ G K h-1 , we get q h = N i=1 E Λ K,θ h-1 E X h , e i |G K h-1 Y K h-1 B (i) h e i = N i=1 E Λ K,θ h-1 AX h-1 , e i |K h-1 B (i) h e i = N i=1 Aq h-1 , e i B (i) h e i =B h Aq h-1 .
Using the state filter process (q h ) computed in Lemma 3.1 above, the recursive equation for the number of jumps, for the occupation time and for the level sum are now stated in Theorem 3.1 and Theorem 3.2 below.

Theorem 3.1. Let J ij h = h l=1
X l-1 , e i X h , e j be the number of jumps from state e i to state e j up to time h. The process γ(J ij h X h ) is computed by the following recursive equation:

γ(J ij h X h ) = B h Aγ(J ij h-1 X h-1 ) + q h-1 , e i B h Ae i , e j e j , h = 1, 2, . . . , (3.19)
with initial condition γ(J ij 0 X 0 ) = 0.

Proof. We have

γ(J ij h X h ) =E(Λ K,θ h J ij h X h |Y K h ) =E Λ K,θ h-1 φ K h -X h ,b X h ,α X h , α φ (K h ) J ij h-1 + X h-1 , e i X h , e j X h Y K h = N l=1 E Λ K,θ h-1 X h , e l J ij h-1 + X h-1 , e i X h , e j |Y K h-1 B (l) h e l = N l=1 E Λ K,θ h-1 X h , e l J ij h-1 |Y K h-1 B (l) h e l + E Λ K,θ h-1 X h , e j X h-1 , e i |Y K h-1 B (j) h e j . (3.20) As Y K h ⊂ G K h , conditional expectation properties imply that E Λ K,θ h-1 X h , e l J ij h-1 |Y K h-1 =E E Λ K,θ h-1 X h , e l J ij h-1 |G K h-1 Y K h-1 =E Λ K,θ h-1 AX h-1 , e l J ij h-1 |Y K h-1 = Aγ(J ij h-1 X h ), e l . (3.21)
Similarly,

E Λ K,θ h-1 X h , e l X h-1 , e i |Y K h-1 =E E Λ K,θ h-1 X h , e j X h-1 , e i |G K h-1 Y K h-1 =E Λ K,θ h-1 AX h-1 , e j X h-1 , e i |Y K h-1 =a ji q h-1 , e i . ( 3
h = B h Ae i , e j , we get γ(J ij h X h ) = N l=1 Aγ(J ij h-1 X h ), e l B (l) h e l + a ji q h-1 , e i B (j) h e j =B h Aγ(J ij h-1 X h-1 ) + q h-1 , e i B h Ae i , e j e j , h = 1, 2, . . . . (3.23) Theorem 3.2. Let T i h (f ) = h l=1 f (K l ) X l-1
, e i be the level sum for state e i , where

f (•) is any bounded function. The process γ(T i h (f )X h ) is determined by: γ(T i h (f )X h ) = B h Aγ(T i h-1 (f )X h-1 ) + f (K h ) q h-1 , e i B h Ae i , (3.24)
with initial condition γ(T i 0 (f )X 0 ) = 0.

Proof. Starting with the definition of γ(T i h (f )X h ) and noticing that Y K h ⊂ G K h , similar arguments as those used in Lemma 3.1 and Theorem 3.1 give :

γ(T i h (f )X h ) =E(Λ K,θ h T i h (f )X h |Y K h ) =E Λ K,θ h-1 φ K h -X h ,b X h ,α X h , α φ (K h ) T i h-1 (f ) + f (K h ) X h-1 , e i X h Y K h = N l=1 E Λ K,θ h-1 X h , e l T i h-1 (f ) Y K h-1 B (l) h e l + N l=1 E Λ K,θ h-1 X h , e l X h-1 , e i Y K h-1 f (K h )B (l) h e l = N l=1 E E Λ K,θ h-1 X h , e l T i h-1 (f ) G K h-1 Y K h-1 B (l) h e l + N l=1 E E Λ K,θ h-1 X h , e l X h-1 , e i G K h-1 Y K h-1 f (K h )B (l) h e l .
It follows that

γ(T i h (f )X h ) = N l=1 E Λ K,θ h-1 AX h-1 , e l T i h-1 (f ) Y K h-1 B (l) h e l + N l=1 E Λ K,θ h-1 X h-1 , e i Y K h-1 a li f (K h )B (l) h e l = N l=1 Aγ(T i h-1 (f )X h-1 ), e l B (l) h e l + N l=1 q h-1 , e i a li f (K h )B (l) h e l =B h Aγ(T i h-1 (f )X h-1 ) + f (K h ) q h-1 , e i B h Ae i .
Now, we apply Theorem 3.2 to get the recursive equations for γ

(O i h X h ), γ(T i h (K h )X h ) and γ(T i h (K 2 h )X h ). Taking f (K h ) ≡ 1 in (3.24), we get: Corollary 3.1. γ(O i h X h ) = B h Aγ(O i h-1 (f )X h-1 ) + q h-1 , e i B h Ae i , (3.25)
with initial condition γ(O i 0 X 0 ) = 0.

Applying Theorem 3.2 with f (K h ) ≡ K h and then with f (K h ) ≡ K 2 h , we obtain the recursive equations for γ(T i h (K h )X h ) and γ(T i h (K 2 h )X h ):

Corollary 3.2. γ(T i h (K h )X h ) =B h Aγ(T i h-1 (K h )X h-1 ) + K h q h-1 , e i B h Ae i , (3.26) γ(T i h (K 2 h )X h ) =B h Aγ(T i h-1 (K 2 h )X h-1 ) + K 2 h q h-1 , e i B h Ae i , (3.27) with the initial conditions γ(T i 0 (K h )X 0 ) = γ(T i 0 (K 2 h )X 0 ) = 0.
3.2 Modelling a clearness index process on a time interval

CTM-k model

We now describe CTM-k, a continuous-time HMM (X t , k t ) t∈[0,T ] for modelling a stochastic process of clearness index (k t ) t∈[0,T ] during a continuous-time interval of a solar day. We define the state process (X t ) t∈[0,T ] as a continuous-time homogeneous Markov chain (Section 2.8): the state space is again S = {e 1 , e 2 , . . . , e N }, where e i denotes the unit N -vector with 1 at the i-th position, and the transition rate matrix

A ∆ = (λ ji ) ∈ R N ×N with λ ii = - j =i λ ji for i = 1, 2, . . . , N .
We define the observation process (k t ), k t representing the random clearness index at time t ∈ [0, T ], as a solution of the following SDE:

dk t = X t , c dt + X t , β dW t , 0 ≤ t ≤ T, (3.28)
where (W t ) is a standard Brownian motion, and

c =(c 1 , c 2 , . . . , c N ) ′ , β =(β 1 , β 2 , . . . , β N ) ′ .
Our model CTM-k (X t , k t ) depends on the parameter vector:

θ = (λ ji , 1 ≤ j = i ≤ N ; c 1 , c 2 , . . . , c N ; β 1 , β 2 , . . . , β N ).
Let G k t = σ{X s , k s ; 0 ≤ s ≤ t} and Y k t = σ{k s ; 0 ≤ s ≤ t} denote the filtrations of the complete data and the incomplete data, respectively.

Change of measure

For each parameter vector θ ∈ Θ, first consider the change of probability measure from the initial measure P θ in the real world (Ω, (F t ), P θ ) to the reference probability measure P by putting We now consider the reverse change, from P to P θ , by putting:

dP dP θ G k t = Λ k,θ t = exp - t 0 X s , c X s , β dW s - 1 2 t 0 X s , c X s , β 2 ds . ( 3
dP θ dP G k t = Λ k,θ t = exp t 0 X s , c X s , β d k t X t , β - 1 2 t 0 X s , c X s , β 2 ds . (3.31)
Then, Ito formula yields

Λ k,θ t = 1 + t 0 Λ θ s X s , c X s , β d k s X s , β , (3.32)
and we have E(Λ k,θ t ) = 1. With this back change, Girsanov theorem implies that under P θ , Wt is a standard Brownian motion if we define

d Wt = d k t X t , β - X t , c X t , β dt.
(3.33) Thus, under P θ , (k t ) satisfies the following equation equivalent to (3.28), due to Ito formula and (3.33):

dk t = X t , c dt + X t , β d Wt , (3.34)
where ( Wt ) is a standard Brownian motion.

That is, under P θ , (k t ) satisfies the real world dynamics. However, as P is a more convenient measure, we will work on the reference probability space (Ω, F t , P ) to compute the unnormalized filters γ(H

t ) = E(Λ k,θ t H t |Y k t )
, where (H t ) is a G k tadopted scalar process (for example H t ≡ J ij t , O i t , . . . ). These results will be related to (and used in) the real world (Ω, (F t ), P θ ) by using a version of conditional Bayes Theorem:

E θ (H t |Y k t ) = E(Λ k,θ t H t |Y k t ) E(Λ k,θ t |Y k t )
.

(3.35)

A complete proof for (3.35) can be found in [Elliott 1995, Theorem 3.2]. We here consider in cases where

E(Λ k,θ t |Y k t ) > 0. Suppose B is any set in G = {ω : E(Λ k,θ t |Y k t ) > 0)} ⊂ Y k t , we must show B E(Λ k,θ t H t |Y k t ) E(Λ k,θ t |Y k t ) dP θ = B E θ (H t |Y k t )dP θ .
Proof. By noting E E ξ|Y k t = E(ξ) for any ξ, we have:

B E(Λ k,θ t H t |Y k t ) E(Λ k,θ t |Y k t ) dP θ =E θ 1 B E(Λ k,θ t H t |Y k t ) E(Λ k,θ t |Y k t ) =E 1 B Λ k,θ t E(Λ k,θ t H t |Y k t ) E(Λ k,θ t |Y k t ) =E E 1 B Λ k,θ t E(Λ k,θ t H t |Y k t ) E(Λ k,θ t |Y k t ) Y k t =E 1 B E[Λ k,θ t |Y k t ] E(Λ k,θ t H t |Y k t ) E(Λ k,θ t |Y k t ) =E 1 B E(Λ k,θ t H t |Y k t ) =E 1 B Λ k,θ t H t = B Λ k,θ t H t dP = B H t dP θ = B E θ (H t |Y k t )dP θ .
Remark 3.1. By Ito formula, the process in (3.29) can be determined by

Λ k,θ t = 1 - t 0 Λ k,θ s X s , c X s , β dW s , (3.36)
and we have Λ k,θ t Λ k,θ t = 1.

Parameter estimation

At each iteration of EM algorithm, the new transition probabilities will be updated by [Elliott 1995, chap. 3, pg. 68-70] and [START_REF] James | [END_REF]]:

λ ′ ji = π(J ji t ) π(O i t ) , 1 ≤ j = i ≤ N, (3.37)
where π(J ji t ) and π(O i t ) are the normalized filter (see Section 2.8.2). Because measures corresponding to Wiener processes with different variances are not absolutely continuous [START_REF] Lipster | Statistics of random processes, parts i general theory[END_REF][START_REF] James | [END_REF], we cannot obtain the ML estimates of parameters β 1 , β 2 , . . . , β N in continuous time. However, we get around this problem using discretization by introducing a Discrete Time Approximate Model (DTAM) in which we can use densities with respect to the Lebesgue measure and compute Radon-Nikodym derivatives of observation processes with different noisy variances.

We now use the EM algorithm mentioned in Section 2.9 to update the ML estimates of the parameter c = (c 1 , c 2 , . . . , c N ) ′ .

Expectation step

We consider the pseudo log-likelihood function defined as:

Q(θ ′ , θ) = E θ log Λ k,θθ ′ t |Y k t ,
where

Λ k,θθ ′ t = dP θ ′ dP θ G k t = exp t 0 X s , c ′ X s , β ′ 2 dk s - 1 2 t 0 X s , c ′ X s , β ′ 2 ds × exp t 0 X s , c X s , β 2 dk s - 1 2 t 0 X s , c X s , β 2 ds -1 = exp t 0 X s , c ′ X s , β ′ 2 - X s , c X s , β 2 dk s - 1 2 t 0 X s , c ′ X s , β ′ 2 - X s , c X s , β 2 ds .
Hence

Q(θ ′ , θ) =E θ t 0 X s , c ′ X s , β ′ 2 - X s , c X s , β 2 dk s - 1 2 t 0 X s , c ′ X s , β ′ 2 - X s , c X s , β 2 ds Y k t (3.38) =E θ N i=1 1 β ′2 i c ′ i t 0 X s , e i dk s - 1 2 c ′2 i t 0 X s , e i ds K t + R(θ, Y k t ), (3.39)
where R(θ, Y k t ) does not depend on θ ′ . Then, by setting θ = θ ′(p) at the p-th iteration of EM algorithm, we get

Q(θ ′ , θ ′(p) ) =E θ ′(p) N i=1 1 β ′2 i c ′ i t 0 X s , e i dk s - 1 2 c ′2 i t 0 X s , e i ds Y k t + R(θ ′(p) , Y k t ). (3.40) Using the notation π(H t ) = E θ ′(p) (H t |Y k t ), for H t ≡ O i t or T i t with O i t = t 0 X s , e i ds, T i t = t 0 X s , e i dk s , i = 1, 2, . . . , N, (3.41)
we rewrite (3.40) as:

Q(θ ′ , θ ′(p) ) = N i=1 1 β ′2 i c ′ i π(T i t ) - 1 2 c ′2 i π(O i t ) + R(θ ′(p) , Y k t ).
(3.42)

3.2.3.2 Maximization step From (3.42), ∂ ∂c ′ i Q(θ ′ , θ ′(p) ) = 0 becomes c ′ i = π(T i t ) π(O i t ) , i = 1, 2, . . . , N. (3.43) Observing that π(H t ) = γ(Ht) γ(1) = γ(XtHt),1 γ(Xt),1
(for

H t ≡ J ′ji t , O i t or T i t ), update from λ ji to λ ′
ji in (3.37) and update from c i to c ′ i in (3.43) up to time T are now given by

λ ′ ji = γ(X T J ji T ), 1 γ(X T O i T ), 1 , 1 ≤ j = i ≤ N, (3.44) c ′ i = γ(X T T i T ), 1 γ(X T O i T ), 1 , 1 ≤ i ≤ N, (3.45)
where the filter processes γ(X T O i T ), γ(X T T i T ) and γ(X T J ji T ) are obtained from Theorem 3.3 and Theorem 3.4 in the next section.

Filtering equations

Let the state process q t be defined as

q t = γ(X t ) = E(X t Λ k,θ t |Y k t )
and let

C = diag c 1 β 1 , c 2 β 2 , . . . , c N β N .
Theorem 3.3. Let O i t and T i t be respectively the occupation time and the level sum defined in (3.41). The recursive equations of filtering processes γ(X t O i t ), γ(X t T i t ) and q t are given by

q t =π 0 + t 0 A ∆ q s ds + t 0 Cq s dk s X t , β , (3.46) γ(X t O i t ) = t 0 A ∆ γ(X s O i s )ds + t 0 Cγ(X s O i s ) dk s X s , β + t 0 q s , e i e i ds, (3.47) γ(X t T i t ) = t 0 A ∆ γ(X s T i s )ds + t 0 Cγ(X s T i s ) dk s X s , β + t 0 c i q s , e i e i ds + t 0
q s , e i e i dk s .

(3.48)

Proof. Firstly, remark that the filtering equations of Theorem 3.3 which involvethe process X t are Y k -measurable, Y k = σ{k s , 0 ≤ s ≤ t}, because the observation process (k t ) is determined by a SDE depending on X t :

dk t = X t , c dt + X t , β dW t >
Consider a general scalar process defined as:

H t = H 0 + t 0 ρ s ds + t 0 δ s dW s , (3.49)
where ρ t and δ t are G t -predictable and square-integrable processes. Ito product rule for semimartingales yields:

X t H t = X 0 H 0 + t 0 X s ρ s ds + t 0 X s δ s dW s + t 0 H s A ∆ X s ds + t 0 H s dV s , (3.50) in which it should be noticed that [X, H] t = X 0 H 0 because ∆X s ∆H s = 0 a.s.
Again Ito product rule applied to X t H t and Λ k,θ t yields:

X t H t Λ k,θ t =X 0 H 0 + t 0 X s ρ s Λ k,θ s ds + t 0 X s δ s Λ k,θ s dW s + t 0 H s A ∆ X s Λ k,θ s ds + t 0 H s Λ k,θ s dV s + t 0 Λ k,θ s X s , c X s , β X s H s dk s X s , β + [XH, Λ k,θ ] t , (3.51)
where

[XH, Λ k,θ ] t = t 0 Λ k,θ s X s , c X s , β X s δ s ds. (3.52)
As X s , e i X s = X s , e i e i , we have

t 0 Λ k,θ s X s , c X s , β X s H s dk s X s , β = N i=1 t 0 Λ k,θ s X s H s , e i c i β i dk s β i e i , (3.53) t 0 X s δ s Λ k,θ s dB s + t 0 Λ k,θ s X s , c X s , β X s δ s ds = t 0 X s δ s Λ k,θ s dk s X s , β = N i=1 t 0 Λ k,θ s δ s X s , e i dk s β i e i . ( 3
X t H t Λ k,θ t =X 0 H 0 + t 0 X s ρ s Λ k,θ s ds + N i=1 t 0 Λ k,θ s δ s X s , e i dk s β i e i + t 0 H s A ∆ X s Λ k,θ s ds + t 0 H s Λ k,θ s dV s + N i=1 t 0 Λ k,θ s X s H s , e i c i β i dk s β i e i .
(3.55)

Appliying conditional expectation E given Y k t to both sides of (3.55) and noticing that E(dV

s |Y k t ) = 0 ([Elliott 2010, Lemma 2.1, page 198]), it follows that γ(X t H t ) =γ(X 0 H 0 ) + t 0 γ(X s ρ s )ds + N i=1 t 0 γ(δ s X s ), e i dk s β i e i + t 0 A ∆ γ(X s H s )ds + N i=1 t 0 γ(X s H s ), e i c i β i dk s β i e i =γ(X 0 H 0 ) + t 0 γ(X s ρ s )ds + t 0 γ(δ s X s ) dk s X s , β + t 0 A ∆ γ(X s H s )ds + t 0 Cγ(X s H s ) dk s X s , β . (3.56)
We now apply (3.56) to obtain the recursive equations of q t , O t i and T i t :

1. Taking H t = H 0 = 1, ρ s = δ s = 0, we obtain the recursive equation for the state process q t = γ(X t ) = E(X t Λ θ t |Y k t ):

q t = π 0 + t 0 A ∆ q s ds + t 0 Cq s dk s X s , β . (3.57) 2. For the occupation time H t ≡ O i t = t 0 X s , e i ds, i = 1, 2, . . . , N , we take H t = O i t , H 0 = 0, ρ s = X s , e i , δ s = 0 and X t ρ s = X s X s , e i = X s , e i e i . Applying (3.56), we obtain γ(X t O i t ) = t 0 A ∆ γ(X s O i s )ds+ t 0 Cγ(X s O i s ) dk s X s , β + t 0
q s , e i e i ds. (3.58)

3. For the level sum T i t = t 0 X s , e i dk s , i = 1, 2, . . . , N , we first observe that

T i t = t 0 c i X s , e i ds + t 0 β i X s , e i dB s . (3.59)
Then, taking H t = T i t , H 0 = 0, ρ s = c i X s , e i , δ s = β i X s , e i and applying (3.56), we get

γ(X t T i t ) = t 0 A ∆ γ(X s T i s )ds + t 0 Cγ(X s T i s ) dk s X s , β + t 0 c i q s , e i e i ds + t 0 q s , e i e i dk s . (3.60) Theorem 3.4. Let J ij t = t 0 X
s , e i dX s , e j be the number of jumps from state e i to state e j , 1

≤ i = j ≤ N , t ∈ [0, T ]. The filtering process γ(J ij t ) is obtained from the following recursive equation γ(X t J ij t ) = t 0 A ∆ γ(X s J ij s )ds+ t 0 Cγ(X s J ij s ) dk s X s , β + t 0
q s , e j λ ji e i ds, (3.61)

where C = diag c 1 β 1 , c 2 β 2 , . . . , c N β N .
Proof. Using dX s = AX s ds + dV s , we first decompose

J ij t = t 0 X s , e i dX s , e j = t 0 X s , e i e ′ j dX s = t 0 X s , e i a ij ds + t 0 X s , e i e ′ j dV s . (3.62)
Applying Itô product rule and noticing that X s X s , e i = X s , e i e i , we have

X t J ij t = t 0 X s , e i a ji dse i + t 0 X s , e i e ′ j dV s e i + t 0 AX s J ij s ds + t 0 J ij s dV s + [X, J ij ] t , (3.63)
where [X, J ij ] t is the quadratic variation computed as follows Plugging into (3.63), the product X t J ij t becomes: Applying Ito product rule to X t J ij t and Λ k,θ t , we obtain

[X, J ij ] t = 0≤s≤t ∆J ij s ∆X s = 0≤s≤t ( X s , e i e ′ j ∆X s )(∆X s ) = t 0 X s , e i e ′ j dX s -X s , e i e ′ j X s -(dX s -X s -) = N l,
X t J ij t = t 0 X s , e i a ji dse i + t 0 X s , e i e ′ j dV s e i + t 0 AX s J ij s ds + t 0 J ij s dV s + N l,
X t J ij t Λ k,θ t = t 0 X s Λ k,θ s , e i a ji dse i + t 0 X s Λ k,θ s , e i e ′ j dV s e i + t 0 A ∆ X s J ji s Λ k,θ s ds + t 0 J ij s Λ k,θ s dV s + N l,m=1 t 0 X s Λ k,θ
s , e i e ′ j e l e i -X s Λ k,θ s , e i e ′ j e m e i , e m a lm ds(e le m )

+ N l,m=1 t 0 X s Λ k,θ
s , e i e ′ j e l e i -X s Λ k,θ s , e i e ′ j e m e i , e m dV s , e l (e le m )

+

N i=1 t 0 X s Λ k,θ s , e i J ij s c i β i dk s β i . (3.64)
The recursive equation of γ(X t J ij t ) is then obtained taking conditional expectation E given Y k t of both sides of (3.64):

γ(X t J ij t ) = t 0 q s , e i a ji dse i + t 0 A ∆ γ(X s J ij s )ds + N l,m=1 t 0
q s , e i e ′ j e l e iq s , e i e ′ j e m e i , e m a lm ds(e le m )

+ t 0 Cγ(X s J ij s ) dk s X s β = t 0 q s , e i a ji dse i + t 0 A ∆ γ(X s J ij s )ds + t 0 Cγ(X s J ij s ) dk s X s β + N l,m=1 t 0
q s , e i e ′ j e l e iq s , e i e ′ j e m e i , e m a lm dse l -N l,m=1 t 0 q s , e i e ′ j e l e iq s , e i e ′ j e m e i , e m a lm dse m

= t 0 q s , e i a ji dse i + t 0 A ∆ γ(X s J ij s )ds + t 0 Cγ(X s J ij s ) dk s X s β t 0
q s , e i a ji dse j -t 0 q s , e i a ji dse i .

Discrete-Time Approximate Model DTAM-k

Since the data that we deal with have been sampled at regular time intervals, the observation range [0, T ] is split into a regular partition

0 = t 0 < t 1 < t 2 < • • • < t h-1 < t h < . . . t M = T, (3.65)
where ∆ = t ht h-1 is such that M ∆ = T and M is the size of the observed sample k 1 , k 2 , . . . , k M used to estimate the model parameter vector.

Components of DTAM-k

Define the following discrete-time observations

u h = 1 ∆ (k h -k h-1 ), h = 1, 2, . . . .

From (3.28) u h verifies:

u h = X h , c + 1 √ ∆ X h , β w h , h = 1, 2 . . . , M, (3.66)
where c = (c 1 , c 2 , . . . , c N ) ′ , β = (β 1 , β 2 , . . . , β N ) ′ and w h i.i.d.

∼ N (0, 1), h = 1, 2 . . . . Equation (3.66) determines an observation sequence (u h ) h=1,2,... which is driven by state process (X h ) h=0,1,2,... , a discrete time Markov chain with state space S = {e 1 , e 2 , . . . , e N } and transition probability matrix:

P = I + ∆A ∆ = (ǫ ji ) ∈ R N ×N .
Note that, for a small enough time step ∆, P is a stochastic matrix. The parameter vector of this model is

θ = (ǫ ji , 1 ≤ j = j ≤ N ; c 1 , c 2 , . . . , c N ; β 1 , β 2 , . . . , β N ).
For i = 1, 2, . . . , N , let us define the level sums of DTAM-k by

T i h (u) h l=1 X l , e i u l , T i h (u 2 ) h l=1 X l , e i u 2 l .
Let Y u h = σ(u 1 , u 2 , . . . , u h ) denote the filtration of incomplete data. We will use the following filter processes to obtain ML estimates of the model parameters:

π(H h ) E θ (H h |Y u h ), γ(H h ) E(H h Λ u,θ h |Y u h ), γ(X h H h ) E(X h H h Λ u,θ h |Y u h ),
where

H h ≡ J ij h , O i h , T i h (u) or T i h (u 2 ).

Discrete-time approximate filtering equations

Our purpose here is to approximate the continuous-time filtering equations in order that the resulting equations are locally Lipschitz continuous given the observation path {k 1 , k 2 , . . . , k M }.

Approximation of state filter equation

For studying the approximate scheme of the state filtering equation, a particular case of the Duncan-Mortenson-Zakai equation, we follow Clark's robust filter approach [Clark 1977]. We first introduce the following processes:

φ (i) t = exp t 0 c i β i dk s β i - 1 2 t 0 c i β i 2 ds , i = 1, 2, . . . , N, (3.67) Φ t = diag φ (1) t , φ (2)
t , . . . , φ (N ) t , (3.68)

q t = Φ -1 t q t . (3.69)
The process q t is called the robust filter of the state process q t . We will use this process to provide an approximate equation for the state equation (3.46).

Applying Ito product rule to the two processes (φ

(i) t ) -1 = 1 - t 0 (φ (i) t ) -1 c i β i
dW s , by (3.31) and (3.36),

q (i) t = π (i) 0 + t 0 (λ i1 q (1) s + • • • + λ iN q (N ) s )ds + t 0 q (i) s c i β i dk s β i , from (3.46), for i = 1, 2, . . . , N , we get (φ (i) t ) -1 q (i) t = π (i) 0 + t 0 (φ (i) s ) -1 (λ i1 q (1) s + • • • + λ iN q (N ) s )ds (3.70)
and

q t = π 0 + t 0 (Φ s ) -1 A ∆ q s ds. (3.71)
From (3.71), we can determine (q t ) by the following ordinary differential equation:

d dt q = (Φ t ) -1 A ∆ Φ t q t , (3.72)
with initial condition q 0 = π 0 . Integrating (3.72) between sampling times t h-1 and t h of the partition (3.65) yields

q t h = q t h-1 + t h t h-1 (Φ s ) -1 A ∆ Φ s q s ds, (3.73)
and the Euler-Maruyama method [Higham 2001] shows that

q t h ≃ q t h-1 + (Φ t h-1 ) -1 A ∆ Φ t h-1 q t h-1 ∆. (3.74)
Multiplying both sides of (3.74) by Φ t h , we get

q t h ≃ Φ t h (Φ t h-1 ) -1 [I + ∆A ∆ ]q t h-1 = Φ t h (Φ t h-1 ) -1 P q t h-1 , (3.75)
where P = I + ∆A ∆ , with I is the N × N unit matrix.. We now consider the product

Φ t h (Φ t h-1 ) -1 .
The definition of φ

(i) t in (3.67) implies that φ (i) t h (φ (i) t h-1 ) -1 = exp t h t h-1 c i β 2 i dk s - 1 2 t h h-1 c i β i 2 ds, , (3.76)
for i = 1, 2, . . . , N . As the partition (3.65) is such that

∆ = t h -t h-1 , k h -k h-1 = ∆u h , (3.76) becomes φ (i) t h (φ (i) t h-1 ) -1 = exp ∆ 2β 2 i (2c i u h -c 2 i ) , i = 1, 2, . . . , N. (3.77)
We then have

Φ t h (Φ t h-1 ) -1 = diag B1 h , B2 h , . . . BN h = Bh , (3.78) where B(i) h = exp ∆ 2β 2 i (2c i u h -c 2 i )
, for i = 1, 2, . . . , N . Our numerical approximation of q t h will be denoted by qh . From (3.75) and (3.78), an approximation of the Duncan-Mortenson-Zakai equation is now given by qh = Bh P qh-1 , q0 = π 0 .

(3.79)

Clark [Clark 1977] has showed that the robust filter (q t ) determined by q t = Φ -1 t q t is a locally Lipschitz continuous function given the observed path {k s , 0 ≤ s ≤ t}. Therefore, applying (3.79), it is seen that the process (q m ) determined by equation (3.79) can be used to approximate the state filter (q t ) which inherits this continuity property.

Approximate filter equation of the number of jumps, of the occupation time and of the level sums

Using similar arguments as those used to obtain a robust filter of the state process q s , robust versions γ(X t H t ) of the filters γ(X t H t ) for the processes

H t ≡ J ij t , O i t , T i t can be obtained from d dt γ(X t J ij t ) =Φ -1 t A ∆ γ(X t J ij t )Φ t γ(X t J ij t ) + q, e j Φ -1 t A ∆ γ(X t J ij t )Φ t e j , e i e i , (3.80)
d dt γ(X t O i t ) =Φ -1 t A ∆ γ(X t O i t )Φ t γ(X t O i t ) + q, e i e i , (3.81) dγ(X t T i t ) =Φ -1 t A ∆ γ(X t T i t )Φ t γ(X t T i t )dt + q, e i e i dk t , (3.82) with initial conditions γ(X 0 J ij 0 ) = γ(X 0 O i 0 ) = γ(X 0 T i 0 (u 0 )) = 0, 1 ≤ i = j ≤ N . For H t ≡ J ij t , O i t or T i t , let γ(X t h H t h
) be an approximation of continuous-time filter processes γ(X t H t ) between sampling times t h-1 and t h in partition (??), using robust filtering equations (3.80), (3.81) and (3.82).

The filter processes γ(X t h J ji t h ), γ(X t h O i t h ) and γ(X t h T i t h ) will be approximated in the DTAM-k model by ∆-time filters γ(X h J ji h ), γ(X h O i h ) and γ(X h T i h (u h )), respectively, where: [James 1996, Theorem 3.3 and Corollary 3.4]:

γ(X h J ij h ) = Bh P γ(X h-1 J ij h-1 ) + qh-1 , e j Bh P e j , e i e i , (3.83) γ(X h O i h ) = Bh P γ(X h-1 O i h-1 ) + qh-1 , e i Bh P e i , (3.84) γ(X h T i h (u h )) = Bh P γ(X h-1 T i h-1 (u h )) + u h qh-1 , e i Bh P e i , (3.85) with initial conditions γ(X 0 J ij 0 ) = γ(X 0 O i 0 ) = γ(X 0 T i 0 (u 0 )) = 0, 1 ≤ i = j ≤ N . The error between γ(X h H h ), 1 (for H h ≡ J ij h , O i h or T i h , h = 1, 2, . . . , M ) and γ(X t h H t h ), 1 (for H t h ≡ J ij t h , O i t h or T i t h , 0 ≤ t h ≤ T ) is evaluated by
| γ(X t h H t h ), 1 [k] -γ(X h H h ), 1 [k]| ≤ C k |∆ + ω ∆ (k)|, (3.86)
where the constant C k depends on k sup

0≤t≤T |k t | and ω ∆ (k) = max {|k t -k s | : 0 ≤ s, t ≤ T, |t -s| ≤ ∆} .
(3.87)

Updating parameter

EM algorithm is also used to obtain numerical ML estimates for DTAM-k. At each iteration of algorithm, ML estimates of transition probabilities are updated by

ǫ ′ ji = γ(X M J ij M ), 1 γ(X M O i M ), 1 , 1 ≤ j = i ≤ N. (3.88)
For updating remaining parameters, we consider the following conditional expectation:

Q(θ ′ , θ) = E θ        log     h l=1 1 √ ∆ X l , σ φ u l -X l ,c ′ 1 √ ∆ X l ,σ ′ 1 √ ∆ X l , σ ′ φ u l -X l ,c 1 √ ∆ X l ,σ     Y u h        , (3.89) where c ′ = (c ′ 1 , c ′ 2 , . . . , c ′ N ) ′ , β ′ = (β ′ 1 , β ′ 2 , . . . , β ′ N ) ′ and φ(•) is the N (0, 1) density function.
At the p-th iteration of EM algorithm, vector θ = θ (p) and transform Q(θ ′ , θ (p) ) as follows

Q(θ ′ , θ (p) ) = N i=1 π(O i h ) log 1 β ′ i - ∆ β ′2 i π(T i h (u 2 h )) -2b ′ i π(T i h (u h )) + b ′2 i π(O i h ) + R(θ ′(p) , Y u h ),
where the function R(θ ′(p) , Y u h ) does not depend on θ ′ . Maximizing the function Q(θ ′ , θ (p) ) above yields

c ′ i = π(T i h (u)) π(O i h ) , i = 1, 2, . . . , N, (3.90)
β ′2 i = ∆ π(O i M ) π(T i M (u 2 )) -2c ′ i π(T i M (u)) + c ′2 i π(O i M ) . (3.91) For H h ≡ J ij h , O i h or T i h (u), by plugging π(H h ) = γ(H h ) γ(1) = γ(X h H h ,1 ) qh ,1
, into (3.90), (3.91), we get:

c ′ i = γ(X M T i M (u M )), 1 γ(X M O i M ), 1 , i = 1, 2, . . . , N (3.92) β ′2 i =∆ γ(X M O i M ), 1 -1 γ(X M T i M (u 2 )), 1 -2c ′ i γ(X M T i M (u)), 1 + c ′2 i γ(X M O i M ), 1 , (3.93)
where the discrete-time filter process γ(X h T i h (u 2 )) h=1,2,...,M is determined by: Iqbal 1983 proposed utilizing the magnitude of the daily clearness index K h to define sky conditions, see Table 3.2. Note that there are several criterias for solarday classification which depend on geographical location and climatic conditions. For expample, 4 classes were identified in Guadeloupe [START_REF] Soubdhan | [END_REF]]: clear sky days, intermittent clear sky days, intermittent cloudy sky days and cloudy sky days while in La Réunion 5 classes were found [Delsaut 2013]: cloudy days, intermittent bad days, disturbed days, intermittent good days and clear sky days. For simplicity, we here use the classification according to Iqbal 1983.

γ(X h T i h (u 2 h )) = Ψ h P γ(X h-1 T i h-1 (u 2 h )) + (u h ) 2 p h-1 , e i Ψ h P e i , (3.94) with initial condition γ(X 0 T i 0 (u 2 )) = 0, i = 1,

Estimations

CTM-k will be considered by estimating DTAM-k parameters with ∆ = 1/3600, using1-s mean CISs: k-DATA-I.1, k-DATA-II.2 and k-DATA-III.3.

Observing the histograms of these data The histograms have incited to take N = 4 states, the parameter vector is

θ = {(a ji ), 1 ≤ j = i ≤ 4; b, α}, with b = (b 1 , b 2 , b 3 , b 4 ) ′ and α = (α 1 , α 2 , α 3 , α 4 ) ′ .
The transition probability matrix is

A =    
a 11 a 12 a 13 a 14 a 21 a 22 a 23 a 24 a 31 a 32 a 33 a 34 a 41 a 42 a 43 a 44     .

Note that a ii = 1j =i a ji , for i = 1, 2, 3, 4. After 50 iterations, we obtain the following estimates:

1. From K-DATA-0809:

A (50) =     0.
= (0.0586, 0.0473, 0.0167, 0.0154) ′ .

2. From K-DATA-0310: Different limit estimations (θ 1, * and θ 2, * ) are obtained from the same observation data when using different initial parameters (θ (1,0) and θ (2,0) , respectively). However, simulated data obtained from the models with these limit vectors do not have many differences. Simulated values generated by these models have the same distribution as observation data. The mean of simulated paths converge to the mean of observation data. These remarks confirm the properties of EM algorithm used in our estimatiion setting.

A (50) =     0
We checked this problem in several cases. For instance, to estimate the parameters of DTM-K with N = 4 states by using the daily CISs observed in October 2011 at La Réunion ( Figure 3.8(a)), we started with two initial parameter vectors θ (1,0) , θ (2,0) as follows:

The vector θ (1,0) includes :

The vector θ (2,0) includes : Results after 100 iterations:

A (1,0) =     0.3 0.3 0.2 0.3 0.4 0.4 0.4 0.4 0.2 0.1 0.3 0.1 0.1 0.2 0.1 0.2     , A (
1. The limit vector θ (1, * ) obtained from θ (1,0) :

A (1,100) =    
0.2192 0.3714 0.0000 0.1646 0.4008 0.2282 0.5722 0.5197 0.3800 0.4005 0.1956 0.0000 0.0000 0.0000 0.2322 0.3158

    ,
b (1,100) =(0.2793, 0.4312, 0.6392, 0.7004) ′ , α (1,100) =(0.0522, 0.0985, 0.0343, 0.0050) ′ .

2. The limit vector θ (2, * ) obtained from θ (2,0) :

A (2,100) =    
0.2223 0.0000 0.0000 0.5443 0.0000 0.1095 0.7774 0.1528 0.7777 0.0000 0.2226 0.0000 0.0000 0.8905 0.0000 0.3029 2,100) =(0.2860, 0.6459, 0.3818, 0.5567) ′ , α (2,100) =(0.0428, 0.0381, 0.1173, 0.1139) ′ .

    , b (
The graphs in Figure 3.12, Figure 3.13 and Figure 3.14 show the evolution of estimates from θ (1,0) and θ (2,0) . Figure 3.8(b) shows a simulated path from the model with limit vector parameter θ (1, * ) . Convergence of averages and estimated cumulative distribution function of simulated values generated by the models with these limit vector parameter θ (1, * ) , θ (2, * ) are shown in Figure 3.15. For CTM-k, we start with

A ∆(0) =     -13 3 2 3 4 -9 4 3 2 3 -10 1 7 3 4 -7     , c (0) =(0.25, 0.35, 0.55, 0.85) ′ , β (0) =(1, 2, 3, 4) ′ .
then estimates of transition matrix in its DTAM-k will be started with:

P = I + ∆A ∆ =    
0.9964 0.0008 0.0006 0.0008 0.0011 0.9975 0.0011 0.0008 0.0006 0.0008 0.9972 0.0003 0.0019 0.0008 0.0011 0.9981

    .
After 50 iterations of EM algorithm, we obtain 1. From k-DATA-I.1: 3. From k-DATA-III.3:

P (50) =     0.
P (50) =    
0.9941 0.0041 0.0073 0.0000 0.0046 0.9940 0.0000 0.0133 0.0013 0.0000 0.9927 0.0000 0.0000 0.0019 0.0000 0.9867 Chapter 4

    , c (
A Stochastic model for the total solar radiation process Résumé Dans ce chapitre, nous proposons un HMM à temps continu, noté CTM-y, pour modéliser le rayonnement solaire total (y t ) t≥0 considéré comme un processus aléatoire en milieu aléatoire.

Comme pour CTM-k, le processus détat (X t ) est une chaîne de Markov à temps continu dont les états représentent les divers régimes météorologiques mais maintenant, le processus d'observation (y t ) des mesures du rayonnement solaire total est solution de l'équation différentielle stochastique en milieu aléatoire suivante :

dy t = [g(X t )I t -y t ]dt + σ(X t )y t dW t ,
où I t désigne le rayonnement extraterrestre à l'instant t et (W t ) est un mouvement Brownien standard.

Nous détaillons la méthode de probabilité de référence avec la technique de changement de mesure et nous établissons les équations de filtrage pour la mise à jour de l'estimation du vecteur de paramètres dans l'algorithme EM. Nous faisons de même pour DTAM-y, un modèle approché à temps discret.

Les applications numériques, faites à partir de mesures échantillonnées à 1Hz, permettent la. Les travaux de ce chapitre ont été présentés à [START_REF][END_REF]] et dans [Tran b].

Abstract

In the present chapter, we propose a continuous-time HMM, denoted CTM-y, for modelling the total solar radiation (y t ) t≥0 considered as a stochastic process in random medium.

As for CTM-k, the state process (X t ) is a continuous-time Markov chain whose states represent various meteorological regimes but now, the total solar radiation observation process (y t ) is solution of the following stochastic differential equation in random medium:

dy t = [g(X t )I t -y t ]dt + σ(X t )y t dW t ,
where I t denotes the extraterrestrial radiation at time t and (W t ) a standard Brownian motion.

We detail the reference probability method with the change of measure technic and we establish the filter equations for updating the parameter vector estimation in EM algorithm. We do the same for DTAM-y, an approximating discrete-time model.

Experiments are done using measurements sampled at 1Hz. The works of this chapter have been presented in [START_REF][END_REF]] and in [Tran b].

CTM-y

State process

Assuming that there are N ≥ 1 regimes (states) in meteorological dynamics, we will model the state dynamics by a continuous-time homogeneous Markov chain (X t ) 0≤t≤T with state space S = {e 1 , e 2 , ..., e N }, transition rate matrix A ∆ = (λ ji ) ∈ R N ×N and verifying the following equation

X t = X 0 + t 0 A ∆ X s ds + V t , (4.1)
where (V t ) is a (F t )-martingale.

Pseudo-clearness index

Recall that clearness index k t is a dimensionless indicator of total solar radiation G t : k t = Gt It , where I t is the extraterrestrial radiation. We now assume that each meteorological regime i = 1, . . . , N induces a level of pseudo-clearness index, say g i . Let g = (g 1 , g 2 , . . . , g N ) ′ and assume that the g i 's are distinct.

The dynamics of such levels is descibed through the pseudo-clearness index function:

g(X t ) = N i=1 X t , e i g i = X t , g . (4.2)

Observation process

During regime i, we assume that the process (y t ) is solution of a SDE with an almost affine drift depending on both pseudo-clearness index and extraterrestrial radiation, the change dy t generally being likely to be negative but not lower than a threshold: The proposed model CTM-y that we propose for the total solar radiation process is the pair of processes (X t , y t ) (Figure 4.1) defined by equations (4.1) and (4.5), with parameter vector:

X t ,
θ := {λ ji , 1 ≤ i = j ≤ N ; g 1 , g 2 , . . . , g N ; σ 1 , σ 2 , . . . , σ N } . (4.6)
The number of states N will be suggested by the user and our aim is to estimate the parameter vector θ given the observed data.

Filtrations

Let G y t = σ{X s , y s ; 0 ≤ s ≤ t} and Y y t = σ{y s ; 0 ≤ s ≤ t} denote the filtrations of the complete data and the incomplete data, respectively. These filtrations contain all the available information of observation history up to time t and Y y t will be used to estimate the parameters of the model.

Change of measure

Let us now detail the reference probability method with the change-of-measure technique.

By Ito formula

d log y t σ(X t ) = g(X t )I t -y t σ(X t )y t dt + dW t . (4.7)
Working on the reference probability space (Ω, (F t ), P ) with a fix probability measure P , we will frst establish the equations of filter processes. Then we will translate these results in the "real world" (Ω, (F t ), P θ ) with a change of probability measure using Girsanov theorem.

Using Ito formula in the "real world" (Ω, (F t ), P θ ) the observation process equation (4.5) has the equivalent form:

d log y t σ(X t ) = g(X t )I t -y t σ(X t )y t dt + dW t . (4.8)
Consider a change of measure from P θ to P such that

P P θ G y t = Λ y,θ t = exp - t 0 g(X s )I s -y s σ(X s )y s dW s - 1 2 t 0 g(X s )I s -y s σ(X s )y s 2 ds . (4.9)
By Girsanov theorem and (4.8), log yt σ(Xt) is a standard Brownian motion under P , a more convenient measure to work with.

We can come back to the "real world" with the initial probability measure P θ by inverse change:

P θ P G y t = Λ y,θ t = exp t 0 g(X s )I s -y s σ(X s )y s d log y s σ(X s ) - 1 2 t 0 g(X s )I s -y s σ(X s )y s 2 ds .
(4.10) Girsanov theorem implies that under P θ , ( Wt ) is a standard Brownian motion if we define

d Wt = d log y t σ(X t ) - g(X t )I t -y t σ(X t )y t dt. (4.11)
It follows that under P θ ,

dy t = [g(X t )I t -y t ]dt + σ(X t )y t d Wt , (4.12)
where ( Wt ) is a standard Brownian motion. That is, P θ is the distribution law of the observation process (y t ) in the "real world". Note that:

1. By Ito formula and (4.10), the process Λ θ t is solution of the SDE:

Λ y,θ t = 1 - t 0 Λ s g(X s )I s -y s σ(X s )y s dW s . (4.13) 2. Similarly, Λ y,θ t =1 + t 0 Λ θ s g(X s )I s -y s σ(X s )y s d log y s σ(X s ) (4.14) =1 + t 0 Λ θ s g(X s )I s -y s σ(X s )y s 2 ds + t 0 Λ θ s g(X s )I s -y s σ(X s )y s dW s . (4.15) 3. We have Λ y,θ t Λ y,θ t = 1.

Parameter estimations in continuous time

As for CTM-k, the estimates of transition probabilities λ ji (1 ≤ i = j ≤ N ) are given by [Elliott 1995, chap. 3, pg. 68-70] and [START_REF] James | [END_REF]]:

λ ′ ji = π(J ji t ) π(O i t ) , 1 ≤ j = i ≤ N (4.16)
where π(J ji t ) and π(O i t ) are the normalized filter (see Section 2.8.2). We now proceed to the computations in EM algorithm for getting ML estimates of the parameter vector g = (g 1 , g 2 , . . . , g N ) ′ .

Expectation Step

Update from parameter θ to new parameter θ ′ will be obtained by maximizing the pseudo log-likelihood function defined as

Q(θ ′ , θ) = E θ log Λ y,θ ′ θ t |Y y t , (4.17)
where

Λ y,θ ′ θ t = dP θ ′ dP θ G y t .
First, we have:

Q(θ ′ , θ) =E θ t 0 X s , g ′ I s -y s ( X s , σ ′ y s ) 2 - X s , g I s -y s ( X s , σ y s ) 2 dy s - 1 2 t 0 X s , g ′ I s -y s X s , σ ′ y s 2 - X s , g I s -y s X s , σ y s 2 ds Y y t , (4.18)
where

g ′ = (g ′ 1 , g ′ 2 , . . . , g ′ N ) ′ , σ ′ = (σ ′ 1 , σ ′ 2 , . . . , σ ′ N ) ′ . At the p-th iteration of the algorithm Expectation Step, we set θ = θ ′(p) in (4.18) to obtain Q(θ ′ , θ ′(p) ) = E θ ′(p) N i=1 1 σ ′2 i g ′ i t 0 I s y 2 s X s , e i dy s - t 0 1 y s X s , e i dy s - 1 2 g ′2 i t 0 I 2 s y 2 s X s , e i ds + g ′ i t 0 I s y s X s , e i ds - 1 2 t 0 X s , e i ds Y y t + R(θ ′(p) , Y y t ), (4.19)
where R(θ ′(p) , Y y t ) does not depend on θ ′ .

Maximization Step

The maximum of Q(θ ′ , θ ′(p) ) is achieved by solving

∂ ∂θ ′ Q(θ ′ , θ ′(p) ) = 0.
We will then find a ML estimate g ′ i for updating g i , i = 1, 2, . . . , N as follows. We first see that

∂ ∂g ′ i Q(θ ′ , θ ′(p) ) = 1 σ ′2 i E θ ′(p) t 0 I s y 2 s X s , e i dy s Y y t -g ′ i E θ ′(p) t 0 I 2 s y 2 s X s , e i ds Y y t + E θ ′(p) t 0 I s y s X s , e i ds Y y t . (4.20)
Let

M i t = t 0 I s y s X s , e i ds, (4.21)
N i t = t 0 I 2 s y 2 s X s , e i ds, (4.22)
T i t = t 0 I s y 2 s X s , e i dy s . (4.23) Then ∂ ∂g ′ i Q(θ ′ , θ ′(p) ) = 0 ⇐⇒ g ′ i = E θ ′(p) (M i t |Y y t ) + E θ ′(p) (T i t |Y y t ) E θ ′(p) (N i t |Y y t ) . ( 4
g ′ i = π(M i T ) + π(T i T ) π(N i T ) , i = 1, 2, . . . , N, (4.25)
Now, we have to compute π(H

t ) = E θ ′(p) (H t |Y y t ), for H t ≡ J ij t , O i t , M i t , N i t and T i t .

Equation of continuous time filters

As mentioned briefly in Chapter 2 and as applied in Chapter 3, we can obtain the normalized filter π(H

t ) = E θ ′(p) (H t |Y y t ) from the unnormalized filter γ(H t ) = E(H t Λ y,θ t |Y y t ) by noticing that π(H t ) = γ(H t ) γ(1) (4.26) Now our aim is to compute γ(H t ) for H t ≡ 1, J ij t , O i t , M i t , N i t or T i t .
Recall that it is not generally possible to obtain directly the equations for γ(H t ), see [START_REF] James | [END_REF]], but, it is possible to obtain the equations for γ(H t X t ) and then get γ(H t ) as γ(H t ) = γ(H t X t ), 1 , where 1 denotes the column N -vector of ones.

In particular, for H t ≡ 1, p t = γ(X t ), and γ(1) = γ(X t ), 1 = p t , 1 so that (4.26) yields

π(H t ) = γ(H t X t ), 1 p t , 1 (4.27)
A recursive equation for γ(X t H t ) for H t ≡ 1, J ij t , O i t , M i t , N i t and T i t is stated in the following theorem.

Theorem 4.1. Let H t a scalar process of the form

H t = H 0 + t 0 µ s ds + t 0 ρ ′ s dV s + t 0 δ s dW s , (4.28)
where the real µ s , δ s and N -dimensional vector ρ s are (F t )-predictable, squareintegrable processes.

The unnormalized estimate γ(H t X t ) is given by

γ(H t X t ) =γ(H 0 X 0 ) + t 0 A ∆ γ(H s X s )ds + t 0 γ(µ s X s )ds + N i,j=1 t 0 γ(ρ j s X s -ρ i s X s ), e i a ji ds(e j -e i ) + t 0 C s γ(X s H s ) dy s σ(X s )y s + t 0 γ(δ s X s ) dy s σ(X s )y s . (4.29)
Proof. Applying Ito product rule for semimartingales, compute the product H t X t (where the process X t is given by (2.36)) as follows:

H t X t =H 0 X 0 + t 0 H s A ∆ X s ds + t 0 H s dV s + t 0 µ s X s ds + t 0 ρ ′ s X s dV s + t 0 δ s X s dB s + [X, H] t (4.30)
Since ∆H t = ρ ′ t ∆V t + ρ ′ t ∆W t , ∆X t = ∆V t and ∆X t ∆W t = 0 a.s., we have:

[X, H] t = 0<s≤t (ρ ′ s ∆X s )∆X s = t 0 ρ ′ s dX s -ρ ′ s X s (dX s -X s ) = N i,j=1 t 0 (ρ j s -ρ i s ) X s , e i dX s , e j (e j -e i ).
Then, using dX s = AX s ds + dV s , we get:

[X, H] t = N i,j=1 t 0 (ρ j s -ρ i s ) X s , e i dV s , e j (e j -e i ) + N i,j=1 t 0 ρ j s X s -ρ i s X s a ji ds(e j -e i )
Plugging in (4.30) yields

H t X t =H 0 X 0 + t 0 H s A ∆ X s ds + t 0 H s dV s + t 0 µ s X s ds + t 0 ρ ′ s X s dV s + t 0 δ s X s dW s (4.31) + N i,j=1 t 0 (ρ j s -ρ i s ) X s , e i dV s , e j (e j -e i ) + N i,j=1 t 0 ρ j s X s -ρ i s X s , e i a ji ds(e j -e i ) (4.32)
Now, Ito product applied to the processes Λ y,θ t and H t X t yields:

Λ y,θ t H t X t =H 0 X 0 + t 0 Λ y,θ s H s A ∆ X s ds + t 0 Λ y,θ s H s dV s + t 0 µ s Λ y,θ s X s ds + t 0 ρ ′ s Λ y,θ s X s dV s + t 0 δ s Λ y,θ s X s dW s (4.33)
+ N i,j=1 t 0 (ρ j s -ρ i s ) Λ y,θ s X s , e i dV s , e j (e j -e i ) + N i,j=1 t 0 ρ j s Λ y,θ s X s -ρ i s Λ y,θ s X s , e i a ji ds(e j -e i ) + t 0 Λ y,θ s X s H s g(X s )I s -y s σ(X s )y s dy s σ(X s )y s + [HX, Λ y,θ ] t (4.34)
where the quadratic variation is given by

[XH, Λ y,θ ] t = t 0 δ s Λ y,θ s X s g(X s )I s -y s σ(X s )y s ds (4.35)
Plugging (4.35) into (4.34), we get

t 0 Λ y,θ s X s H s g(X s )I s -y s σ(X s )y s dy s σ(X s )y s = N i=1 t 0 Λ y,θ s X s H s , e i g i I s -y s σ i y s dy s σ i y s e i , (4.36)
t 0 δ s Λ y,θ s X s dB s + t 0 δ s Λ y,θ s X s g(X s )I s -y s σ(X s )y s ds = t 0 δ s Λ y,θ s X s dy s σ(X s )y s = N i=1 t 0 Λ y,θ s δ s X s , e i dy s σ i y s e i . (4.37) Hence Λ y,θ t H t X t =H 0 X 0 + t 0 Λ y,θ s H s A ∆ X s ds + t 0 Λ y,θ s H s dV s + t 0 µ s Λ y,θ s X s ds + t 0 ρ ′ s Λ y,θ s X s dV s + N i,j=1 t 0 (ρ j s -ρ i s ) Λ y,θ s X s , e i dV s , e j (e j -e i ) + N i,j=1 t 0 ρ j s Λ y,θ s X s -ρ i s Λ y,θ s X s , e i a ji ds(e j -e i ) + N i=1 t 0 Λ y,θ s X s H s , e i g i I s -y s σ i y s dy s σ i y s e i + N i=1 t 0 Λ y,θ s δ s X s , e i dy s σ i y s e i . (4.38)
Taking conditional expectation E given Y y t of both sides of (4.39) we obtain

γ(H t X t ) =γ(H 0 X 0 ) + t 0 A ∆ γ(H s X s )ds + t 0 γ(µ s X s )ds + N i,j=1 t 0 γ(ρ j s X s -ρ i s X s ), e i a ji ds(e j -e i ) + N i=1 t 0 γ(X s H s ), e i g i I s -y s σ i y s dy s σ i y s e i + N i=1 t 0 γ(δ s X s ), e i dy s σ i y s e i . (4.39)
As

N i=1 t 0 γ(δ s X s ), e i dy s σ i y s e i = t 0 γ(δ s X s ) dy s σ(X s )y s , (4
t = γ(X t ) = E(X t Λ y,θ t |Y y t ) is: p t = π 0 + t 0 A ∆ p s ds + t 0 C s p s dy s σ(X s )y s (4.40)
As in (3.62), we decompose the number of jumps from e i to e j in the time interval [0, t] as follows:

J ij t = t 0 X s , e i a ji ds + t 0 X s , e i e ′ j dV s
Then, applying Theorem 4.1 with H t = J ij t , H 0 = 0, µ s = X s , e i a ji , ρ = X s , e i e j , δ s = 0, we get:

Corollary 4.2. γ(J ij t X t ) = t 0 A ∆ γ(J ij t X s )ds + t 0 a ji p s , e i e j ds + t 0 C s γ(J ij s X s ) dy s σ(X s )y s (4.41) Note here that N i,j=1 t 0 γ(ρ j s X s -ρ i s X s
), e i a ji ds(e je i ) = t 0 a ji p s , e i e j ds -t 0 a ji p s , e i e i ds. because ρ j s = X s , e i , ρ i s = 0. Recalling that X s X s , e i = X s , e i e i and applying Theorem 4.1 to the processes

O i t = t 0 X s , e i ds, M i t = t 0 Is ys X s , e i ds and N i t = t 0 I 2 s y 2 s
X s , e i ds, i = 1, 2, . . . , N , we obtain the recursive equations for γ(O i t X t ), γ(M i t X t ) and γ(N i t X t ), respectively, as stated in Corollary 4.3, Corollary 4.4 and Corollary 4.5 below.

Take

H t = O i t , H 0 = 0, µ s = X s , e i , δ s = 0, ρ = 0, to get Corollary 4.3. γ(O i t X t ) = t 0 A ∆ γ(O i s X s )ds + t 0 C s γ(O i s X s ) dy s σ(X s )y s + t 0 p s , e i e i ds (4.42) Take H t = M i t , H 0 = 0, µ s = Is ys X s , e i , δ s = 0, ρ = 0, to get Corollary 4.4. γ(M i t X t ) = t 0 A ∆ γ(M i s X s )ds + t 0 C s γ(M i s X s ) dy s σ(X s )y s + t 0 I s y s p s , e i e i ds (4.43) Take H t = N i t , H 0 = 0, µ s = I 2 s y 2 s X s , e i , δ s = 0, ρ = 0, to get Corollary 4.5. γ(N i t X t ) = t 0 A ∆ γ(N i s X s )ds + t 0 C s γ(N i s X s ) dy s σ(X s )y s + t 0 I 2 s y 2 s p s , e i e i ds (4.44) Consider now T i t = t 0 Is y 2 s
X s , e i dy s , i = 1, 2, . . . , N and observe that

T i t = t 0 I s y 2 s (I s g i -y s ) X s , e i ds + t 0 I s y s σ i X s , e i dB s . (4.45) Applying Theorem 4.1 with H t = T i t , H 0 = 0, µ s = Is y 2 s (I s g i -y s ) X s , e i , ρ = 0, δ s = Is ys σ i X s , e i , we obtain γ(T i t X t ) = t 0 A ∆ γ(T i s X s )ds + t 0 C s γ(T i s X s ) dy s σ(X s )y s + t 0 I s y 2 s (I s g i -y s ) p s ,
γ(T i t X t ) = t 0 A ∆ γ(T i s X s )ds + t 0 C s γ(T i s X s ) dy s σ(X s )y s + t 0 I s y 2 s (I s g i -y s ) p s , e i e i ds + t 0 I s y 2 s p s , e i e i dy s .
(4.47)

Discrete-time approximating model DTAM

In this section we consider partition (3.65

) of time interval [0, T ] (0 = t 0 < t 1 < t 2 < • • • < t m-1 < t m < • • • < t M ≤ T with equal width ∆ = t m -t m-1
) to obtain an approximating model of Continuous-Time Model (CTM) and the computable approximations of the continuous-time filter equations described in Section 4.2.

Components of the model

From partition (??), we first define the following discrete-time sampled observations:

z m = 1 ∆ (y m -y m-1 ), m = 1, 2 . . . (4.48)
The CTM-y discrete-time approximating model is the pair of processes (X m , z m ) m=1,2,... defined as follows:

1. the state process is a discrete time Markov chain (X h ) h=0,1,2,... with state space S = {e 1 , e 2 , . . . , e N } and transition probability matrix:

P = I + ∆A ∆ = (ǫ ji ) ∈ R N ×N ;
2. the state process (X m ) is related to the observation sequence (z m ) m=1,2,..., through the equation:

z m = [ X m , g I m -y m ] + 1 √ ∆ X m , σ y m w m , m = 1, 2 . . . , (4.49)
where w m i.i.d.

∼ N (0, 1), m = 1, 2 . . . ;

3. the parameter vector of the model is:

θ = {ǫ ji , 1 ≤ j = i ≤ N ; g 1 , g 2 , . . . , g N ; σ 1 , σ 2 , . . . , σ N }.

Robust approximation of filter equations

We propose now some computable equations approximating the continuous-time equations described in Section 4.2.

For the approximating schema of the state filter equation, we follow again Clark robust filter approach [Clark 1977] using the robust forward state filter density introduced by [START_REF] James | [END_REF][START_REF] Krishnamurthy | [END_REF].

We first introduce the following processes:

φ (i) t = exp t 0 g i I s -y s σ i y s dy s σ i y s - 1 2 t 0 g i I s -y s σ i y s 2 ds , i = 1, 2, . . . , N, (4.50)
Φ t = diag φ (1) t , φ (2)
t , . . . , φ

(N ) t . (4.51) Definition 4.1. For t ∈ [0, T ], define p t (Φ t ) -1 p t . (4.52)
The process (p t ) is called the robust filter of the state filter process (p t ).

Considering similarly as in (4.13) and (4.14), we have (φ

(i) t ) -1 = 1 - t 0 (φ (i) t ) -1 g i I s -y s σ i y s dW s , i = 1, 2, . . . , N. (4.53)
From the state filter equation,

p (i) t = π (i) 0 + t 0 (λ i1 p (1) s + • • • + λ iN p (N ) s )ds+ t 0 p (i) s g i I s -y s σ i y s dy s σ i y s , i = 1, 2, . . . , N. (4.54)
By Itô rule, we have (φ

(i) t ) -1 p (i) t = π (i) 0 + t 0 (φ (i) s ) -1 (λ i1 p (1) s + • • • + λ iN p (N ) s )ds, i = 1, 2, . . . , N. (4.55)
An equivalent form of (4.52) is therefore:

p t = (Φ t ) -1 p t = π 0 + t 0 (Φ s ) -1 A ∆ p s ds. (4.56)
By (4.56), the process (p t ) will be determined by the ordinary differential equation:

d dt p t = (Φ t ) -1 A ∆ Φ t p t , p 0 = π 0 . (4.57)
An approximating schema to get the equation of (p t ) from (p t ) will be considered below.

From (4.57), an approximation of p t between the sampling times t m-1 and t m of the partition (3.65) is given by

p tm = p t m-1 + tm t m-1 (Φ s ) -1 A ∆ Φ s p s ds.
(4.58)

Euler-Maruyama method (see e.g. [Higham 2001]) yields

p tm ≃ p t m-1 + (Φ t m-1 ) -1 A ∆ Φ t m-1 p t m-1 ∆. (4.59)
Multiplying both sides of (4.59) by Φ tm , we obtain

p tm ≃ Φ tm (Φ t m-1 ) -1 [I + ∆A ∆ ]p t m-1 , (4.60)
where I is the N × N unit matrix.

Plugging à = I + ∆A ∆ into (4.60) yields

p tm ≃ Φ tm (Φ t m-1 ) -1 Ãp t m-1 . ( 4
φ (i) tm (φ (i) t m-1 ) -1 = exp tm t h-1 g i I s -y s (σ i y s ) 2 dy s - 1 2 tm m-1 g i I s -y s σ i y s 2 ds . (4.62)
As the partition (3.65) is such that ∆ = t mt m-1 , y my m-1 = ∆z m , (4.62) becomes

φ (i) tm (φ (i) t m-1 ) -1 = exp ∆ (g i I m -y m ) z m (σ i y m ) 2 - 1 2 g i I h -y m σ i y m 2 = ψ (i) m i = 1, 2, . . . , N. (4.63)
We then have

Φ tm (Φ t m-1 ) -1 = diag ψ (1) m , ψ (2) m . . . ψ (N ) m = Ψ m , (4.64)
where

ψ (i) m = exp ∆ (g i I m -y m ) z m (σ i y m ) 2 - 1 2 g i I m -y m σ i y m
Similarly to Section 3.3.2.2, the continuous-time filtering processes γ(J ji t X t ), γ(O i t X t ), γ(M i t X t ), γ(N i t X t ) and γ(T i t X t ) in CTM-y will be approximated by the processes γ(J ji m X m ), γ(O i m X m ), γ(M i m X m ), γ(N i m X m ) and γ(T i m X m ), respectively, if they are determined as follows:

γ(J ji m X m ) = Ψ m Ãγ(J ji m-1 X m-1 ) + α ji q m-1 , e i ψ j m e j , (4.66) γ(O i m X m ) = Ψ m Ãγ(O i m-1 X m-1 ) + q m-1 , e i Ψ m Ãe i , (4.67) γ(M i m X m ) = Ψ m Ãγ(M i m-1 X m-1 ) + I m y m q m-1 , e i Ψ m Ãe i , (4.68) γ(N i m X m ) = Ψ m Ãγ(N i m-1 X m-1 ) + I m y m 2 q m-1 , e i Ψ m Ãe i , (4.69) γ(T i m X m ) = Ψ m Ãγ(T i m-1 X m-1 ) + z m I m y 2 m q m-1 , e i Ψ m Ãe i , (4.70) with initial conditions γ(J ji 0 X 0 ) = γ(O i 0 X 0 ) = γ(M i 0 X 0 ) = γ(N i 0 X 0 ) = γ(T i 0 X 0 ) = 0,
where 0 denotes the column-vector of N zeros and (q m ) is the state filter approximation obtained from (4.65). By using the discrete-time approximate filter processes above, updates, in our DTAM-y model, from ǫ ji to ǫ ′ ji (1 ≤ j = i ≤ N ) and from g i to g i (i = 1, 2, . . . , N ) are given by

ǫ ′ ji = γ(J ji M X M ), 1 γ(O i M X M ), 1
, (4.71)

g ′ i = γ(T i M X M ), 1 + γ(M i M X M ), 1 γ(N i M X M ), 1 . (4.72)

Estimation of the noise variance

In DTAM-k model, ML estimation of the noise variance can be obtained by using Radon-Nikodym derivative [START_REF] James | [END_REF]]:

dP θ ′ dP θ = M h=1 1 1 √ ∆ y h X h , σ ′ √ 2π exp - (z h -I h X h , g ′ + y h ) 2 2 (y h X h ,σ ′ ) 2 ∆ × 1 1 √ ∆ y h X h , σ √ 2π exp - (z h -I h X h , g + y h ) 2 2 (y h X h ,σ ) 2 ∆ -1 . (4.73)
The conditional expectation of the log-likelihood function is determined by

Q(θ ′ , θ) =E θ log dP θ ′ dP θ |Y z m =E θ m h=1 -log X h , σ ′ - 1 2 ∆ (y h X h , σ ′ ) 2 z 2 h + I 2 h X h , g ′ 2 + y 2 h -2z h I h X h , g ′ + 2z h y h -2y h I h X h , g ′ Y z m + R(θ, Y z m ) =E θ m h=1 - N i=1 X h , e i log σ ′ i - ∆ 2 N i=1 X h , e i 1 σ ′2 i z 2 h y 2 h + I 2 h y 2 h g ′2 i + 1 -2 z h I h y 2 h g ′ i + 2 z h y h -2 I h y h g ′ i Y z m + R(θ, Y z m ) = N i=1 -π(O i m ) log σ ′ i -
where the discrete-time filter processes γ

(O i M X M ), γ(M i M X M ), γ(N i M X M ), γ(T i M X M )
(L i M X M ), γ(Q i M X M ) are computed from γ(L i m X m ) =Ψ m ÃγL i m-1 X m-1 ) + z m y m 2 q m-1 , e i Ψ m Ãe i , (4.77) γ(Q i m X m ) =Ψ m Ãγ(Q i m-1 X m-1 ) + z m y m q m-1 , e i Ψ m Ãe i , m = 1, 2, . . . , M, (4.78) with initial condition γ(L i 0 X 0 ) = γ(Q i 0 X 0 ) = 0.

Experiments with real data

We proceed now to parameter estimations with all the total solar radiation sequences mentioned in Table 3.3. Various shapes of data histograms being observed, we tried different numbers of states: N = 2, 4, 6, . . . . The sequence of ML estimates {θ (p) , p = 1, 2, . . . } generated by EM algorithm fastly converges to a stability point when choosing reasonable initial parameter vectors and correct number of states.

Here are some illustrations when dealing with two 4 hours-long data sequences:

• y-DATA-I.4, a total solar radiation sequence of type I (clear day) recorded from 06:00 to 10:00 during the 262th-day, 2006, in Guadeloupe. Clearness index mean and variance are 0.7564, 0.0014, repectively

• y-DATA-III.4, a total solar radiation sequence of type III (cloudy day) recorded from 13:00 to 17:00 during the 180th-day, 2006, in Guadeloupe. Clearness index mean and variance are 0.1456, 0.0018, respectively.

Similarly, we proceed with two 1 hour-long sequences of type II (partially cloudy day): y-DATA-II.1, y-DATA-II.2 (see Table 3.3).

Observing these data histograms ( The parameters of the CTM-y will be approximated by the estimates of its DTAM-y from the initial parameter vector {P (0) , g (0) , σ (0) }, where

P (0) = I + ∆A (0) =    
0.9917 0.0036 0.0019 0.0022 0.0031 0.9892 0.0019 0.0028 0.0025 0.0036 0.9944 0.0025 0.0027 0.0036 0.0017 0.9925

    ,
and ∆ = 1/3600.

The estimated parameters after 100 iterations are (from the 30-th iteration, the estimated parameters are quite similar, see 1. For y-DATA-I.4:

P (100) =    
0.9661 0.0000 0.0000 0.0000 0.0339 0.9718 0.0000 0.0000 0.0000 0.0282 0.9812 0.0000 0.0000 0.0000 0.0188 1.0000 (20.5595, 15.2412, 12.6386, 1.5645) ′ , σ (100) = (0.5107, 0.4960, 0.4993, 0.1346) ′ .

    , g (100) =
2. For y-DATA-III.4:

P (100) =    
0.9905 0.0000 0.0000 0.0007 0.0000 0.9833 0.0013 0.0000 0.0000 0.0167 0.9934 0.0013 0.0095 0.0000 0.0053 0.9980 100) = (-0.7715, 1.7319, 0.6673, 0.0638) ′ , σ (100) = (0.0923, 0.0881, 0.0782, 0.0551) ′ . We now consider numerical estimations based on the data sequences y-DATA-II.1 (Figure 4.9a1)) and y-DATA-II.2 (Figure 4.9b1)). We started with the following initial parameters: 0) = (0.15, 0.45, 0.65, 0.95) ′ , σ (0) = (1.5, 2.5, 3.5, 4.5) ′ .

    , g (

The graphs in

A (0) =     -7 4 3 2 3 -13 3 4 3 2 -9 4 1 7 3 -10     , g (
and the matrix:

P (0) = I + ∆A (0) =    
0.9981 0.0011 0.0008 0.0006 0.0008 0.9964 0.0008 0.0011 0.0008 0.0006 0.9975 0.0011 0.0003 0.0019 0.0008 0.9972

    ,
where ∆ = 1/3600. The estimates obtained after 100 iterations:

1. From y-DATA-II.1:

P (100) =    
0.9418 0.1219 0.0000 0.1186 0.0212 0.8781 0.0060 0.0000 0.0070 0.0000 0.9843 0.0000 0.0300 0.0000 0.0096 0.8814 100) = (-30.7737, 92.3213, 3.5591, 94.4960) ′ , σ (100) = (3. 4343, 3.6536, 1.3474, 3.6932).

    , g (
2. From y-DATA-II.2: The evolutions of these estimates are shown in Figure 4.10, Figure 4.11 and Figure 4.12.

P (100) =     0.
Using the estimated parameters, we generate some simulations for y-DATA-II.1, y-DATA-II.2 shown in Figure 4.9a2) and Figure 4.9b2), respectively. Chapter 5

Some applications using our models

Résumé

Ce chapitre présente deux applications utilisant des trajectoires simulées par nos modèles : amélioration de l'estimation de la distribution de l'indice journalier de clarté et prédiction du rayonnement solaire total. Pour la première application, nous utilisons DTM-K estimé avec les données mesurées à l'île de La R éunion, et npEM un algorithme d'estimation de mélanges.

Pour la seconde application, nous utilisons CTM-k et CTM-y estimés par des données mesurées à 1Hz à l'île de La Guadeloupe et nous proposons une prédiction à un horizon d'une heure.

Abstract

This chapter presents two applications using paths simulated by our models: improvement of the daily clearness index distribution estimation and prediction of total solar radiation.

For the first application, we use DTM-K estimated from La R éunion island data and npEM a mixture estimation algorithm.

For the second application, we use CTM-k and CTM-y estimated from data sampled at 1Hz in Guadeloupe and we propose a prediction with horizon one hour.

Estimating the distribution of daily clearness index

Estimating the distribution of daily clearness index over a month or over a specific period can be of interest for deciding wether our model estimated over this period still works on a longer period or not. It can also be used for clustering daily CISs observed on various periods.

Indeed, using DTM-K with its parameter estimated from a sample of daily CISs, of one-month-length say, (Table 3.1), we can simulate a much larger n-sample of K h , say K ′ 1 , . . . , K ′ n , over this period and get a smooth estimation of the pdf over this month. Doing the same with another month and getting another n-sample, a KS (Kolmogorov-Smirnov) test can be performed to reject or not the hypothesis that both pdf are the same. If the hypothesis is rejected (w.r.t. a p-value), we can reject the hypothesis that both models are the same. On the other hand, KS distance between two sequences, computed from the two n-samples, can be used for clustering CISs by performing some standard clustering methods.

Kernel estimators

The Gaussian kernel estimator of the density is the function f δ defined as

f δ (x) = 1 δn n h=1 φ x -K ′ h δ , (5.1)
where δ > 0 is a bandwidth (a smoothing parameter) and φ(•) denotes the N (0, 1) density function kernel. This estimator is of course much smoother than the uniform kernel estimator, that is the empirical pdf f , defined as follows: divide [0, 1] interval (the range of K h ) into L sub-intervals ]x l-1 , x l ] of equal length ∆x = 1 L with x 0 = 0 and x l = l∆x, l = 1, 2, . . . , L, then

f (x) = n l n 1 ∆x , (5.2)
where n l is the number of observed values in the interval ]x l-1 , x l ], l = 1, 2, . . . , L.

Mixtures of nonparametric densities

The estimator of the preceding section can be improved by assuming that the density is a finite mixture of nonparametric densities, which means that the components of the mixture are completely unspecified. While such a mixture is identifiable and can be estimated by an EM-type algorithm in the parametric case (for example in the case of a Gaussian mixture, Figure 5.1), this model, in the nonparametric case, is not necessarily identifiable but can be estimated by the npEM (nonparametric EM) algorithm [Benaglia 2009] which is implemented in the npEM fuction of the R package mixtools.

Prediction

In this Section we implicitely assume that the model estimated on a time interval [0, T ] still works on the time interval [T, T ′ ], with T ′ > T . Using the model with its estimated parameters and the extraterrestrial radiation computed on [T, T ′ ], we generate a large number of paths in the period [T, T ′ ]. Prediction of total solar radiation on [T, T ′ ] then follows.

Confidence region and prediction error for hourly total solar radiation

First, CTM-k and CTM-y parameters have been estimated from the 1 Hz measurements of total solar radiation observed during T = 1h. Extraterrestrial radiation has been computed on [T, T ′ ] with T ′ = 1h. Using this, L = 500 paths have been simulated on [T, T ′ ], providing a confidence region of prediction for the total radiation.

The error of prediction is usually evaluated by the MSPE (Mean Square Prediction Error), that is the expectation (estimated by a mean) of the squared difference between predicted values and observed values. MSPE has also to be averaged over the L paths, that can also provide a confidence interval. We prefer here to simplfy the prediction evaluation by computing, for each simulated path, the relative error between the mean of predicted values and the mean of obseved values and then provide the mean of these relative errors.

Observed datasets are mentioned in Table 3.

Predictions and real value:

• Measurements: 07h-10h (DATA-07h10h), (k t ) average: m k = 0.76, (k t ) variance σ 2 k = 0.0014, 95% confidence interval for AVERAGE-10h11h: 95%C.I.07h10h = (942.4, 948.6) W/m 2 .

• Measurements: 09h-10h (DATA-09h10h), m k = 0.79, σ 2 k = 0.0001, 95%C.I.09h10h = (935.7, 942.6) W/m 2 .

• Forecasting period 10h-11h, Observed: m k = 0. Comments:

• AVERAGE-10h11h does not belong to 95%C.I.07h10h, the relative error between the average of simulated data and the average of observed data is 0.91%), probability distribution/histogram, average, variance as well as the disposition of classes of DATA-07h10h and that of observed data in 10h-11h are different (see the figure and the table above).

• In contrast, 95%C.I.09h10h contains AVERAGE-10h11h (relative error is 0.29%) and similar statistical characteristics are observed.

Case 2: Measurements: day 285, type III (cloudy day).

Predictions and observed values: (see table below)

• Measurements: DATA-07h10h (m k = 0.06, σ 2 k = 0.0015). 95% confidence interval for AVERAGE-10h11h: (118.3, 120.0) W/m 2 .

• Measurements: DATA-09h10h (m k = 0.11, σ 2 k = 0.0004). 95% confidence for AVERAGE-10h11h: (162.2, 163.8) W/m 2 .

• Observed: AVERAGE-10h11h 249.5 W/m 2 (m k = 0.22, σ 2 k = 0.0017). 5.16 and Figure 5.17) Our simulated data are distributed as the data used for estimating the models. So, if the meteorological conditions in the forecasting period change significantly (more clouds for instance), predictions (from CTM-y as well as CTM-k) do not work well.

Therfore we propose in our illustrations a short-horizon prediction interval [T, T ′ ], T ′ = T + 1 (hour) and this horizon is also of interest for solar energy providers.

Conclusion

In order to understand and model the behaviour of solar radiation and clearness index, besides the two other approaches commonly used ("physical modelling" and "statistical solar climatology", as mentioned in the introduction), we have proposed a new approach, a HMM-type stochastic model for taking in account the influence of meteorological regimes.

The parameters of the model are estimated from real data using filtering equations, Girsanov change of measure theorem in stochastic calculus and the celebrated EM algorithm. The method has been tested and illustrated on several different types of observational data. The simulated data has been used to estimate the distribution of the daily clearness index and to predict the total solar radiation in a short horizon.

As a conclusion, we emphasize some elements that compare our model to physical models and other statistic models. We then take some notes in the problem of model construction, the problem of parameter estimation as well as the problem of simulated data application. From this, we determine the work which can be continued to study in the future.

Comparison with physical models

Our models are mainly data-driven as parameters of models are estimated from observed data and the number of regimes is chosen by observing the distribution of those data but they also include a physical model part through the extraterrestrial radiation.

In the physical models, direct beam (I b ) and diffuse (I d ) radiation components are obtained as a function of the specific atmospheric transmittances. They require several physical parameters as inputs: water vapor absorption (T w ), Rayleigh scattering (T r ), uniformly mixed gases absorption (T g ), ozone absorption (T o ), aerosol total extinction (T a ), . . . For instance, Psiloglou et al. [Psiloglou 2000] proposed a clear-sky radiation model with the following components: However, our models also have a physical factor, namely the local extraterrestrial radiation I t . The observation equation (4.4) of CTM-y depends on I t . For DTM-K and CTM-k, the observations used in the parameter estimation are CISs computed using I t and measurements of total solar radiation G t :

(1.10) :

k t = G t I t ,
or (1.11) : K ∆t = ∆t G s ds ∆t I s ds .

As a consequen our models generate data having meteorological characteristics of areas where data used for estimation were collected.

Comparison with other statistical models

Statistical models and methods have played a wide range of applications in solar radiation and the climatology research.

As far as we know, most of statistical models used for solar radiation and clearness index hinge on regression models using the correlative relation between statistical variables (including main solar radiation components, clearness index and meteorological parameters such as sunshine duration, cloudiness, temperature, etc). The parameters of the statistical models are estimated from complete observation data of all these statistical variables.

For instance, Angstrom (1924) suggested the following linear expression for the relationship between daily clearness index K h and sunshine duration ratio S h [Tovar-Pescador 2008]:

K h = a + bS h ,
where a, b are model parameters which are estimated by regression technique from observation data of statistical variables K h and S h (complete data of model). This linear expression was used in practical applications for many years to estimate the daily, monthly and annual total solar radiation from the comparatively simple measurements of sunshine duration. [START_REF] Ogelmen | [END_REF][START_REF] Akinoglu | [END_REF]] proposed a quadratic expression by adding a non-linear term.

Our method is also based on several statistic techniques (ML parameter estimation, EM algorithm, filtered estimate, . . . ). However, in comparison with usual statistical models, our proposed models have three original features:

• Our models take in account the hidden (unobserved) dynamic of environmental regimes as well as observed variables (extraterrestrial, direct and diffuse radiations)

• Our model parameters can be estimated from incomplete data. This is important for solar radiation research as well as for climatological research because all climatological quantities and meteorological components are not always available.

• Our models having a random component, simulated data can be used to study the probability distribution of clearness index and total radiation.

Future works

Problems in the model construction

Our models depend on measurements of the total solar radiation and the extraterrestrial radiation I t with astronomical calculations depending on the geographical location. A better modelling should consider more information on observed data and meteorological parameters, dealing with multidimensional vector as it is done in many physical models.

Problems in parameter estimation

The technique we used in the problem of parameter estimation is the filtering estimates, also called forward estimates, which are based on the observation history up to time h in discrete time, Y K h σ{K 1 , K 2 , . . . , K h }, h = 1, 2, . . . (or up to time t in continuous time, Y Y t σ{Y s : 0 s t}, t ∈ [0, T ], respectively). To achieve better estimates, we think of backward estimates, an estimating technique based on observations in the future, Y K h:M σ{K h , K h+1 , . . . , K M } (or Y Y t:T σ{Y s : t s T }, resp.). Backward estimates are calculated as a backward recursion from the end of the batch of observations. The forward-backward estimate is termed smoothing estimation, based on the past, present and future of observation data, Y K h ∨ Y K h:M (or Y Y t ∨ Y Y t:T , resp.) [START_REF] Elliott | [END_REF][START_REF] James | [END_REF].

Problems in prediction

First, our prediction results should be compared to recent prediction results obtained by EDF (French electricity company) using some data analysis techniques. Next, our 1h prediction horizon was a short-term one and we need to refine our models to deal with very short-term (10mn) prediction and also long-term prediction (next day, next month, next year). Moreover, as seen in the examples of chapter 5, prediction is connected to classification problems (classification of sequences, of days, of several months periods) and this classification aspect is not considered in our present models. Conversely it can be thought that estimated parameters can be used for classification purposes.

Van Ly TRAN

Modèles Stochastiques des Processus de Rayonnement Solaire

Résumé :

Les caractéristiques des rayonnements solaire dépendent fortement de certains événements météorologiques non observés comme fréquence, taille et type des nuages et leurs propriétés optiques (aérosols atmosphériques, albédo du sol, vapeur d'eau, poussière et turbidité atmosphérique) tandis qu'une séquence du rayonnement solaire peut être observée et mesurée à une station donné. Ceci nous a suggéré de modéliser les processus de rayonnement solaire (ou d'indice de clarté) en utilisant un modèle Markovien caché (HMM), paire corrélée de processus stochastiques. Notre modèle principal est un HMM à temps continu (Xt, yt) t≥0 est tel que (yt), le processus observé de rayonnement, soit une solution de l'équation différentielle stochastique (EDS) : dyt = [g(Xt)Ityt]dt + σ(Xt)ytdWt, où It est le rayonnement extraterrestre à l'instant t, (Wt) est un mouvement Brownien standard et g(Xt), σ(Xt) sont des fonctions de la chaîne de Markov non observée (Xt) modélisant la dynamique des régimes environnementaux. Pour ajuster nos modèles aux données réelles observées, les procédures d'estimation utilisent l'algorithme EM et la méthode du changement de mesures par le théorème de Girsanov. Des équations de filtrage sont établies et les équations à temps continu sont approchées par des versions robustes. Les modèles ajustés sont appliqués à des fins de comparaison et classification de distributions et de prédiction. Mots de clé : rayonnement solaire, indice de clarté, HMM, EDS, algorithme EM, Thérème de Girsanov, filtrage.

Stochastic Models of Solar Radiation Processes

Abstract : Characteristics of solar radiation highly depend on some unobserved meteorological events such as frequency, height and type of the clouds and their optical properties (atmospheric aerosols, ground albedo, water vapor, dust and atmospheric turbidity) while a sequence of solar radiation can be observed and measured at a given station. This has suggested us to model solar radiation (or clearness index) processes using a hidden Markov model (HMM), a pair of correlated stochastic processes. Our main model is a continuous-time HMM (Xt, yt) t≥0 is such that the solar radiation process

Figure 1 :

 1 Figure 1: Measurements of total solar radiation and extraterrestrial radiation (a).Corresponding clearness index (b).[START_REF] Soubdhan | [END_REF]

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Extraterrestrial solar radiation . . . . . . . . . . . . . . . . . 2 1.2.1 Extraterrestrial normal radiation . . . . . . . . . . . . . . . . 2 1.2.2 Extraterrestrial horizontal radiation . . . . . . . . . . . . . . 4 1.3 Zenith angle calculation . . . . . . . . . . . . . . . . . . . . . 4 1.3.1 Equation of time . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3.2 Apparent solar time . . . . . . . . . . . . . . . . . . . . . . . 5 1.3.3 Hour angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3.4 Declination . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3.5 Zenith angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Total solar radiation . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4.1 Direct solar radiation . . . . . . . . . . . . . . . . . . . . . . 7 1.4.2 Diffuse solar radiation . . . . . . . . . . . . . . . . . . . . . . 7 1.5 Clearness index . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.6 Solar radiation measurement . . . . . . . . . . . . . . . . . . 8 1.6.1 Solar radiometers . . . . . . . . . . . . . . . . . . . . . . . . . 8

Figure 1 . 3 :

 13 Figure 1.3: Relative positions of the Sun (photo: [Tingilinde 2006]).

Figure 1 . 4 :

 14 Figure 1.4: The declination angles.

Figure 1

 1 Figure 1.5: G t = I b + I d .

  .6a). Pyranometer A Pyranometer is used to perfomed the horizontal total radiation G t or the diffuse radiarion I d . Pyranometers have a 180 0 field of view. The horizontal total radiation G t is measured by a Pyranometer with a horizontal sensor (Figure 1.6c) while the diffuse radiarion I d is measured by a shaded Pyranometer under a tracking ball (Figure 1.6b).

Figure 1

 1 Figure 1.6: Typical instruments for measuring solar radiation components: (a) Pyrheliometer, (b) shaded Pyranometer, (c) Pyranometer with a horizontal sensor.

Contents 2 . 1

 21 Conditional expectations . . . . . . . . . . . . . . . . . . . . . 2.1.1 Radon-Nikodym derivative . . . . . . . . . . . . . . . . . . . 2.1.2 Jensen inequality . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.3 Conditional Bayes formula . . . . . . . . . . . . . . . . . . . . 2.2 Martingale difference sequence . . . . . . . . . . . . . . . . . 2.3 Binary vector representation of a Markov chain . . . . . . . 2.4 Hidden Markov models . . . . . . . . . . . . . . . . . . . . . . 2.5 Discrete-time HMM . . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Filtrations, number of jumps, occupation time and level sums 2.5.2 Reference Probability Method of measure change . . . . . . . 2.5.3 Normalized and unnormalized filters . . . . . . . . . . . . . . 2.6 Some recalls on stochastic calculus . . . . . . . . . . . . . . . 2.6.1 Ito product rule . . . . . . . . . . . . . . . . . . . . . . . . . 2.6.2 Ito formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6.3 Girsanov theorem . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Continuous-time homogeneous Markov chain . . . . . . . . . 2.8 Continuous-time HMM . . . . . . . . . . . . . . . . . . . . . . 2.8.1 Filtrations, number of jumps, occupation time and level sums 2.8.2 Change of measure . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . 2.9.1 Likelihood function . . . . . . . . . . . . . . . . . . . . . . . . 2.9.2 Pseudo log-likelihood function . . . . . . . . . . . . . . . . . . 2.9.3 EM Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . .

Contents 3 . 1

 31 Modelling a daily clearness index sequence . . . . . . . . . . 3.1.1 State process . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Observation Process and model parameters . . . . . . . . . . 3.1.3 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . 3.1.4 Filtering equations . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Modelling a clearness index process on a time interval . . . 3.2.1 CTM-k model . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Change of measure . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . 3.2.4 Filtering equations . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Discrete-Time Approximate Model DTAM-k . . . . . . . . . 3.3.1 Components of DTAM-k . . . . . . . . . . . . . . . . . . . . . 3.3.2 Discrete-time approximate filtering equations . . . . . . . . . 3.3.3 Updating parameter . . . . . . . . . . . . . . . . . . . . . . . 3.4 Experiments with real data . . . . . . . . . . . . . . . . . . . 3.4.1 Real data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Estimations . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Figure 3.1: Histogram of observed data: (a) daily CIS measured in October 2010, La Réunion; (b) k-DATA-II.1 (see Table3.3).

Figure 3

 3 Figure 3.2: An illustration for the DTM-K: (X h , K h ) h=1,2,... .

m=1 t 0 XXXXX

 0 s , e i e ′ j e l -X s , e i e ′ j e m X s -, e m dX s , e l (e le m ) s , e i e ′ j e l -X s , e i e ′ j e m X s -, e m a lm ds(e le m ) s , e i e ′ j e l -X s , e i e ′ j e m X s -, e m dV s , e l (e le m ) s , e i e ′ j e l e i -X s , e i e ′ j e m e i , e m a lm ds(e le m ) s , e i e ′ j e l e i -X s , e i e ′ j e m e i , e m dV s , e l (e le m ).

m=1 t 0 XX

 0 s , e i e ′ j e l e i -X s , e i e ′ j e m e i , e m a lm ds(e le m ) s , e i e ′ j e l e i -X s , e i e ′ j e m e i , e m dV s , e l (e le m ).

  Figure 3.3, Figure 3.4 and Figure 3.5, we have decided that N = 4 states. 3.4.2.1 DTM-K parameter estimations DTM-K parameter are estimated when using three daily CISs observed in August 2009 (K-DATA-0809, variance 0.0080, Figure 3.6(a)), March 2010 (K-DATA-0310, variance 0.0484, Figure 3.7(a)) and October 2011 (K-DATA-1011, variance 0.0268, Figure 3.8(a)).

Figure 3 . 4 :Figure 3

 343 Figure 3.3: Histogram of daily CIS measured in: (a) August 2009, (b) March 2010, La Réunion.

Figure 3 . 6 :

 36 Figure 3.6: Graph of daily CIS measured in August 2011, La Réunion-France, (a), and its simulation, (b).

Figure 3 . 7 :

 37 Figure 3.7: Graph of daily CIS measured in March 2010, La Réunion-France, (a), and its simulation, (b).

  Figure 3.8: Graph of daily CIS measured in October 2011, La Réunion-France, (a), and its simulation, (b).

  = (0.0620, 0.4268, 0.5816, 0.6396) ′ , α (50) = (0.0425, 0.1097, 0.0250, 0.0084) ′ .The graphs in Figure3.9, Figure3.10 and Figure3.11 show the evolution of these estimates.Models with their estimated parameters were used to simulated some paths:Figure 3.6(b) and Figure 3.7(b).

Figure 3

 3 Figure 3.9: Estimation of probability transition matrix A from daily CIS measured in: (a) August 2009, (b) March 2010, La Réunion-France.

Figure 3

 3 Figure 3.10: Estimation of parameter vector b from daily CIS measured in: (a) August 2009, (b) March 2010, La Réunion-France.

Figure 3

 3 Figure 3.11: Estimation of parameter vector α from daily CIS measured in: (a) August 2009, (b) March 2010, La Réunion-France.

Figure 3 Figure 3 Figure 3

 333 Figure 3.12: Estimation of the transition matrix A from K-DATA-1011 and the initial parameter vector: (a) θ (1,0) ; (b) θ (2,0) .

  50) = (0.2189, -0.1635, 1.1028, -1.0794) ′ , β (50) = (0.0068, 0.0068, 0.0096, 0.0100) ′ . The graphs in Figure 3.16, Figure 3.17 and Figure 3.18 show the evolution of estimates from k-DATA-III.3 and k-DATA-II.2. The models with estimated parameters are used to simulate some paths: Figure 3.19(b), Figure 3.20(b) and Figure 3.21(b).

Figure 3 Figure 3

 33 Figure 3.16: Estimation of transition probability matrix P from: (a) k-DATA-III.3; (b) k-DATA-II.2.

Figure 3

 3 Figure 3.18: Estimation of parameter vector β from: (a) k-DATA-III.3; (b) k-DATA-II.2.

Figure 3

 3 Figure 3.19: Graph of observed data k-DATA-I.1, (a), and its simulation, (b).

Figure 3

 3 Figure 3.20: Graph of observed data k-DATA-III.3, (a), and its simulation, (b).

Figure 3

 3 Figure 3.21: Graph of observed data k-DATA-II.2, (a), and its simulation, (b).

Figure 4 . 1 :

 41 Figure 4.1: An illustration for the CTM-y: (X t , y t ) t∈[0,T ] .

  .39) completes the proof. Let us specify Theorem 4.1. Take H t = H 0 = 1, µ s = δ s = 0, ρ = 0. Then Corollary 4.1. The recursive equation for the state process p

Figure 4 .

 4 Figure 4.2: Histograms of observed data: (a) y-DATA-I.4; (b) y-DATA-III.4.

Figure 4 .

 4 Figure 4.3: Histograms of observed data: (a) y-DATA-II.1; (b) y-DATA-II.2.

Figure 4 . 4 (

 44 2 and Figure 4.3), we have decided that the number of states is N = 4. Consider the two data sequences: y-DATA-I.4 (Figure 4.7(a)), y-DATA-III.Figure 4.8(a)) and start with initial parameter vector θ (0) :

Figure 4 . 4 ,

 44 Figure 4.5 and Figure 4.6):

Figure 4 . 4 ,

 44 Figure 4.5 and Figure 4.6 show the estimates evolution for y-DATA-I.4 and y-DATA-III.4. Using the estimated parameters, we generate some simulations for y-DATA-I.4, y-DATA-III.4 which are shown in Figure 4.7(b) and Figure 4.8(b), respectively.

Figure 4 . 4 :

 44 Figure 4.4: Estimation of transition probability matrix P : (a) from y-DATA-I.4; (b) from y-DATA-III.4.

Figure 4 . 5 :

 45 Figure 4.5: Estimation of parameter vector g: (a) from y-DATA-I.4; (b) from y-DATA-III.4.

Figure 4 . 6 :

 46 Figure 4.6: Estimation of parameter vector σ: (a) from y-DATA-I.4; (b) from y-DATA-III.4.

Figure 4 .

 4 Figure 4.7: (a) Graph of observed data y-DATA-I.4; (b) A simulation for y-DATA-I.4.

Figure 4 .

 4 Figure 4.8: (a) Graph of observed data y-DATA-III.4; (b) A simulation for y-DATA-III.4.

Figure 4 .

 4 Figure 4.9: (a) Graph of observed data y-DATA-II.1, a1), and its simulation, a2); (b) Graph of observed data y-DATA-II.2, b1), and its simulation, b2).

Figure 4 .

 4 Figure 4.10: Estimation of transition probability matrix P : (a) from y-DATA-II.1; (b) from y-DATA-II.2.

Figure 4 .

 4 Figure 4.11: Estimation of parameter vector g: (a) from y-DATA-II.1; (b) from y-DATA-II.2.

Figure 4 .

 4 Figure 4.12: Estimation of parameter vector σ: (a) from y-DATA-II.1; (b) from y-DATA-II.2.

4. 6

 6 Simulations of total solar radiation day Applying similar methods, we have obtained the simulations shown in Figure 4.13(b), Figure 4.14(b) and Figure 4.15(b) for observed data y-DATA-day-I, y-DATA-day-II and y-DATA-day-III whose graphs are plotted in Figure 4.13(a), Figure 4.14(a) and Figure 4.15(a), respectively.

Figure 4 .

 4 Figure 4.13: (a) Graph of observed data y-DATA-day-I; (b) A simulation for y-DATA-day-I.

Figure 4 .

 4 Figure 4.14: (a) Graph of observed data y-DATA-day-II; (b) A simulation for y-DATA-day-II.

Figure 4 .

 4 Figure 4.15: (a) Graph of observed data y-DATA-day-III; (b) A simulation for y-DATA-day-III.

Figure 5

 5 Figure 5.1: pdf of La Réunion K h in March, estimated as a mixture by npEM algorithm

Figure 5

 5 Figure 5.6: pdf of K h for May in La Réunion: (a) observed values with ∆x = 1/200; (b) kernel estimation with δ = 1/200.

  3. The results are shown in Table 5.1. The results obtained from CTM-y seem better than that from CTM-k. Some cases of simulated paths are shown in Figure 5.7, Figure 5.8, Figure 5.9, Figure 5.11, Figure 5.13, Figure 5.15 and Figure 5.

  17

  .

Figure 5 Figure 5 Figure 5

 555 Figure 5.7: Paths simulated from 14h to 15h (a) using CTM-k with parameters estimated from k-DATA-I.2 (left), observed data (right), (b) using CTM-y with parameters estimated from y-DATA-I.2.

Figure 5 Figure 5 Figure 5

 555 Figure 5.10: (a) Simulated paths in 10h-11h generated by CTM-y with parameters estimated from observations in 07h-10h, day 262 (type I, clear day), 2006, Guadeloupe. (b) Data observed during 10h-11h: blue solid line.

Figure 5 Figure 5 Figure 5 Figure 5

 5555 Figure 5.14: (a) Simulated paths in 10h-11h generated by CTM-y with parameters estimated from observations in 07h-10h, day 118 (type II, partially cloudy day), 2006, Guadeloupe. (b) Observed data: blue solid line.

•

  total radiation : G t = I b + I d ,• direct-beam radiation :I b = I 0 cos θ z T w T r T g T a ,• diffuse radiation :I d = I 0 cos θ z T w T g T o T aa 1 -T a T aa T r /2 + I dm ,where T aa is the absorption aerosol broadband transmittance function, T dm is a multiple-scattering component, θ z is the zenith angle and I 0 is the extraterrestrial normal solar radiation in the n l -th day of the year:I 0 = I SC (1.00011 + .034221 cos Γ + .00128 sin Γ + .000719 cos 2Γ + .000077 sin 2Γ), here I SC = 1373W/m 2 is the solar constant and Γ (in radians) is the day angle which is represented by:

  (yt) t≥0 is a solution of the stochastic differential equation (SDE) : dyt = [g(Xt)Ityt]dt + σ(Xt)ytdWt, where It is the extraterrestrial radiation received at time t, (Wt) is a standard Brownian motion and g(Xt), σ(Xt) are functions of the unobserved Markov chain (Xt) modelling environmental regimes. To fit our models to observed real data, the estimation procedures combine the Expectation Maximization (EM) algorithm and the measure change method due to Girsanov theorem. Filtering equations are derived and continuoustime equations are approximated by robust versions. The models are applied to pdf comparison and classification and prediction purposes. Keywords : solar radiation, clearness index, HMM, SDE, EM algorithm, Girsanov theorem, filtration. MAPMO UMR 7349, Fédération Denis Poisson Université d'Orléans, UFR Science Bâtiment de mathématiques -Rue de Chartres B.P. 6759 -45067 Orléans cedex 2 FRANCE

Table 3

 3

	Data set DATA-I.1 DATA-I.2 DATA-I.3 DATA-II.1 DATA-II.2 DATA-II.3 DATA-III.1 DATA-III.2 DATA-III.3	Mean of CIS Variance of CIS time -dayth (in 2006) 0.7904 0.0069 09h10h -262th 0.7536 0.0062 13h14h -109th 0.7205 0.0097 13h14h -139th 0.6162 0.2715 09h10h -118th 0.4172 0.2234 13h14h -234th 0.6178 0.1730 09h10h -234th 0.2975 0.0308 13h14h -339th 0.1052 0.0201 09h10h -285th 0.1229 0.0269 13h14h -180th

.3: Mean clearness index (k) sequences computed from measurements of total solar radiation sample at 1-s in Guadeloupe-France in 2006. DATA-I.1 stands for sequences belonging to type I in Table

3

.2, and similarly for others.

  e i dy t = X t , e i [(g i I ty t )dt + σ i y t dW t ], (4.3) where I t is the extraterrestrial radiation received at time t, (W t ) is a standard Brownian motion and σ 2 i is a constant noise variance during regime i. , e i [(g i I ty t )dt + σ i y t dW t ].

	Markov chain X t ∈ {e 1 , e 2 , e 3 }
	← X t =e 3
	X t =e 2
	← X t =e 1 →
	We then have N i=1 and since N i=1 X t , e i = 1, (4.4) can be written as X t , e i dy t = N i=1 X t (4.4) 09:00 12:00 15:00 17:48 06:32 09:00 12:00 15:00 17:48 0 200 400 600 800 1000 time (hh:mm) y t (W/m 2 ) dy 06:32 1200 Observation process (y t )

t = [g(X t )I ty t ]dt + σ(X t )y t dW t , (4.5)

where σ(X t ) = N i=1 X t , e i σ i = X t , σ with σ = (σ 1 , σ 2 , . . . , σ N ) ′ .

  e i e i ds +

	Noticing that
	0	t	I s y s	σ i p s , e i e i	dy s σ(X s )y s	=	0	t	I s y s	σ i p s , e i e i	dy s σ i y s
	we then get
	Corollary 4.6.
										0	t	I s y s	σ i p s , e i e i	dy s σ(X s )y s	(4.46)

  are obtained from recursive equations (4.67), (4.68), (4.69), (4.70), respectively, and the processes γ

  Contents 5.1 Estimating the experimental distribution of K h . . . . . . . 5.1.1 Kernel estimators . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.2 Mixtures of nonparametric densities . . . . . . . . . . . . . . 5.1.3 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Prediction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.2.1 Confidence region and prediction error for hourly total solar radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Discussion on the prediction results . . . . . . . . . . . . . .

Histogram of k t of uesd DATA: from 07h to 10h

  80, σ 2 k = 0.0001, Average total radiation: 936.9 W/m 2 .

	0.04
	0.02
	0.64 0	0.66	0.68	0.7	0.72	0.74	0.76	0.78	0.8	0.82	0.84
	0.04
	0.02
	0.775 0	0.78		0.785	0.79		0.795	0.8		0.805	0.81

Histogram of k t of uesd DATA: from 09h to 10h

	0.04
	0.02
	0.79 0	0.795	0.8	0.805	0.81	0.815	0.82	0.825

Histogram of k t of observed DATA: from 10h to 11h

		k t		95% C.I./	err.
	Time	m k	σ 2 k	Mean obs. data (%)
	07h-10h .76 .0014	(942.4,948.6)	0.91
	09h-10h .79 .0001	(935.7,942.6)	0.29
	10h-11h .80 .0001	936.9 W/m 2	0.00
	(simulated paths are plotted in Figure 5.10 and
			Figure 5.11)

Histogram of k t of uesd DATA: from 07h to 10h

	0.04
	0.02
	0	0	0.02	0.04	0.06	0.08	0.1	0.12	0.14	0.16
	0.1
	0.05
	0.04 0	0.06		0.08	0.1	0.12		0.14	0.16

Histogram of k t of uesd DATA: from 09h to 10h

	0.1
	0.05
	0	0.16	0.18	0.2	0.22	0.24	0.26	0.28	0.3	0.32

Histogram of k t of observed DATA: from 10h to 11h

  52.3%) neither to 95%C.I.09h10h (error: 34.7%). Distribution characteristics of the data used for estimations are not similar to that of data in the prediction interval 10h-11h. Measurements: Day 118, type II (partially cloudy day). Distribution characteristics of the data used for estimations (07h-10h, 09h10h) are similar to those of data observed in the prediction interval 10h-11h. The obtained predictions are quite nice (see the figure and the table below).

								k t	95% C.I./	err.
							Time	m k	σ 2 k	Mean obs. data (%)
							07h-10h .06 .0015	(118.3,120.0)	52.3
							09h-10h .11 .0004	(162.2,163.8)	34.7
							10h-11h .22 .0017	249.5 W/m 2	0.00
								Simulated paths are plotted
								in Figure 5.12 and Figure 5.13
	Comments:		AVERAGE-10h11h does not belong to 95%C.I.07h10h
	0.2 (error: Case 3: 0.1 0.3 0.4 0 0.02 0.04	0.5	0.6	0.7	0.8	0.9	1

Histogram of k t of uesd DATA: from 07h to 10h

	0.04
	0.02
	0.1 0	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1

Histogram of k t of uesd DATA: from 09h to 10h

	0.1
	0.05
	0	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1

Histogram of k t of observed DATA: from 10h to 11h

  Case 4: Measurements Day 234, type II.The distribution characteristics of data used for estimation (07h-10, 09h-10h) are not similar to that of data in the prediction interval 10h-11h. Predictions present large errors (see table and figure below).

	0.04
	0.02
	0.2 0	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
										k t	95% C.I./	err.
									Time	m k	σ 2 k	Mean obs. data (%)
									07h-10h .63 .0643	(787.4,806.2)	0.06
									09h-10h .62 .0472	(780.5,798.7)	0.75
									10h-11h .61 .0627	797.2 W/m 2	0.00
										(Simulated paths are plotted
										in Figure 5.14 and Figure 5.15)

Histogram of k t of uesd DATA: from 07h to 10h

	0.04
	0.02
	0.2 0	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1

Histogram of k t of uesd DATA: from 09h to 10h

	0.1
	0.05
	0.2 0	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1

Histogram of k t of observed DATA: from 10h to 11h

		k t		95% C.I./	err.
	Time	m k	σ 2 k	Mean obs. data (%)
	07h-10h .54 .0397	(656.7,671.5)	36.0
	09h-10h .61 .0471	(646.3,661.9)	33.9
	10h-11h .42 .0499	488.5 W/m 2	0.00
		(Simulated paths are plotted
		in Figure

Remerciements

Experiments with real data

The parameters will be estimated from the observed data, the number of states being chosen after examining the data histograms. We deal with data coming from some tropical humid areas, but our method can also be tested on other types of climate.

Real data

Using standard formulas of the extraterrestrial radiation (Chapter 1), we compute the CISs from the total solar radiation measurements performed in two French islands: Guadeloupe and La Réunion which are tropical areas (as mentioned in Section 1.6.2).

To apply DTM-K model, we have used one month length samples of daily clearness index (Table 3 We have applied DTAM-k model of CTM-k with total radiations sampled at 1-s (Section 1.6.2) in Guadeloupe in 2006 (see Table 3.3) by Dr. T. Soubdhan, University of Antilles-Guyane. Time interval [0, T ] with T = 1 hour was split using (3.65) with a small enough ∆ = 1/3600.

Mean

Discussion on the prediction results

Our prediction method by simulating a large number of paths hinges on the strong law of large numbers. Predictions quality highly depends on the pdf, average and variance of the data observed during the learning time interval [0, T ] and used to estimate the parameters of the model. As it will be observed, prediction quality also depends on the number of states N of the model and on the distribution of the forecasting period data.

We illustrate this considering 4 examples of data recorded in Guadeloupe in 2006 and using CTM-y model. Predictions on 10h-11h are done from data observed in the interval 07h-10h (resp. 09h-10h).