Benoît Jean-François Scheid

Madame Sabrina

Madame Elodie

Anis Amri

Benjamin Alvarez

Clèmence Karmann

Hassan Mohsen

Mattieu Brachet

Nassim Sahki

Exopsé sur les " Longitudinal multicategorical processes : Generalized

Keywords: prediction, ⇓: filter, ←: smoothing, →: increasing horizon prediction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

We consider state space models where the observations are multicategorical and longitudinal, and the state is described by CHARN models. We estimate the state by generalized Kalman recursions, which rely on a variety of particle filters and EM algorithm. Our results are applied to estimating the latent trait in quality of life, and this furnishes an alternative and a generalization of existing methods. These results are illustrated by numerical simulations and an application to real data in the quality of life of patients surged for breast cancer. -Generalized Kalman recursions -Generalized state space models -Multicategorical longitudinal data -Latent variables -Particle filters -EM algorithm vii Je tiens tout d'abord à remercier le directeur de cette thèse, M. Joseph NGATCHOU-WANDJI, pour m'avoir fait confiance malgré les connaissances plutôt légères que j'avais en octobre 2013 sur le modèle espace-état, puis pour m'avoir guidé, encouragé, conseillé, et pour avoir corrigé erreurs scientifiques et fautes d'orthographe pendant presque quatre ans. Je remercie Messieurs Didier MAQUIN (directeur de l'école doctorale), Xavier AN-TOINE (directeur du laboratoire), Olivier Garet (directeur du département de mathématiques) pour la confiance qu'ils m'ont accordée. Je voudrais remercier les rapporteurs

Résumé

Les modèles espace d'état, aussi connus sous le nom de modèles dynamiques, relient des séries chronologiques observées à des séries chronologiques non observées, par un système de deux équations dont l'une décrit en général une relation entre les termes de la série observée et ceux de la série non observée (équation d'observation), et l'autre une relation autorégressive entre les valeurs de la série non observée (équation d'état). Les valeurs de la série non observée sont appelées états et celles de la série observée observations. Nous considérons dans cette thèse, des modèles espace d'état très généraux dans lesquels les fonctions associées peuvent être non-linéaires et les bruits des modèles non-gaussiens. Nous estimons les états de ces modèles en utilisant les récursivités de Kalman généralisées, les filtres particulaires et l'algorithme EM. Notre travail est motivé par l'envie d'estimer le trait latent en qualité de vie, et d'autres types de variables latentes en économie ou dans le monde industriel. Des exemples concrets de variables latentes sont : la santé des patients, la confiance dans les entreprises, le moral des clients d'une firme, le niveau d'anxiété d'utilisateurs de machines ou de robots dans les usines. Ici, les questionnaires sont considérés comme les observations et les variables latentes comme les états. Notre approche est à la fois une alternative aux travaux existant dans la littérature et leur généralisation. Plus précisément, dans notre travail, nous nous intéressons aux variables latentes X i (t) produites par un individu i, (i = 1, • • • , n), au temps t, (t = 1, • • • , T ). Les X i (t) peuvent être la santé du patient, un trait latent, etc. Nous observons seulement Y i (t) au lieu de X i (t). Les Y i (t) sont les réponses des individus aux questionnaires. Conformément aux objectifs et à la portée de l'étude, le contenu de cette thèse est structuré en six chapitres, que nous résumons ci-dessous.

Chapitre un : Dans ce chapitre nous faisons un survol de la littérature. Nous y présentons notamment quelques-uns des travaux où l'équation d'état est un modèle linéaire.

Chapitre deux : Ce chapitre présente brièvement les principaux outils mathématiques utilisés. Premièrement, nous rappelons la définition du modèle espace-état donnée par [START_REF] Fahrmeir | Multivariate statistical modelling based on generalized linear models[END_REF]. Aussi connu sous le nom de modèle dynamique, il relie des observations de séries chronologiques ou données longitudinales {Y t } à des "états" non i observés {X t } par un modèle d'observations donnant {Y t } sachant {X t }. Les états sont supposés suivre un modèle de transition stochastique. Dans ce chapitre, nous présentons différents modèles espace-état : linéaires avec des bruits gaussiens, nonlinéaires avec bruits gaussiens, et non-linéaires avec bruits non-gaussiens. Les récursivités de Kalman sont des algorithmes qui utilisent une série de mesures (variables) observées au fil du temps pour produire des estimations des variables latentes de trois manières différentes : prédiction pour t > T, filtrage pour t = T, et lissage pour t < T. Dans les modèles espace-état linéaire ou non-linéaire, si les bruits sont gaussiens, la distribution à posteriori est gaussienne. Par conséquent, les récursivités de Kalman linéaires sont utilisées pour estimer la moyenne à posteriori et la matrice de variancecovariance. Dans le modèle espace-état non-linéaires avec bruits gaussiens, la distribution à posteriori est encore approximativement gaussienne.

Nous utilisons dans ce cas le développement de Taylor pour linéariser les fonctions non-linéaires impliquées dans les modèles étudiés. L'algorithme obtenu s'appelle "récursivités de Kalman étendues". Le modèle espace-état non-linéaire et non-gaussien applique la perspective bayésienne en utilisant les densités conditionnelles. Les filtres particulaires permettent d'approximer les distributions a posteriori permettant d'estimer les valeurs latentes par la prédiction, le filtrage et le lissage. Plusieurs algorithmes de filtres particulaires ont été proposés. Nous présentons les plus connus.

1. L'algorithme du filtre particulaire de [START_REF] Kitagawa | Monte Carlo filter and smoother for non-Gaussian nonlinear state space models[END_REF]. L'idée fondamentale est d'utiliser un grand nombre d'échantillons (particules) auxquels nous attribuons des poids, pour approximer la probabilité a posteriori des états.

2. L'algorithme des filtres particulaires auxiliaires : L'algorithme précédent présente un problème de dégénérescence. Pitt et Shephard (1999) ont proposé "l'algorithme des filtres particulaires auxiliaires" pour résoudre ce problème. L'idée de base est d'introduire une variable auxiliaire ς dans la distribution d'importance pour générer les particules.

3. L'algorithme du filtre particulaire auxiliaire de Kalman étendu itéré.

: i = 1, • • • , n, k = 1, • • • , q, t = 1, • • • , T
), où n est le nombre d'individus, q le nombre d'items et t un instant. À chaque instant, nous observons des réponses binaires au lieu de variables latentes X i (t). Les probabilités conditionnelles sont données par le modèle de Rasch qui est largement utilisé dans divers domaines psychométriques comme la recherche en éducation ou l'analyse de la santé, et en particulier en qualité de vie. Ils ont considéré une étude longitudinale, où les patients répondent aux questions d'un formulaire à des dates régulières de visite, afin que soit mesurée leur santé du moment. Il est reconnu que la santé est un concept multidimensionnel latent. En pratique, chaque dimension est généralement évaluée par une ou plusieurs questions. Les auteurs se sont concentrés sur un cas particulier d'options de réponse dichotomique pour chaque question (oui, non, d'accord, en désaccord, etc.).

Chapitre trois : Dans ce chapitre, nous présentons une nouvelle classe de processus longitudinaux multivariés multicatégorielles. Les données proviennent d'une étude longitudinale, où les patients participent à une entrevue. L'entrevue vise à mesurer la santé des patients à intervalles réguliers, dont les dates sont déterminées avant l'étude. Il s'agit souvent de remplir un questionnaire avec des questions à choix multiples. Ce type questionnaire est utilisé pour des études en qualité de vie pour estimer le bien être des patients, en économie pour estimer la confiance des entreprises ou le moral des clients, et enfin dans le domaine industriel pour estimer le niveau d'anxiété des utilisateurs de machines ou de robots dans les usines. Notre approche est une alternative aux travaux de Bousseboua et Mesbah (2010). C'est aussi une généralisation de ceux-ci, ainsi que de ceux de [START_REF] Bartolucci | Longitudinal analysis of selfreported health status by mixture latent auto-regressive models[END_REF], [START_REF] Bartolucci | Modeling Longitudinal Data by Latent Markov Models with Application to Educational and Psychological Measurement[END_REF], Fahrmeir et Wagenpfiel (1997) et [START_REF] Fahrmeir | Multivariate statistical modelling based on generalized linear models[END_REF]. Comme nous l'avons déjà dit, nous nous intéressons aux variables latentes X i (t) produites par un individu i, (i = 1, • • • , n), au temps t, (t = 1, • • • , T ). Nous observons seulement Y i (t) au lieu de X i (t). Dans les travaux précédents, les variables latentes X i (t) sont décrites par le modèle autorégressif du premier ordre AR(1) avec un bruit gaussien. Dans ce chapitre, ils sont décrits par des modèles conditionnellement hétéroscédastiques non-linéaires (CHARN) du premier ordre avec un bruit gaussien. Ici nous faisons l'hypothèse que les observations proviennent d'une distribution multinomiale et que les variables latentes sont gaussiennes. Puisque la distribution a posteriori n'est pas symétrique, le mode a posteriori est une estimation ponctuelle de la variable latente. Dans ce chapitre, nous présentons deux approches pour estimer le mode postérieur. La première est basée sur l'utilisation de la récursivité de Kalman étendue. Dans l'appendice A.2, nous rapellons l'équation du filtre de Kalman étendu. La deuxième approche est basée sur l'utilisation du"maximum a posteriori (MAP)" issu de l'algorithme du filtre particulaire auxiliaire de Kalman étendu itéré. Les paramètres des modèles sont estimés par la méthode du maximum de vraisemblance via l'algorithme EM. Nous rappellons que l'algorithme EM est un algorithme itératif qui génère une séquence d'estimations du paramètre étudié. Chaque itération se décompose en deux étapes. La première étape E "Expectation", c'est-à-dire "Espérance" calcule une vraisemblance à partir de la formule d'un modèle espace-état. La seconde étape M "Maximisation" consiste à rechercher un jeu de paramètres maximisant la vraisemblance estimée à l'étape E. Nous calculons la matrice de variance-covariance des paramètres en exécutant l'"Identité d [START_REF] Oakes | Direct calculation of the information matrix via the EM[END_REF]. L'estimation du vecteur des paramètres de notre modèle est basé sur la distribution a posteriori p(X i | Y i ). dont la densité peut être calculée par une approche bayésienne. Dans le modèle état-espace linéaire gaussien, la distribution posteriori est gaussienne. Par conséquent, la moyenne conditionnelle et la matrice de variance-covariance du vecteur d'état estimé sont calculées par les récursivités du filtre de Kalman. Avec le modèle état-espace généralisé non-gaussien comme celui que nous étudions, la distribution de X i (t) sachant Y i (t) est généralement non-gaussienne. Par conséquent, il est nécessaire d'utiliser les méthodes des filtres particulaires pour trouver l'approximation de la distribution a posteriori. L'algorithme que nous utilisons ici est celui du filtre particulaire auxiliaire de Kalman étendu itéré proposé par [START_REF] Xi | The auxiliary iterated extended Kalman particle filter[END_REF]. Nous développons cet algorithme avec notre modèle et nous obtenons les équations du mode a posteriori et celle de la covariance a posteriori. La méthode utilisée est présentée à l'appendice A.2.

Chapitre quatre : Dans ce chapitre, nous généralisons notre modèle au cas où le bruit du modèle d'état est supposé issu d'une famille exponentielle. En pratique, les lois parmi les plus courantes de ces familles sont la loi normale, la loi exponentielle, la loi Gamma, la loi du khi-deux, la loi Beta, la loi de Dirichlet, la loi de Bernoulli, la loi multinomiale, la loi de Poisson, la loi de Wishart, la loi inverse Wishart et plusieurs autres. Comme au chapitre 3, nous trouvons aussi le mode a posteriori par deux approches : celle basée sur l'utilisation de la récursivité de Kalman étendue puis celle basée sur le "maximum a posteriori (MAP)" calculé par l'algorithme du filtre particulaire auxiliaire de Kalman étendu itéré. Comme précédemment, les paramètres des modèles sont estimés par la méthode du maximum de vraisemblance par l'algorithme EM. Nous calculons aussi la matrice de variance-covariance des estimateurs en utilisant l'"Identité d [START_REF] Oakes | Direct calculation of the information matrix via the EM[END_REF]. La loi a posteriori p(X i | Y i ) est approchée par l'algorithme du filtre particulaire auxiliaire de Kalman étendu itéré.

Chapitre cinq : Dans ce chapitre, nous présentons et commentons les résultats des simulations numériques que nous avons effectuées et l'application de nos résultats aux données réelles de santé. Il y a deux parties. La première, concerne les simulations numériques où les données longitudinales sont engendrées avec des variables latentes décrites par un modèle CHARN dont le bruit est de loi appartenant à une famille exponentielle. Dans la deuxième partie, nous appliquons nos résultats aux données réelles sur la qualité de vie des femmes ayant subi une opération pour cause de cancer du sein. Dans les deux parties, des routines R sont créées à partir des méthodes décrites dans les chapitres 3 et 4. Notre objectif dans ce chapitre est d'estimer les variables latentes par mode a posteriori via les récursivités du filtre de Kalman étendu.

Application aux données réelles :

Rotonda et al. (2011) ont conduit une étude sur la qualité de vie. Elle porte sur les facteurs corrélés à la fatigue liée au cancer du sein, chez les femmes ayant subi une opération pour cause de cette maladie. Des patients en nombre 502 ont été recrutés de septembre 2008 à septembre 2010 dans trois centres cancéreux français : le centre anticancéreux Alexis Vautrin de Lorraine, le centre anticancéreux Georges-Francois Leclerc de Bourgogne et le centre anticancéreux Paul Strauss d'Alsace, France. Ils ont rempli le questionnaire plusieurs fois, soit lors de leurs visites en clinique, soit à leur domicile après avoir reçu une enveloppe affranchie pour retourner leurs réponses. Le questionnaire utilisé est "State-Trait Anxiety Inventory" dont la version française (STAI-B) comporte vingt items. Pour chaque item, les réponses des patients sont classées en 4 catégories (non, plutôt non, plutôt oui, oui). Ici, la variable latente est la fatigue du patient après la chirurgie. Cette variable est supposée être quantitative et varie dans le temps autour d'une valeur moyenne supposée nulle dans notre travail. Nous estimons la fatigue du patient par le mode a posteriori en utilisant deux modèles pour la variable latente : les modèles AR [START_REF] Agresti | Categorical Data Analysis[END_REF] et CHARN [START_REF] Agresti | Categorical Data Analysis[END_REF][START_REF] Agresti | Categorical Data Analysis[END_REF].

Chapitre six : Ce chapitre présente les contributions de l'étude, ses conclusions et quelques recommandations pour les études futures. v 1. Conclusion : Nous avons appliqué notre approche sur différentes études de simulation, qui ont été conçues avec les différents types de questionnaires et les différentes formes d'équations d'état. Nous l'avons aussi appliquée directement aux données réelles sur la qualité de vie des femmes ayant subi une opération pour le cancer du sein. Notre approche a fourni une estimation de la variable latente à partir de données sur la qualité de vie. Ce qui permet donc d'avoir une idée plus précise de l'intensité, à chaqu'instant, de la fatique de chaqu'une des patientes ayant participé à l'expérience. Ce qui peut aussi permettre une meilleure prise en charge et même dans une certaine mesure de faire de la prévention dans le traitement du cancer du sein.

Perspectives :

• Notre approche, qui permet d'estimer une variable latente dans le domaine de la santé, peut aussi s'appliquer dans des domaines tels que le domaine économique, pour estimer la confiance des entreprises ou le moral des clients, le domaine industriel, pour estimer le niveau d'anxiété des utilisateurs de machines ou de robots en usine.

• Dans le chapitre 5, nous avons utilisé la première approche des récursivités de Kalman pour la création de routines R. Mais celles-ci peuvent aussi être créées pour la méthode du "maximum a posteriori (MAP)" basée sur le filtre particulaire auxiliaire de Kalman étendu itéré. On pourrait comparer les deux approches.

• Dans l'analyse des données réelles, nous avons constaté que les données étaient incomplètes : il y manquait des covariables ou des réponses des individus. Analyser ces données fut difficile. Nous avons donc choisi d'ignorer les individus pour lesquels il manquait des données. Dans nos travaux à venir, nous allons étudier la possibilité d'étendre notre approche au cas où les données sont incomplètes.

• Dans notre étude, nous avons supposé que la fonction de lien est un prédicteur linéaire η s ik (t) = u i (t)β s k + X i (t). Mais Hastie et [START_REF] Tibshirani | Generalized additive models[END_REF] ont introduit une classe des modèles additifs généralisés qui remplace la forme linéaire par une somme de fonctions lisses j s j (u j ), où les sj(.) sont des fonctions non spécifiées. De la même manière nous pouvons développer notre modèle avec η s ik (t) = s i (u i (t)) + X i (t).

Dédicaces

• A la mémoire de mon père qui me dit, peu avant sa mort: "mon fils, je veux te voir professeur".

• A la mémoire de mon frère, décédé le 23/01/2015 alors que j'étais en France. Je me souviens encore de ses paroles lorsqu'il m'a appelé la nuit avant sa mort:"Tu me manques tellement, je voudrais te serrer dans mes bras".

• A ma respectueuse mère, qui a prié pour moi tout au long de mes études doctorales.

• A ma chère épouse pour sa contribution, sa patience et sa compréhension tout au long de mes études doctorales. Elle m' a incroyablement soutenu et sans elle, rien de tout cela n'aurait été possible pour moi.

• Pour mes filles, Zahraa, Tuqa, Jannah et Safa qui n'ont cessé de m'accompagner au travail par leur amour, qui a toujours été ma plus grande source d'inspiration.

Chapter 1 Introduction

Introduction

State space models, also known as dynamic models, relate time series observations or longitudinal data {Y t } to unobserved "states" {X t } by an observation model for {Y t } given {X t }. The states are assumed to follow a stochastic transition model [START_REF] Fahrmeir | Multivariate statistical modelling based on generalized linear models[END_REF]). We consider generalized state space models where the functions involved may be nonlinear and the noises non-gaussian. We aim at estimating the state observations by using generalized Kalman recursions, which involve the use of particle filters and the EM algorithm. Our work is motivated by the desire to estimate the latent trait in quality of life, where we consider questionnaires as the observations and the latent variables as the states.

Our approach is an alternative to [START_REF] Bousseboua | Processus de markov longitudinal latent rasch observable[END_REF], but derived in a more general context. It is also a generalization of [START_REF] Bartolucci | Longitudinal analysis of selfreported health status by mixture latent auto-regressive models[END_REF], Bartolucci (2014), Fahrmeir and Wagenpfiel (1997) and [START_REF] Fahrmeir | Multivariate statistical modelling based on generalized linear models[END_REF]. [START_REF] Bousseboua | Processus de markov longitudinal latent rasch observable[END_REF] give a natural extension of the Rasch model to longitudinal data, but they studied a special case of a dichotomous response option for each question (yes-no, agree-disagree, etc.) They defined their model and estimated the parameters by maximizing the likelihood function via E-M algorithm. They illustrated their model with two classes of latent processes: i) the general Markov latent processes, ii) the auto-regressive AR(1) latent processes. They assumed the noise in the state equation has a gaussian distribution. They determined the distribution of latent trait and estimated the parameters of their models. The monthly questionnaire contained questions on the tendency of realization and expectations of variables like production, orders in hand, and demand. Answers were categorical, most of them trichotomous with categories like "increase", "decrease", and "no change". Thus, for each firm, the data form a categorical time series. Considering all firms within a specific branch, they had a categorical panel or longitudinal data. [START_REF] Durbin | Time series analysis of non-Gaussian observations based on state space models from both classical and Bayesian perspectives[END_REF] discussed the non-gaussian state space models from both classical and Bayesian perspectives with real non-gaussian time series data. They proposed an approach based entirely on importance sampling and antithetic variables.

They gave advantage to their approach: first, they entirely avoided the convergence problems that were associated with MCMC algorithms. Second, they easily computed error variances due to simulation as a routine part of the analysis, where they confirmed that the investigator could attain any predetermined level of simulation accuracy by increasing the sample size, where necessary, by a specific amount. . Almost all these above cited works considered linear models for the latent variable. In our work, we considered a non-linear model for this variable. We first assume the noise in state equation from gaussian distribution. We develope the work of Fahrmeir and Wagenpfiel (1997) where the Extended Kalman filter and smoother is used.

Our work provides a new class of multicategorical longitudinal multivariate processes of observed outcomes. Its context is that of a longitudinal study, where participant patients were interviewed about their health at regular intervals. The dates of which established before the study. The interviews typically involved filling out a questionnaire in which they were asked multiple choice questions. These questionnaires were constructed to measure the patient's perceived health at the time of the visit. It is recognized that health is a latent multidimensional concept.

More precisely, we are interested in latent variables X i (t) produced by an individual i,

(i = 1, • • • , n), at a time t, (t = 1, • • • , T ).
The X i (t) s might be the patients health, a latent trait. Our work can be applied in health field for studying the patients health.

It can be applied in economics field for studying the business confidence or the morale of customers. It can be applied in industrial field for the study of the level of anxiety of workers from machines or robots in factories. We only observe Y i (t) instead of X i (t). The Y i (t) s are the responses of the individuals to the questionnaire. The latent processes are described by Conditional Heteroscedastic Autoregressive Nonlinear (CHARN) models in two cases: i) gaussian noise processes ii) noise processes from exponential family distributions. We estimate the parameters of the model by the EM algorithm, after we use the Auxiliary Iterated Extended Kalman Particle Filter( AIEKPF) algorithm for determining the posterior likelihood. Next, the latent variables are estimated by the penalized likelihood and working extended Kalman filtering recursions and the maximum a posteriori (MAP) approach.

Overview of the dissertation :

In conformity with the objectives and the scope of the study, the contents of this dissertation is structured in six chapters. The dissertation chapters are organized so that the study objectives are apparent and conducted in the sequence outline.

Chapter Two : This chapter briefly presents the main mathematical tools used: State space models: linear, non-linear, gaussian, and non-gaussian. Particle filter methods. The Rasch models. The EM algorithm.

Chapter Three : This chapter presents a new class of longitudinal data through the state space models, where the observations are multicategoricals and the state variables are described by CHARN model with gaussian noise process. We discuss the estimation of the model parameters by the maximum likelihood via the EM algorithm, and the consistency and asymptotic normality to these estimators. Estimation of the parameters needed the posterior distribution which is found by using the Auxiliary Iterated Extended Kalman filter(AIEKPF). The optimality of the particles is presented. In this chapter, the objective is achieved through estimating the latent trait, which is estimated by the posterior mode via two methods. First, the penalized likelihood and working extended Kalman filtering. Second, the maximum a posteriori estimation (MAP).

Chapter Four : The objective of this chapter is to generalize our approach of chapter 3 to the case of noise processes from exponential family distribution. As in chapter 3, we find the estimation of the latent trait and model parameters. The consistency and asymptotic normality of the parameters is discussed. We present the algorithms to some of the methods explained in chapter 3 to avoid the prolongation .

Chapter Five : In this chapter, the numerical results are carried out to estimating latent variables. This chapter contains two parts. In the first, we generate longitudinal data with the latent trait from exponential family distributions. In the second, the real data are from a longitudinal study in women aged 18 years and older, who were surged for breast cancer in France. R-packages were created for the algorithms were mentioned in chapters 3.

Chapter Six : This chapter provides with the summary and detailed discussions of the dissertation conclusions. Areas for future research are also invoked.

Chapter 2

Mathematical tools 2.1 Introduction

This chapter focuses on the theoretical side of the state space models through the linear and nonlinear state space models. Kalman Filter and Smoother recursions are widely used to estimate the state variables and their variances. Therefore, it is presented for linear and nonlinear state space models. Two applications are given. First, the exponential family state space models discussed by Fahrmeir Wagenfiel (1997). Second, the Longitudinal Rasch process studied by Bousseboua and Mesbah (2010).

State space representations

X t = F t X t-1 + H t ε t , t = 1, 2, • • • (2.2.2)
where

{ε t , t = 1, 2, • • • } is the state noise, that is, a k-dimensional white noise. ε t for all t = 1, 2, • • • , ε t ∼ N (0; R t ), R t is a k × k matrix variance-covarince of state noise.
F t is k × k a matrix called state matrix, and H t is k × m a matrix called input matrix. {ε t } is uncorrelated with {ξ t } (i.e., E(ξ t ε s ) = 0 for all s and t). To complete the specification, it is assumed that the initial state X 0 is uncorrelated with all of the noise terms {ε t } and {ξ t }, and X 0 ∼ N (x 0 ; R 0 ). Basically, the system matrices G t , F t , Σ t , R t and x 0 , R 0 (the mean and variance for initial state) are assumed to be deterministic and known.

Non-linear state space models

The observation process (Y t ), Y t ∈ R , is described by an equation of the form

Y t = G(u t , X t , λ) + ξ t . (2.2.3)
where G(.) is a non-linear function, ξ t is the observation noise, λ is a parameter, (u t ) ∈ R r is the covariate process and (X t ) is the state process described by an equation of the form

X t = F (u t , X t-1 , γ) + H(u t , X t-1 , δ) • ε t , (2.2.4)
where F (.) and H(.) are non-linear functions, (ε t ) is the state noise, and γ is a parameter.

Gaussian State space models make the assumption that:

• {ξ t , t = 1, 2, • • • } is assumed to have a gaussian distribution: ξ t ∼ N (0; Σ t ), Σ t
the variance-covariance matrix of the observations.

• {ε t , t = 1, 2, • • • } has also a gaussian distribution for all, t = 1, 2, • • • , ε t ∼ N (0; R t ), R t ∈ R k
, where R t is the matrix variance-covariance of state noise.

• (ξ t ) and (ε t ) are uncorrelated that is, for all t, j

E[ξ t ] = 0 , E[ξ t ξ t ] = Σ t , E[ξ t ξ j ] = 0 for t = j. E[ε t ] = 0, E[ε t ε t ] = R t , E[ε t ε j ] = 0 for t = j, E[ξ t ε j ]
= 0 for all t and j.

Since this state space model is nonlinear and gaussian, the posterior distribution p(X t | Y t ) is also approximately gaussian. For that, we use the extended Kalman filter and smoother to compute the posterior mode and covariance matrix.

Non-gaussian state space models make the assumption that, (ξ t ) and (ε t ) have density functions r(ξ) and q(ε) respectively both of known forms. We can easily check that

E[X t | X t-1 , u t ] = F [u t , X t-1 , γ] Var[X t | X t-1 , u t ] = H 2 [u t , X t-1 , δ]R t .
Non-linear state space models are called Generalized state space models (GSSM).

There are two important types of these models "parameter driven" and "observation driven." Both of which are frequently used in time series analysis. In a parameter-driven model, the state equation is constructed by using the previous state. In an observationdriven model, the state equation is constructed by using the past observation.

Parameter-Driven models

The observation and state equations (2.2.3) and (2.2.4) are replaced by the conditional probability densities of the observation and the state variables. The distribution of Y t given (X t , X t-1 , Y t-1 ) is independent of (X t-1 , Y t-1 ), then the observation equation is:

p(Y t | X t ) := p(Y t | X t , X t-1 , Y t-1 ), t = 1, 2, • • • . (2.2.5)
where

Y t = (Y 1 , Y 2 , • • • , Y t ) and X t = (X 1 , X 2 , • • • , X t ). The distribution of X t given (X t-1 , X t-2 , Y t-1 ) is independent of (X t-2 , Y t-1
), then the state equation yields :

p(X t | X t-1 ) := p(X t | X t-1 , X t-2 , Y t-1 ), t = 1, 2, • • • . (2.2.6)
Note : x := y, y =: x or x def = y means x is defined to be another name for y, under certain assumptions taken in context.

Example Assume Y t | µ t ∼ Poisson (µ t ), this means: Pr(Y t = y t | µ t ) = exp(-µ t )µ Yt t µ t ! , t = 1, 2 • • • and log(µ t ) = u t β + X t ,
where β is a parameter vector, (u t ) is the covariate process and {X t } is a stationary gaussian process, e.g. (AR(1) process)

(X t + σ 2 /2) = φ(X t-1 + σ 2 /2) + ε t , ε t ∼ N (0; σ 2 (1 -φ 2 )),
with σ > 0, σ 2 an φ standing for unknown parameters. This example shows that the state variables are a function of previous states.

Observation-Driven models

In an observation-driven model, it is again assumed that Y

t conditional to (X t , X t-1 , Y t-1 ) is independent of (X t-1 , Y t-1
). The model is specified by the conditional probability densities

p(Y t | X t ) = p(Y t | X t , Y t-1 ), t = 1, 2, • • • (2.2.7) p(X t | X t-1 , Y t ) = p(X t | Y t ), (2.2.8) where p(X 1 | Y 0 ) =: p 1 (X
(Y t = y t | µ t ) = exp(-µ t )µ yt t µ t ! , t = 1, 2 • • • and log(µ t ) = u t β + X t ,
where β is a parameter vector, (u t ) is the covariate process and -Filtering for t = T, -Smoothing for t < T, -Prediction for t > T.

{X t , t = 1, 2, • • • } is a function of the past observations Y s , s < t. X t = φ 1 Y t-1 + φ 2 Y t-2 + • • • + φ p Y t-p , with φ i , i = 1
Under the gaussian assumption, the optimal solution to the filtering problem is given by the conditional or posterior mean of X t given Y t

x t|t := E(X t | Y t ).
Since the model is linear and gaussian, the posterior distribution of X t is also gaussian,

X t | Y t ∼ N (x t|t ; V t|t ),
the posterior covariance matrix

V t|t := E[(X t -x t|t )(X t -x t|t ) ].
The famous linear Kalman filter and smoother calculates the posterior means and covariance matrices in an efficient recursive way. The usual derivations of the Kalman filter and smoother take advantage of the fact that the posterior distributions are gaussian.

In Figure (2.1) the Kalman filter algorithm for t = 1, • • • , T, can be performed as follows: start with initial value x 0 . Compute x 1|0 by using x 0 in prediction step, then using x 1|0 to find the x 1|1 via the filter step. Thereafter, compute x 2|1 by using x 1|0 in prediction step and in filter step find x 2|2 by using x 2|1 . Thus, the prediction and filter steps continue until T where each instant t has a prediction value and a filter value. The smoothing step consists of backward recursions for t = T,

• • • , 1. Namely for each instant t = 1, • • • , T -1 computing x t|T depend to the value of instant T.
Therefor, in Figure (2.1), we used the lines to illustrate the three algorithms filtering, prediction and smoothing, another paths can be produced by using these algorithms.

In Figure (2.1), it is important to appoint that we have one filtering step after the prediction step.

Figure 2.1: The paths of recursive computation by the Kalman filter and smoothing algorithm. ⇒: prediction, ⇓: filter, ←: smoothing, →: increasing horizon prediction.

x 3|0 x 2|0 x 1|0 x 4|0 x 5|0 x 3|1 x 2|1 x 1|1 x 4|1 x 5|1 x 3|2 x 2|2 x 1|2 x 4|2 x 5|2 x 3|3 x 2|3 x 1|3 x 4|3 x 5|3 x 3|4 x 2|4 x 1|4 x 4|4 x 5|4 ⇒ ⇓ ⇓ ⇒ ⇓ ⇒ ⇓ ⇒ ⇓ 2.3.1.1 Linear Kalman filter
• The predictive step is x t|t-1 and V t|t-1 is the covariance matrix of predictive step.

• The filter step is x t|t and V t|t is the covariance matrix of filter step.

• The smoother step is x t|T and V t|T is the covariance matrix of smoother step.

The Kalman algorithm for the linear state space models in equations ( 2.2.1, 2.2.2 ) is as follows

Initialization: x 0|0 = x 0 , V 0|0 = R 0 . For t = 1, • • • , T :
Prediction step:

x t|t-1 = F t x t-1|t-1 , V t|t-1 = F t V t-1|t-1 F t + H t R t H t , (2.3.1)
Correction or filtering step:

x t|t = x t|t-1 + K t (Y t -G t x t|t-1 ), V t|t = (I -K t G t )V t|t-1 . (2.3.2)
Kalman gain:

K t = V t|t-1 G t [G t V t|t-1 G t + Σ t ] -1 .
The optimal Kalman gain gives the formula for the updated estimate variance-covariance matrix V t|t . Usage of other gain values leads to a non optimal variance-covariance matrix.

Smoother step: The smoother for X t given Y T where t < T, is as follows :

x t|T := E(X t | Y T ),
and again the posterior is normal,

X t | Y T ∼ N (x t|T ; V t|T ), with V t|T := E[(X t -x t|T )(X t -x t|T ) ].
The ("fixed-interval") smoother consists of backward recursions for t = T, • • • , 1 :

x t|T = x t|t + A t (x t+1|T -x t+1|t ), V t|T = V t|t + A t (V t+1|T -V t+1|t )A t , (2.3.3)
with

A t = V t|t F t+1 V -1
t+1|t . The derivatives of linear Kalman equations are presented in Appendix (A.1).

Extended Kalman Filter and smoother (EKF) recursions

In many practical applications, especially in engineering, most systems are nonlinear. Therefore, some attempts immediately applied the filtering methods to nonlinear systems. Most of this work was done at NASA Ames by McElhoe (1966) 1971) extended an approximate posterior mode in a nonlinear system from conditionally gaussian to exponential family observations. They derived an algorithm that depends on linearizing the non-linear functions by the Taylor expansion about the points recursively. Under the gaussian assumption, the optimal solution to the filtering problem is given by the conditional or posterior mode

x t|t := E(X t | Y t ) of X t given Y t .
Since the model is non-linear and gaussian, the posterior distribution of X t is also approximately gaussian,

X t | Y t ∼ N (x t|t ; V t|t ),
with posterior covariance matrix

V t|t := E[(X t -x t|t )(X t -x t|t ) ].
The derivation of Extended Kalman filter and smoother algorithm is given in Appendix (A.2).

Initialization:

x 0|0 = x 0 , V 0|0 = R 0 .

For t = 1, • • • , T :

Prediction step:

x t|t-1 = F (u t , X t-1|t-1 , λ) V t|t-1 = A t V t-1|t-1 A t + C t R t C t , (2.3.4)
where

A t = ∂F (u t , x, γ) ∂x | x=x t-1|t-1 C t = H t (u t , X t-1|t-1 , δ).
Correction or filtering step:

x t|t = x t|t-1 + K t (Y t -G(u t , x t|t-1 , λ), (2.3.5) V t|t = V t|t-1 -K t B t V t|t-1 , (2.3.6)
where

B t = ∂G(u t , x, λ) ∂x | x=x t-1|t-1 .
Kalman gain:

K t = V t|t-1 B t [B t V t|t-1 B t + Σ t ] -1 .
Smoothing step : EKF's fundamental idea is that linearizing the non-linear functions by the Taylor expansion about the points recursively, several attempts to linearizing the system non-linear implemented as Backward-Smoothing Extended Kalman Filter by Mark (2005). We found challenging to deriving the smoothing step of Extended Kalman Filter with equations ( 2.2.3 and 2.2.4).

General Filtering and smoothing

Nonlinear non-gaussian state space model stratifies the Bayesian perspective by using the conditional densities probability. The prediction, filtering, and smoothing distribution can be approached using the relations p(

X t | X t-1 , Y t-1 ) = p(X t | X t-1 ) and p(Y t | X t , Y t-1 ) = p(Y t | X t ), where Y t = (Y 1 , • • • , Y t ).
p(X t | Y t-1 ) = ∞ -∞ p(X t , X t-1 | Y t-1 )dX t-1 = ∞ -∞ p(X t | X t-1 , Y t-1 )p(X t-1 | Y t-1 )dX t-1 = ∞ -∞ p(X t | X t-1 )p(X t-1 | Y t-1 )dX t-1 . (2.3.7)
The formula (2.3.7) is an extension of the one-step-ahead prediction of the Kalman filter. Here, p(X t | X t-1 ) is the density function of the state X t when the past state X t-1 is given, which is specified by the state equation (2.2.4). Consequently, if the filter p(X t-1 | Y t-1 ) of X t-1 is given, the one-step-ahead predictor p(X t | Y t-1 ) can be estimated.

Filtering distribution

p(X t | Y t ) = p(X t | Y t , Y t-1 ) = p(Y t | X t , Y t-1 )p(X t | Y t-1 ) p(Y t | Y t-1 ) = p(Y t | X t )P (X t | Y t-1 ) p(Y t | Y t-1 ) (2.3.8)
where

p(Y t | Y t-1 ) = ∞ -∞ p(Y t | X t )p(X t | Y t-1 )dX t .
In fact, the filter formula (2.

(X t | X t+1 , Y T ) = p(X t | X t+1 , Y t ),
with the state-space models of (2.2.3) and (2.2.4), then holds as follows

p(X t , X t+1 | Y T ) = p(X t+1 | Y T )p(X t | X t+1 , Y T ) = p(X t+1 | Y T )p(X t | X t+1 , Y t ) = p(X t+1 | Y T ) p(X t | Y t )p(X t+1 | X t , Y t ) p(X t+1 | Y t ) = p(X t+1 | Y T ) p(X t | Y t )p(X t+1 | X t ) p(X t+1 | Y t ) . (2.3.9)
Integrating both sides of (2.3.9), yields the following sequential formula for the smoothing problem:

The smoothing distribution

p(X t | Y T ) = ∞ -∞ p(X t , X t+1 | Y T )dX t+1 = p(X t | Y t ) ∞ -∞ p(X t+1 | Y T )p(X t+1 | X t ) p(X t+1 | Y t ) dX t+1 . ( 2
= T -1, • • • , 1, the smoothing distribution p(X T -1 | Y T ), • • • , p(X 1 | Y T ) is obtained successively.

The sequential Monte Carlo filter

The sequential Monte Carlo filter and smoother is also known as a particle filter or bootstrap filter. Gordon et al. (1993) used the Sequential Monte Carlo in advanced Signal processing and Bayesian inference. Their work was considered as the first application of a Monte Carlo resampling algorithm in Bayesian statistical inference. The authors named their algorithm "the bootstrap filter", and proved that compared to other filtering methods, their algorithm does not require any assumption about the state-space or the noise of the system. Several articles were published on a related "Monte Carlo filter", or " particle filters" by [START_REF] Kitagawa | Monte Carlo filter and smoother for non-Gaussian nonlinear state space models[END_REF], Del Moral (1996 ), and Carvalho, Del Moral, Monin and Salut (1997) .

The famous methods of Monte Carlo filter and smoother are:

1. The Particle Filter (PF).

2. The Auxiliary Particle Filters (APF).

3. The Auxiliary Iterated Extended Kalman Particle Filter (AIEKPF).

Particle Filter (PF)

Kitagawa (1996), and Kitagawa and Gersch (1996) proposed an algorithm that can be applied to general nonlinear state space model. Their algorithm is based on the approximation of successive prediction and filtering density functions by many of their realizations. The difference between this algorithm and other Monte Carlo-Gibbs sampling methods as [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] and [START_REF] Tierney | Markov chains for exploring posterior distributions[END_REF] is that the Monte Carlo method is used for the entire filtering and smoothing procedures, whereas the others are used only for numerical integration. The advantage of this algorithm is that it can be applied to an extensive class of nonlinear non-Gaussian of higher dimensional state space models. They succeeded in estimating the posterior density of the state variables given the observations by using their algorithm.

The sequential Monte Carlo method can characterize the true distribution by using many particles. Each particle is considered as a realization drawn from the true distribution.

The predictive, the filter and the smoothing distributions are approximated by m particles. The number of particles m is usually set between 1000 and 100, 000. The particles number is chosen based on the complexity of the distribution and the necessary accuracy. A true cumulative distribution function can be approximated by an empirical distribution function defined using m particles. The particles {p

(1) t , • • • , p (m)
t } where m is the particles number, that follow the predictive distribution, are generated from the particles {f

(1) t , • • • , f (m) t
} used for the approximation of the filter distribution of the previous state. Then, the realizations {f

(1) t , . . . , f (m) t
} of the filter can be generated by re-sampling the realizations {p

(1) t , • • • , p (m)
t } of the predictive distribution. The table (2.1) shows the distributions and their approximations in the sequential Monte Carlo filter and smoother algorithm.

t , • • • , p (m) t } Filter distribution p(X t | Y t ) {f (1) t , • • • , f (m) t } Smoothing distribution p(X t | Y T ) {s (1) t|T , • • • , s (m)
t|T } Distribution of state noise p(ε t ) {ε

(1)

t , • • • , ε (m)
t }

One-step-ahead prediction

The particles {f

(1) t , • • • , f (m) t
} can be considered as realizations of the filter step, which can be generated from the filter distribution p(X t-1 | Y t-1 ) of the previous step X t-1 by m particles. The particles ε

(1) t , • • • , ε (m) t
can be considered as independent realizations of the state noise ε t . For j = 1, • • • , m, we can write

f (j) t-1 ∼ p(X t-1 | Y t-1 ), ε (j) t ∼ q(ε t ). (2.3.11) Therefore, p (j)
t is the jth particles defined by using the state equation (2.2.4)

p (j) t = F (u (j) t , f t -1 (j) , γ) + H(u (j) t , f (j)
t-1 , δ)ε

(j) t , (2.3.12)
where {p

(1)

t , • • • , p (m)
t } can be considered as independent realization of the one step ahead predictor distribution of the state X t .

Filtering

In the filtering step, the Bayes factor (or likelihood) α (j) t of the particle p (j) t can be computed with respect to the observation Y t (t). For j t can be considered as a weighting factor representing the importance of the particle p

= 1, • • • , m, α (j) t = r t (S t (Y t , p (j) t )) × ∂S t ∂Y t , (2.3
(j) t . Then, m particles f (1) t , • • • , f (m) t are obtained by re-sampling p (1) t , • • • , p (m) t
with probabilities proportional to the "likelihoods" α

(1) t , • • • , α (m) t . That is, a new particle f (j) t , j = 1, • • • , m is obtained according to f (j) t =      p (1) t probability α (1) t /(α (1) t + • • • + α (m) t ) . . . p (m) t probability α (m) t /(α (1) t + • • • + α (m) t ). (2.3.14) Then, f (1) t , • • • , f (m) t
can be considered as the realizations generated from the filter distribution p(X t | Y t ). The derivation of particle filter algorithm is presented in Appendix (A.3 ).

Algorithm for the Monte Carlo filter [START_REF] Kitagawa | Monte Carlo filter and smoother for non-Gaussian nonlinear state space models[END_REF] described the Monte Carlo filter algorithm as follows:

1. A k-dimensional random number Generate from f (j) 0 ∼ p 0 (X) for j = 1, • • • , m. k is the dimensional for state vector X t 2. For t = 1, • • • , T, repeat the following steps (a) Generate -dimensional random numbers ε (j) t ∼ q(ε t ), j = 1, • • • , m. (b) Find the new particles: for j = 1, • • • , m. p (j) t = F (u (j) t , f (j) t-1 , γ) + H(u (j) t , f (j) t-1 , δ)ε (j) t , (c) Calculate the Bayes factor : for j = 1, • • • , m. α (j) t = r t (S t (Y t , p (j) t )) × ∂S t ∂Y t , (d) Generate {f (1)
t , • • • , f (m) t
} by applying re-sampling method (sampling with replacement) m times from {p

(1) t , • • • , p (m)
t } with probabilities proportional to {α

(1) t , • • • , α (m) t }.

Re-sampling method

The re-sampling method is essentially based on random sampling as described below For j = 1, • • • , m, repeat the following step (a) -(c).

(a) Generate uniform random number z

(j) t ∈ U [0; 1]. (b) Search for i that achieves C -1 i-1 l=1 α (l) t < z (j) t ≤ C -1 i l=1 α (l) t where C = m l=1 α (l) t .
(c) Get a particle which approximates the filter by setting f

(j) t = p (j)
t . The objective of re-sampling is to re-express the distribution function specified by the particles {p

(1) t , • • • , p (m) t } with weights {α (1) t , • • • , α (m)
t } through representing the empirical distribution function determined by re-sampled particles with equal weights. As consequence, it is not essential to perform exact random sampling.

Monte Carlo smoothing method

The Monte Carlo smoothing method is created by extending the Monte Carlo filtering method through preserving past particles. In the following, the vector of the particles (s

(j) 1|t , • • • , s (j) t|t ) indicates the j-th realization from the t-dimensional joint distribution function p(X 1 , • • • , X t | Y t ).
It should be only modifying step 2(d) of the Monte Carlo filtering algorithm to achieve the objective of smoothing as follows.

(d-S) For j = 1, • • • , m by re-sampling the t-dimensional vector (s

(j) 1|t-1 , • • • , s (j) t|t-1 , p (j) t ) , generate (s (j) 1|t , • • • , s (j) t|t , s (j) t|t ) .
In this modification, by re-sampling {(s

(j) 1|t-1 , • • • , s (j) t-1|t-1 , p (j) t ) , j = 1, • • • , m}
with the same weights as in step 2 -d, fixed interval smoothing for non-linear non-Gaussian state space model can be achieved [START_REF] Kitagawa | Monte Carlo filter and smoother for non-Gaussian nonlinear state space models[END_REF]). In practical computation, since the re-sampling methods repeated a finite number of particles (m particles), the number of different particles progressively decreases and then the weights become concentrated on a small number of particles. As result, the shape of distribution finally collapses. Accordingly, in the smoothing algorithm, it is necessarily the step (d -S) performs as follows:

(d-L) For j = 1, • • • , m, generate (s (j) t-L|t , • • • , s (j) t-1|t , s (j) t|t ).
Here, L is a fixed integer, usually less than 30, and assume f

(j) t = s (j)
t|t by re-sampling (s

(j) t-L|t-1 , • • • , s (j) t-1|t-1 , p (j) t ).
Then, the modification to algorithm turns out to conform to the L-lag fixed-lag smoothing algorithm. If L is fixed and large, the fixed interval smoother p(X t | Y T ) is approximated by the fixed-lag smoother p(X t | Y t+L ). Furthermore, the distribution is determined by X 1 t|t+L , • • • , X m t|t+L not distant from p(X t | Y t+L ). Subsequently, L should be taken not very large, i.e., L = 20 or 30.

The Auxiliary Particle Filters (APF)

Pitt and Shephard (1999) proposed the Auxiliary Particle Filters (APF). They extended standard particle filtering methods by adding an auxiliary variable which permits the particle filter to be adapted more efficiently. The work of the auxiliary variable ς, is only to improve a simulation performance. The essential idea of the APF is to insert an auxiliary variable ς. It plays an important role of index of the mixture component. The augmented joint distribution p(X t , ς | Y t ) with this additionally auxiliary variable is updated as :

p(X t , ς = m | Y t ) ∝ p(Y t | X t )p(X t , ς = m | Y t ) = p(Y t | X t )p(X t | ς = m, Y t ) ×p(m | Y t ). = p(Y t | X t )p(X t | X m t-1 )w m t-1 , (2.3.15)
where y ∝ x means that y = kx for some constant k. Thus, the APF generates a sample {X m t , ς m } N m=1 from the joint distribution p(X t , ς | Y t ), where ς m (ς m ≡ {ς = m}) refers to the index of particle m s parent. The importance density is used in APF as follows

q(X t , ς | Y t ) ∝ p(Y t | µ m t )p(X t | X m t-1 )w m t-1 , (2.3.16)
where µ m t is some description of X t given X m t-1 , such as the mean, mode or another statistics. If it is the mean, then µ m t = E(X t | X m t-1 ) or a sample µ m t ∼ p(X t | X m t-1 ). Then, by forgetting the auxiliary variable, the sample {X m t } N m=1 is simply obtained. The importance density is written as

q(X t , ς | Y t ) = q(ς | Y t )q(X t | ς, Y t ).
(2.3.17)

Then, defining

q(X t | ς m , Y t ) := p(X t | X m t-
q(ς m | Y t ) ∝ p(Y t | µ m t )w m t-1 . (2.3.19)
Recursively, the importance weights are updated as

w m t ∝ p[X m t , ς m | Y t ] q[X m t , ς m | Y t ] . ( 2
w m t ∝ p(Y t | X m t )p(X t | X m t-1 )w m t-1 p(Y t | µ ς m t )p(X t | X m t-1 )w m t-1 = p[Y t | X m t ] p[Y t | µ ς m t ]
.

(2.3.21)

Importance density: Importance sampling is favourite among the variance reduction techniques that can be used with complex distributions when Monte Carlo sampling is used. The necessary procedure in importance sampling is to choose a density function which "encourages" the important values. If the importance density is applied primary in the simulation, the result is a biased estimator. However, the new importance sampling estimator is unbiased due to the fact that the outputs of a simulation are weighted to correct for the use of the biased distribution. On other words, the importance weight is given by w t = p(x) q(x) , which is called the likelihood ratio. q(x), p(x) are known as the importance density and the nominal density respectively. The good choice of importance density leads to minimize the variance of importance weights [START_REF] Kroese | Monte carlo methods[END_REF]).

Auxiliary Particle Filter (APF) Algorithm

The APF algorithm can be performed as follows:

1. Initialization (t = 0) : for m = 1, • • • , N , generate the states (particles) X m 0 from the prior p(X 0 ), set

x m 0 = E[X m 0 ] = F (u m 0 , x m 0 , δ m ) and V m 0 = E[(X m 0 (0) -x m 0 )(X m 0 -x m 0 ) ] = Var(X m 0 ) = H 2 (u m 0 , x m 0 , δ m )R 0 .
For t = 1, • • • , T, repeat the following steps:

2. For m = 1, • • • , N, generate µ m t ∼ p(X t | X t-1 ). 3. For m = 1, • • • N, compute w m t = q(m | Y t ) ∝ p(Y t | µ m t )w m t-1 .
4. Resample to obtain the index ς m of particle m's parent.

For

m = 1, • • • , N, generate X m t ∼ p(X m t | X m t-1 , ς m ).
6. Update the second-stage weights

w m t = p[Y t | X m t ] p[Y t | µ ς m t ]
.

-Normalize the weights to avoid the Degeneracy problem is as follows

w m t = w m t N m=1 w m t . 7. Output: a set of weighted particles (samples) [{X m t , w m t } N m=1 ]. 2
p(X t | Y t-1 ) ≈ N m=1 w (m) t|t-1 δ(X t -X (i) t ),
where δ(.) indicates the Dirac delta function. The posterior mean are approximated by

X t|t-1 ≈ N m=1 w (m) t|t-1 X (i) t .
Then, the posterior distribution can be approximated as

p(X t | Y t ) ≈ N m=1 w m t δ(X t -x t ) (2.3.22)
where

x t = (x m 1 , • • • , x m t
) is a set of particles, the weights are normalized such that N m=1 w m t = 1. The weight w m t is chosen using the precept of importance sampling, which is given by

w m t (t) ∝ p[X m t | Y t ] q[X m t | Y t ] , (2.3.23)
where q[X m t | Y t ] is known as the importance density. As we mentioned earlier in APF algorithm the auxiliary variable, ς, played a role of index of the mixture component.

The augmented joint distribution

p(X t , ς | Y t ) is updated as: p(X t , ς = m | Y t ) ∝ p(Y t | X t )p(X t , ς = m | Y t ) = p(Y t | X t )p(X t | ς = m, Y t ) ×p(m | Y t ). = p(Y t | X t )p(X t | X m t-1 )w m t-1 . (2.3.24)
Thus, the importance density is

q(X t , ς | Y t ) ∝ p(Y t | µ m t )p(X t | X m t-1 )w m t-1 , (2.3.25)
where µ m t is mean, mode or another statistic. Then, by deleting the auxiliary variable, the sample {X m t } N m=1 is obtained. The importance density yields

q(X t , ς | Y t ) = q(ς | Y t )q(X t | ς, Y t ). (2.3.26) Then, defining q(X t | ς m , Y t ) = p(X t | X m t-1 ), (2.3
q(ς m | Y t ) ∝ p(Y t | µ m t )w m t-1 . (2.3.28)
Recursively, the importance weights are updated as

w m t ∝ p[X m t , ς m | Y t ] q[X m t , ς m | Y t ] . ( 2
w m t ∝ p(Y t | X m t )p(X t | X m t-1 )w m t-1 p(Y t | µ ς m t )p(X t | X m t-1 )w m t-1 = p[Y t | X m t ] p[Y t | µ ς m t ] . (2.3.30)
Occasionally, the likelihood function is also narrow compared to the prior distribution, or it happens to lie in the tails of the prior distribution, then the posterior distribution is poorly approximated by the importance density. To devise the optimal importance density, the IEKF is used to update particles by blending the most current observation with the optimal gaussian approximation of the state. In other words, the IEKF computes the following recursive approximation to the real posterior filtering density

p(X t | Y t ) = N X tj ; V tj (2.3.31) where t = 0, • • • , T, j = 1, • • • , c
and c is the number of the Iterated Extended Kalman Filter (IEKF), X tj is the iterative value of the X t on the jth iteration, and V tj is the covariance matrix of X tj . The derivation of equations to iterated extended Kalman filter for the non-linear non-gaussian as equations (2.2.3, 2.2.4) in Appendix A.2. According to the IEKF update equation, X tj and V tj can be updated as follows :

X t|t,j = X t|t-1,j + K tj [Y t -G(u t , X t , λ)] (2.3.32) K tj = V t|t-1,j B tj [B tj V t|t-1,j B tj + Σ t ] -1 V t|t,j = (I -K tj B tj )V t|t-1,j (2.3.33)
where B tj is the Jacobian of G t , the function appearing in the observation equation (2.2.3), and K tj is a Kalman gain matrix. Within the APF algorithm, the importance density is generated for each particle by the IEKF algorithm as follows

q(X m t , ς | Y t ) = N X ς m tj ; V ς m tj (2.3.34)
where m = 1, 2, • • • N, and N is the number of the particles.

In short, the AIEKPF uses the IEKF to update the equations with the new observations to compute the mean and covariance of the importance density for each particle at time t -1. Next, the mth particle is sampled from the distribution.

Auxiliary Iterated Extended Kalman Filter (AIEKPF) Algorithm

The proposed AIEKPF algorithm can be performed by the following steps:

1. Initialization (t = 0) : For m = 1, • • • , N , generate the states (particles) X m 0 from the prior p(X 0 ), and set

x m 0 = E(u m 0 , X m 0 ) = F (x m 0 , δ m )
and

V m 0 = E[(X m 0 (0) -x m 0 )(X m 0 -x m 0 ) ] = Var(X m 0 ) = H 2 (u m 0 , x0 m , δ m )R 0 For t = 1, • • • , T, repeat the following steps : 2. For m = 1, • • • , N , generate µ m t ∼ p(X t | X t-1 ).
3. For m = 1, • • • , N , using the IEKF algorithm to update the particles -Calculate the Jacobians A m t , C m t of the process model

A m t = ∂F (u t , x, γ) ∂x | x=x m t-1|t-1 C m t = H(u t , x m t-1|t-1 , δ)
-Using the IEKF to predict the particle :

x m t|t-1 ≈ F (u t , x m t-1|t-1 , γ) V m t|t-1 = A m t V m t-1|t-1 A m t + C m t R t C m t -For j = 1, • • • , c (c is the iteration number of the IEKF) a -Compute the Jacobian B m tj of G t function in observation equation B m tj = ∂G(u t , x, λ) ∂x | x=x m t|t-1,j b -

For

m = 1, • • • N, Compute w m t = q(m | Y t ) ∝ p(Y t | X t ) × w m t-1 , (2.3.35)
5. Resample to find the index ς m of particle m's parent.

Find the importance sampling : for

m = 1, • • • , N , -Draw samples X m t ∼ q(X t , ς m | Y t ) = N ( X ς m tj ; V ς m tj ), j = c, (2.3.36)
we recall c is the iteration number of the IEKF.

-Compute the importance weights of particles as follows

w m t = p(Y t | X m t ) p(Y t | µ ς m t ) (2.3.37)
-Normalize the weights to avoid the Degeneracy problem is as follows

w m t = w m t N m=1 w m t . 7. Output: a set of weighted particles (samples) [{X m t , w m t } N m=1 ].
Remark: practically, the iteration to perform the equations (2.3.36, 2.3.37) leads to Degeneracy problem. On other words, all the particles will have negligible weights but a few of the particles will have a significant weight. The term effective sample size is used to measure the degeneracy problem as follows:

N ef f = 1 N m=1 (w m t ) 2
where a smaller N ef f means a larger variance for the weights. One of the techniques to solve this problem is

Resampling with replacement N particles { X m t , m = 1, • • • , N } from the set {X m t , m = 1, • • • , N }.
The importance weights can be normalized as follows :

w m t = w m t N m=1 w m t . 2.3.
√ N ( ϕ t -φt ) ⇒ N (0; σ 2 AP F (ϕ t )), (2.3.38)
where ϕ t is the function of trajectories X t and φt is the expectation of ϕ t with respect to the filtering distribution.

φt = ϕ t p(X t | Y t )dX t and ϕ t = N m=1 W m t ϕ m t , (2.3.39)
where

W m t = w m t [ N m=1 w m t ] -1
, and

(2.3.40)

w m t = p[Y t | X m t ] p[Y t | µ ς m t t ] . At time t = 1, obtain σ 2 AP F (ϕ 1 ) = (p(X 1 | Y 1 ) 2 q 1 (X 1 ) (ϕ 1 -φ1 ) 2 dX 1 , (2.3.41)
whereas for t > 1,

σ 2 AP F (ϕ t ) = [p(X 1 | Y t )] 2 q 1 (X 1 ) ϕ 2 1,t dX 1 + t-1 r=2 [p(X 1:r | Y t )] 2 p(X r-1 | Y r )q r (X r | X r-1 ) × ϕ 2 r,t dX r + [p(X t | Y t ] 2 p(X t-1 | Y t )q t (X t | X t-1 ) ×(ϕ t -φt ) 2 dX t . (2.3.42) where ϕ r,t = ϕ t p(X r+1:t | Y r+1:t )dX r+1:t -φt

Maximum a posteriori estimation (MAP)

The Bayesian approach certainly leads to compute the posterior distribution by using the particle filtering methods. If the posterior distribution is gaussian, the point that maximizes the posterior distribution is posterior mean. In contrast, the mode maximize the non-gaussian posterior distribution which is called maximum a posteriori probability (MAP).

X M AP t = arg max Xt p(X t | Y t-1 ), t = 1, • • • , T. (2.3.43)
The MAP can be used to compute an estimation of state variable X t . We present the MAP via three methods:

1. Viterbi algorithm.

2. Particle filter algorithm.

3. Auxiliary Iterated Extended Kalman Filter (AIEKPF) algorithm.

Particle filters MAP via the Viterbi algorithm

Godsill et al.(2000) developed a technique which was implemented by using maximum a posteriori (MAP) sequence estimation in non-linear non-gaussian dynamic models. They used Monte Carlo filtering and Viterbi algorithm, and focused on estimating the MAP sequence as follows

X M AP t = (X M AP 1 , • • • , X M AP t ) X M AP t ∼ = arg max Xt (p(X t | Y t )) = arg max Xt p(Y t | X t ).p(X t | Y t-1 ) p(Y t | Y t-1 ) (2.3.44)
The dependency upon Y t is indicated by t. Note that p(Y t | Y t-1 ) does not depend on the value of X t . Consequently, it can be equivalently written

X M AP t ∼ = arg max Xt [p(Y t | X t ).p(X t | Y t-1 )] (2.3.45) practically, computing p(X t | Y t-1
) can only be executed in closed form for linear gaussian models by using the Kalman filter-smoother and for finite state space hidden Markov models. In contrast, the numerical techniques must be employed to approximate the posterior distribution, such as the extended Kalman filter, gaussian sum methods and general numerical integration procedures Monte-Carlo filtering (particle filtering) as the methods proposed by Kitagawa (1987), [START_REF] Kitagawa | Monte Carlo filter and smoother for non-Gaussian nonlinear state space models[END_REF]

p(X t | Y t-1 ) = p(X t , X t-1 | Y t-1 )dX t-1 = p(X t | X t-1 , Y t-1 )p(X t-1 | Y t-1 )dX t-1 = p(X t | X t-1 )p(X t-1 | Y t-1 )dX t-1 , (2.3.46)
the approximation can be made at time t as follows

p(X t | Y t ) ≈ N m=1 w (m) t δ(X t -x m t ), (2.3.47)
where δ indicates the Dirac delta function and w (m) t

is the weight joined to particle x

(m) t , w (m) t ≥ 0 and N m=1 w (m) t = 1. The approximation of p(X t | Y t-1 ) is as follows p(X t | Y t-1 ) ≈ N j p(X t | X j t-1 )w j t-1 . (2.3.48)
The Viterbi algorithm calculate the MAP as

X M AP t ∼ = arg max Xt t k=1 [log p(Y k | X k ) + log p(X k | X k-1 )]. (2.3.49)
The Viterbi algorithm is a famous technique to estimating discrete state-space hidden Markov models. [START_REF] Doucet | On sequential Monte Carlo sampling methods for Bayesian filtering[END_REF] developed this algorithm to estimating continuous state-space Markov model via a discrete approximation of the state space as particles.

The Viterbi algorithm is summarized as follows:

Viterbi algorithm

Initialization : For 1 ≤ m ≤ N δ 1 (m) = log p(X (m) 1 | X (m) 0 ) + log p(Y 1 | X (m) 1 ) (2.3.50) Recursion: For 2 ≤ k ≤ t and 1 ≤ j ≤ N δ k (j) = log p(Y k | X (j) k ) + max m δ k-1 (m) + log p(X (j) k | X (m) k-1 ) ψ k (j) = arg max m δ k-1 (m) + log p(X (j) k | X (m) k-1 ) (2.3.51)
Termination:

m t = arg max m δ t (m) X M AP t = X mt t (2.3.52) Backtracking: For k = t -1, t -2, • • • , 1 m k = ψ k+1 (m k+1 ) X M AP k = X (m k ) k (2.3.53) 2.3.5.2
The maximum a posteriori (MAP) with particle filtering and smoothing algorithm [START_REF] Saha | On the Monte Carlo marginal MAP estimator for general state space models[END_REF] proposed a method to estimating MAP to a general state space model, in which they utilised gradient based optimization method. In this method, the MAP locate can be approximated by computing the posterior (filtering) density at the predicted particle {x (m) t } N m=1 , and choosing the particle with the highest density. One of the advantages of this method is that the posterior density is approximate to any specific point, in addition to be approximate to particles forming the clouds.

They compared the particle filters using the Viterbi algorithm and their method (pf-MAP) via numerical simulations. They observed that their method was better, as an estimate of current state space. They concluded that their results were not certainly similar to the Viterbi algorithm results. Furthermore, the Viterbi algorithm goal is maximization of p(x t | y t ), whereas the pf-MAP goal is maximization of p(x t | y t ). They noted that the two estimates would be principally similar in the linear-gaussian state space models. In contrast, the difference would certainly appear for the above two estimators in nonlinear and/or non-gaussian state space models. They investigated the behaviours of the two algorithms through numerical simulations. First, they studied the linear gaussian state space model and compared these two estimators. The criteria is used a root mean square error estimate (RMSE) to a true MAP. The RMSE's for the two methods behaved similarly. Second, they took the nonlinear-gaussian state space model, three different initial values for estimating the RMSE are used. They observed through RMSE estimate that the pf-MAP implements better as an estimate of current state whereas it is also computationally much easy. The filtering and smoothing steps to find the MAP as follows:

Filtering step

The MAP estimate of the filtering density at time t is defined as

X M AP t = arg max Xt p(X t | Y t ), (2.3.54)
where

Y t = (Y 1 , • • • , Y t ) and X t = (X 1 , • • • , X t )
. By using the Baye's rule, the posterior (filtering) density can be written as

p(X t | Y t ) = p(X t | Y t , Y t-1 ) = p(Y t | X t , Y t-1 )p(X t | Y t-1 ) p(Y t | Y t-1 ) = p(Y t | X t )p(X t | Y t-1 ) p(Y t | Y t-1 ) (2.3.55)
Therefore, the MAP estimate of the posterior distribution is

X M AP t = arg max Xt p(Y t | X t )p(X t | Y t-1 ) p(Y t | Y t-1 ) , (2.3.56) the denominator p(Y t | Y t-1 ) is independent of X t , the MAP estimate yields X M AP t arg max Xt [p(Y t | X t )p(X t | Y t-1 )] . (2.3.57)
Since the conditional likelihood p(Y t | X t ) is known for each X t . The main problem for computing MAP is the computation of the predictive density p(X t | Y t-1 ) (the second terms in equation (2.3.57) ), whereas it not available in closed form. Nonetheless, using the formula of prediction density in Bayesian approach yields

p(X t | Y t-1 ) = p(X t , X t-1 | Y t-1 )dX t-1 = p(X t | X t-1 , Y t-1 )p(X t-1 | Y t-1 )dX t-1 = p(X t | X t-1 )p(X t-1 | Y t-1 )dX t-1 , (2.3.58)
Then, using the equation (2.3.58) one obtains

p(X t | Y t-1 ) = p(X t | X t-1 )p(X t-1 | Y t-1 )dX t-1 .
(2.3.59)

Using particle filtering technique, this posterior distribution can be approximated by a set of N weighted particles as

p(X t | Y t ) ≈ N m=1 w m t δ(X t -x m t ) (2.3.60)
where

x m t = (x m 1 , • • • , x m t ) . Now, use equation (2.3.60) to approximate p(X t | Y t-1 ) as p(X t | Y t-1 ) ≈ N j p(X t | X j t-1 )w (j) t-1 . (2.3.61)
Substituting (2.3.61) into (2.3.57) the MAP estimation is obtained by computing the global maxima of the posterior distribution by particles as

X M AP t = arg max X (m) t p(Y t | X (m) t ) × N j p(X (m) t | X (j) t-1 )w (j) t-1 . (2.3.62)
It should be noted that for each time step, the memory necessity of this MAP estimator is O(N ) and the computational complexity is O(N 2 ).

X M AP t|T = arg max Xt p(X t | Y T ) = arg max X (m) t p(X (m) t | Y t ) w (m) t|T w (m) t , (2.3.63)
where for the filtering density p(X t | Y t ) at the particle cloud {X (m) t

} N m=1 the evaluation during the forward filtering step can be defined as

p(X t | Y t ) ≈ p(Y t | X (m) t ) N j p(X (m) t | X (j) t-1 )w (j) t-1 p(Y t | Y t-1 ) , (2.3.64) since p(Y t | Y t-1 ) in equation (2.3.64) is independent of X (m) t . Consequently, obtaining x M AP t|T requires to replace p(X (m) t | Y t ) by the filtered density q(X (m) t | Y t ) = p(Y t | X (m) t ) N j p(X (m) t | X (j) t-1 )w (j) t-1 .
(2.3.65)

The derivation of the Forward-Backward smoothing is given in Appendix A.4.

The MAP via Particle filter Algorithm :

• Given observations Y t = (Y 1 , • • • , Y t ) , • For m = 1, • • • , N
where N is the number of particles Forward filtering step:

• Assume p(X 0 ), draw X (m) 0 from p(X 0 ), set w (m) 0 = 1 N .
• Apply particle filter to generate and store {X

(m) t , w (m)
t } for t = 0, • • • , T.

• Evaluate (un-normalized) filtering distribution for t = 1,

• • • , T, at cloud points m q(X (m) t | Y t ) = p(Y t | X (m) t ) N j p(X (m) t | X (j) t-1 )w (j) t-1 . starting with q(X (m) 0 ) = p(X (m)
0 ) and store Backward smoothing step:

• Set w (m) T |T = w (m) T • For t = T -1, • • • , 1 the smoother importance weights is computed as w (m) t|T = w (m) t N j=1 w (m) t+1|T p(X (j) t+1 | X (m) t ) N k=1 p(X (j) t+1 | X (k) t )w (k) t .
• Compute the approximate smoother MAP as

X M AP t|T = arg max X (m) t q(X (m) t | Y t ) w (m) t|T w (m) t .

Auxiliary Iterated Extended Kalman Particle Filter-MAP(AIEKPF-MAP) Algorithm

We can extend the MAP estimation algorithm by using the weights calculated by Auxiliary Iterated Extended Kalman Particle Filter (AIEKPF) algorithm :

w m t = p[Y t | X m t ] p[Y t | µ ς m i t ]
.

The AIEKPF algorithm can be performed in the following steps:

1. Initialization (t = 0) : For m = 1, • • • , N , generate the states (particles) X m 0 from the prior p(X 0 ), and set

X m 0 = E(X m 0 ) = F (u m 0 , X m 0 , δ m )
and

V m 0 = E[(X m 0 -X m 0 )(X m 0 -X m 0 ) ] = Var(X m 0 ) = H 2 (u m 0 , X m 0 , δ m )R 0 For t = 1, • • • , T, repeat the following steps: 2. For m = 1, • • • , N , generate µ m t ∼ p(X t | X t-1 ).
3. For m = 1, • • • , N , using the IEKF algorithm to update the particles -Compute the Jacobians A m t , C m t of the process model

A m t = ∂F (u t , x, γ) ∂x | x=X m t-1|t-1 C m t = H(u t , X m t-1|t-1
, δ) -Predict the particle with the IEKF:

x m t|t-1 ≈ F (u t , x m t-1|t-1 , γ) V m t|t-1 = A m t V m t-1|t-1 A m t + C m t R t C m t -For j = 1, • • • c (c

For

m = 1, • • • N, calculate w m t = q(m | Y t ) ∝ p(Y t | X t ) × w m t-1 , (2.3.66)
5. Perform the Re-sampling methods to obtain the index ς m of particle m's parent.

6. Find the importance sampling :

for m = 1, • • • , N , -Draw samples X m t ∼ q(X t , ς m | Y t ) = N ( X ς m tj ; V ς m tj ), j = c (2.3.67)
-Compute the importance weights of particles by using

w m t = p(Y t | X m t ) p(Y t | µ ς m t ) (2.3.68)
-Normalize the weights

w m t = w m t N m=1 w m t .
Forward Filtering step :

X M AP t = arg max X (m) t p(Y t | X (m) t ) j p(X (j) t | X (j) t-1 )w (j) t-1
Backward smoothing step :

• Set w (m) T |T = w (m) t • For t = T -1, • • • , 1, compute the smoother importance weights as w (m) t|T = w (m) t N j=1 w (j) t+1|T p(X (j) t+1 | X (m) t ) N r=1 p(X (j) t+1 | X (r) t )w (r) t
• Compute the approximate smoother MAP as

X M AP t|T = arg max X (m) t q(X (m) t | Y t ) w (m) t|T w (m) t .

The state space representation and latent variable models

State space models methodology provides a good technique for analyzing non-gaussian time series and longitudinal data, as a model for discrete longitudinal observations. In this section, we present two applications of non-gaussian state space model. First, the exponential family state space models (generalized linear models) was studied by Fahrmeir and Wagenpfiel (1997)

) : i = 1, • • • , n, k = 1, • • • , q, t = 1, • • • , T ),
where n is the number of individual and q is the number of item and t the instant.

The exponential family state space models

Consider the time series observations {Y t } and the state variables {X t } with dimensions and k, respectively. The exponential family observation equation is

Y t | X t ∼ p(Y t | X t ) = exp θ t Y t -b t (θ t ) φ + c t (Y t , φ), (2.4.1)
where θ t the natural parameter is a function of η t = Z t X t , Z t is a × k matrix. φ is a dispersion parameter, and c t (.) and b t (.) are known functions. By the properties of exponential families, mean and variance functions can be written as

E(Y t | X t ) = µ t = ∂b t (θ t ) ∂θ t (2.4.2) Var(Y t | X t ) = Σ t = ∂ 2 b t (θ t ) ∂θ t ∂θ t . (2.4.3)
As in generalized linear models, the mean µ t is related to the linear predictor η

t = Z t X t by µ t = h(Z t X t ), ( 2
t = F t X t-1 + ε t , (2.4.5)
where F t is the transition matrix and {ε t } is a Gaussian noise with ε t ∼ N (0; R t ). The initial state X 0 ∼ N (x 0 ; R 0 ). Fahrmeir and Wagenpfiel (1997) assumed the following conditional independence assumptions:

(A 1 ) Conditional on X t , the current responses Y t are independent of past states

X t-1 , • • • , X 0 , i.e. p(Y t | X t , X t-1 , • • • , X 0 ) = p(Y t | X t ), t = 1, • • • , T. (A 2 ) The sequence {X t , t = 1, 2, • • • } is Markovian, i.e. p(X t | X t-1 , • • • , X 0 ) = p(X t | X t-1 ).
According to (2.4.5) one obtains

p(X t | X t-1 ) ∼ N (F t X t-1 , R t ).

Penalized likelihood estimation

Fahrmeir and Wagenpfeil (1997) proposed a simple method of inference founded on posterior mode, or equivalently maximum penalized likelihood estimation. Their approach can also be interpreted as a nonparametric method for smoothing time-varying coefficients. They estimated the posterior mode by using Fisher scoring method via the iterative Kalman filtering and smoothing. They used the numerical method to maximum penalized likelihood called "Working Kalman filtering and smoother (WKFS)" with the exponential family distribution. Define

Y t = (Y 0 , Y 1 , • • • , Y t ) , and X t = (X 0 , X 1 , • • • , X t )
, the posterior mode smoother is defined as

a ≡ a 0|T , a 1|T , • • • , a T |T ∈ R m , t = 1, • • • , T.
with m = (T + 1)k, The posterior distribution of X T given Y T is obtained by Bayes' theorem :

p(X T | Y T ) = 1 p(Y T ) T t=1 p(Y t | X t , Y t-1 ) × T t=1 p(X t | X t-1 , Y t-1 ) × p(X 0 ) (2.4.6) p(Y T ) does not depend on X T p(X T | Y T ) ∝ T t=1 p(Y t | X t , Y t-1 ) T t=1 p(X t | X t-1 ).p(X 0 ). (2.4.7)
Taking logarithms and substituting the densities of X t of the state equation (2.4.5) and X 0 , one can obtain the penalized log-likelihood function

P L : R m -→ R, m = (T + 1)k P L(X t ) = T t=1 log p(Y t | X t , Y t-1 ) + T t=1 log p(X t | X t-1 ) + log p(X 0 ), (2.4.8)
which can be written as:

P L(X t ) = T t=1 p(Y t | X t , Y t-1 ) + G 1 + G 2 , (2.4.9)
where

G 1 = - 1 2 (X 0 -a 0 )R -1 0 (X 0 -a 0 ) G 2 = - 1 2 T t=1 (X t -F t X t-1 ) R -1 t (X t -F t X t-1 ).
The densities P (Y t | X t , Y t-1 ) are defined by the exponential family observation model equation (5.17). Thus, the posterior mode smoother is given by

a ≡ (a 0|T , a 1|T , • • • , a T |T ) = arg max X {P L(X)}, (2.4.10)
The Maximization of p(X T | Y T ) is equivalent to the maximization of the penalized log-likelihood (2.4.9). Numerical maximization of the penalized log-likelihood can be achieved by various algorithms. [START_REF] Fahrmeir | Posterior mode estimation by extended Kalman filtering for multivariate dynamic generalized linear models[END_REF] suggested the generalized extended Kalman filter and smoother as an approximative posterior mode estimator in dynamic generalized linear models. Fahrmeir and Kaufmann (1991) developed iterative forward-backwards Gauss-Newton (Fisher-scoring) algorithms. The iterative application of linear Kalman filtering and smoothing to a "working" model give the same results by using Gauss-Newton smoothers.

Gauss-Newton and Fisher-scoring Filtering and smoothing

A maximization of the penalized log-likelihood P L(X t ) with the best performance to approximation can be computed by Gauss-Newton or Fisher-scoring iterations. In other words, this can also be achieved by applying the linear Kalman filtering and smoothing to a "working model" in each Fisher-scoring iteration. Then, the penalized log-likelihood criterion (2.4.8) can be written in a matrix notation as:

P L(X t ) = l(X T ) - 1 2 X t KX T , (2.4.11)
where

l(X t ) = T t=0 l t (X t ), with l t (X t ) = ln p(Y t | X t ), t = 0, • • • , T,
and l 0 = -(X 0 -a 0 ) R -1 0 (X 0 -a 0 )
)/2, the penalty matrix K is symmetric and block-tridiagonal, with blocks that can easily be computed from (2.4.8):

K =         K 00 K 01 • • • • • • 0 K 10 K 11 K 12 • • • . . . . . . K 21 . . . . . . . . . . . . . . . . . . K T -1,T 0 • • • • • • K T,T -1 K T,T         with K t-1,t = K t,t-1 , t = 1, • • • , T K 00 = F 1 R -1 1 F 1 , K tt = R -1 t + F t+1 R -1 t F t+1 , t = 1, • • • , T F T +1 = 0, K t-1,t = -F t R -1 t , t = 1, • • • , T. R t is
Y = (Y 0 , Y 1 , • • • , Y T ) .
Fahrmeir and Wagenpfeil (1997) assumed Y 0 = a 0 . In contrast, the matrices of expectations are defined by adding

µ 0 µ = (µ 0 , µ 1 , • • • , µ T ) .
Fahrmeir and Wagenpfeil (1997) assumed µ 0 = X 0 , where µ t = Z t X t , the blockdiagonal covariance matrix

Σ = diag(R 0 , Σ 1 , • • • , Σ T ), the block-diagonal design matrix Z = diag(I, Z 1 , • • • , Z T ),
with I ∈ R k×k is the unit matrix and the block-diagonal matrix

D = diag(I, D 1 , • • • , D T ),
where for all t = 1, . . . , T, D t = ∂h(η t )/∂η t first-order derivative of the response function h(η) is evaluated at η t = Z t X t . The score function of l(X T ) in (2.4.11) can be written as

S = (S 1 , • • • , S T ) = Z DΣ -1 {Y -µ}, (2.4.12)
with components

S 0 = R -1 0 (a 0 -X 0 ) (2.4.13) S t = Z t D t Σ -1 t {Y t -µ t } t = 1, • • • , T, (2.4.14)
the weight matrix

W = diag(W 0 , W 1 , • • • , W T ) := DΣ -1 D , (2.4.15)
with diagonal blocks

W 0 = R -1 0 W t = D t Σ -1 t D t , t = 1, • • • , T. (2.4.16)
Denote the first-order derivative of P L(X) in (2.4.11) by

M = ∂P L(X)/∂X = S -KX, (2.4.17)
and the (expected) information matrix of l(X) by

I = -E ∂ 2 P L(X) ∂XX = S + K = Z WZ + K, (2.4.18)
with diagonal blocks

I 0 = R -1 0 (2.4.19) I t = Z t W t Z t , t = 1, • • • , T.
(2.4.20)

A single Fisher scoring to the next iterate X 1 ∈ R m , with m = (T + 1)k is given by

I 0 X 1 -X 0 = M 0 .
where X 0 ∈ R m , is the current iterate, and I 0 , M 0 the first and second derivatives of the current iterate. This can be rewritten as

X 1 = Z W(X 0 )Z + K -1 Z W(X 0 ) Y (X 0 ), (2.4.21) the "working" observation Y(X 0 ) = ( Y 1 (X 0 1 ), • • • , Y T (X 0 T )) can be computed as Y(X 0 ) := D -1 (X 0 ) Y -µ(X 0 )) + ZX 0 , (2.4.22)
with components

Y 0 = a 0 Y t = D -1 t (X 0 t ) Y t -µ(X 0 t ) + Z t X 0 t , (2.4.23)
where µ(X 0 ) = h(ZX 0 ) and X 0 is the value of the state variable of current iterate. Fahrmeir and Wagenpfeil (2013) mentioned one can apply their approach to a linear Gaussian state space model. Then, one can achieve above by assume D is the identity matrix and the score function yields

S = ZΣ -1 (Y -ZX), (2.4.24)
with

Σ = diag(R 0 , Σ 1 , • • • , Σ T ), where Σ t = cov(Y t | X t ).
The weight matrix W t reduces to Σ -1 t , and one uses the actual observations Y instead of the "working" observations Y, since D = I, µ(X 0 ) = ZX 0 , then

a = (Z Σ -1 Z + K) -1 Z Σ -1 Y, (2.4.25)
where a = (a 0|t , a 1|t , • • • , a T |T ) is the vector of smoother estimates. As noted earlier, the classical linear Kalman filter and smoother solves (2.4.25) efficiently, without the expression of the block-tridiagonal matrix (Z Σ -1 Z+K), for more details see Fahrmeir and Tutz (2013), chapter (8).

Working Kalman Filter and Smoother (WKFS)

In the following algorithm, the predictive, filter and smoother values a t|t-1 , a t|t , a t|T respectively are numerical approximations of X t , and V t|t-1 , V t|t and V t|T are numerical approximations of error covariance matrices of predicted, filtered and smoothed values respectively.

Initialization:

a 0|0 = a 0 , V 0|0 = R 0 . (2.4.26) For t = 1, • • • , T :
Prediction:

a t|t-1 = F t a t-1|t-1 V t|t-1 = F t V t-1|t-1 F t + R t (2.4.27)
Filtering:

a t|t = a t|t-1 + K t ( Y t -Z t a t|t-1 ) K t = V t|t-1 Z t {Z t V t|t-1 Z t + (W -1 t ) } -1 V t|t = V t|t-1 -K t Z t V t|t-1
(2.4.28)

Smoothing : For smoothing one uses the classical fixed interval smoother for t = T, • • • , 1 :

a t-1|T = a t-1|t-1 + B t (a t|T -a t|t-1 ) V t-1|T = V t-1|t-1 + B t (V t|T -V t|t-1 )B t where B t = V t-1|t-1 F t V t|t-1 .

Longitudinal Rasch process : Binary data

Bousseboua and Mesbah (2010) introduced a class of longitudinal latent processes where, at any time, a set of categorical binary observations are observed instead of the latent variables. The conditional probabilities are given by Rasch model which is widely used in various psychometric scopes such as educational research and the health analysis, especially in the quality of life. They considered a longitudinal study, where participant patients are interviewed at regular dates of visit. In this interviews, the patients are asked to answer questions in a questionnaire. This questionnaire was constructed to measure their health at the time. It is recognized that the health is a latent multidimensional concept. In practice, each dimension is usually assessed by one or more questions. They focused on the measurement of a single dimension, and in a particular case of dichotomous response options for each item (yes-no, agree-disagree, etc.).

A binary state space model consists of two processes. First, the observed process

Y i (t) = (Y i1 (t), • • • , Y iq (t)) , where the conditional distribution of Y i (t) given the n- dimensional state variable X i (t) is Bernoulli. That is Y i (t) ∼ Bernoulli (π i (t)).
Second, the state process {X i (t)} is assumed to follow a n-dimensional Markov process.

Description of the Rasch latent process

The model defined in Bousseboua and Mesbah (2010) consists in a finite trajectory of a multivariate process:

{(Y ik (t), X i (t)) : 1 ≤ i ≤ n, 1 ≤ k ≤ q, 1 ≤ t ≤ T }
where Y ik (t) represents the longitudinal process of observations. For all i, X i (t) is an unobservable (latent) process. The variable Y ik (t) is the response of the person i at instant t to the item k. This questionnaire is administered on different occasions to the same individuals. These models are characterized by the fact that each these response variable Y it depends only on the corresponding latent trait variable X i (t). The path of the model can be clarified at any time t in Figure (2.2).

{Y 11 (t) • • • • • • Y 1q (t)} Y 1 (t) • • • • • • • • • • • • {Y n1 (t) • • • • • • Y nq (t)} Yn(t) ↑ ↑ X 1 (t) • • • • • • • • • • • • X n (t)
The process of observation

Y i (t) = (Y i1 (t), • • • , Y iq (t)) has values in {0, 1}, where i = 1, • • • , n, t = 1, • • • , T, k = 1, • • • , q.
The value 1 denotes a correct response of the i individual to the item k . Practically, in the education studies, the latent variable X i (t) can be interpreted as a measure of the ability of an individual. Whereas, in the quality of life, in the health field, the latent variable can represent a level of the health of an individual. In general, it is a measure of the individual position in scalar axis corresponding to the latent unidimensional trait measured by the questionnaire. Consequently, the latent variable X i (t) depends only on the individual i. The observation and state equations of the model are written as follows.

The observation equation

They assumed that:

1. The conditional probability of Y ik (t) given X i (t) is Pr(Y ik (t) = y ik (t) | X i (t)) = π ik (t) y ik (t) (1 -π ik (t)) 1-y ik (t) , (2.4.29)
where

Pr(Y ik (t) = 1 | X i (t)) = π ik (t) = exp{y ik (t)(X i (t) -β k )} 1 + exp{X i (t) -β k } , ( 2
X i (1) = x i (1), • • • X i (T ) = x i (T )) : Pr [Y ik (1) = y ik (1), • • • , Y ik (T ) = y ik (T ) | X i (1) = x i (1), • • • , X i (T ) = x i (T ) ] = T t=1 Pr [Y ik (t) = y ik (t) | X i (t) = x i (t) ],
that means the response variables at an instant t is based only on the corresponding latent variable.

3. The response variables Y i1 , • • • , Y ik , are assumed to be conditionally independent to the latent variable X i (t), where for all k = 1,

• • • , q, Y ik = (Y ik (1), • • • , Y ik (T )) : Pr [Y i1 (t) = y i1 (t), • • • , Y ik (t) = y ik (t) | X i (t) = x i (t) ] = q k=1 Pr [Y ik (t) = y ik (t) | X i (t) = x i (t) ].
This means that the q items are conditionally independent relatively to the corresponding latent variable.

State equation

Bousseboua and Mesbah (2010) proposed two classes of latent process: a general firstorder Markov latent process and a first-order autoregressive latent process (AR(1) latent process).

First-order Markov latent process

This model considers the latent process (X i (t) : 1 ≤ t ≤ T ) as being a Markov chain of order one. They assumed the chain (X i (t) : 1 ≤ t ≤ T ) has a real support and a Gaussian distribution with variance σ 2 i . Then:

p(X i (t) | X i (t -1), X i (t -2), ..., X i (1)) = p(X i (t) | X i (t -1)) ∼ N (µ i ; σ 2 i ). (2.4.31)
The vector

X i = (X i (1), • • • , X i (T )
) is gaussian and has probability density g i defined by :

g i (X i ) = 1 σ T i (2π) T exp - 1 2σ 2 i X 2 i (1) + T t=2 (X i (t) -X i (t -1)) 2 (2.4.32)
Autoregressive latent process (AR(1) Process)

{X i (t)} is a stationary gaussian process satisfying a first-order autoregressive AR(1) equation

X i (t) = ρ i X i (t -1) + ε i (t), ε i (t) ∼ N (0; σ 2 i ) (2.4.33) p(X i (t) | X i (t -1), X i (t -2), ..., X i (1)) = p(X i (t) | X i (t -1)) ∼ N (µ i ; σ 2 i 1 -ρ 2 i ). (2.4.34)
The joint distribution of the variables (X i (1), • • • , X i (t)) can be written as :

g i (X i ) = ( 1 -ρ 2 i ) T (2π) T σ T i exp - (1 -ρ 2 i ) 2σ 2 i X 2 i (1) + T t=2 (X i (t) -X i (t -1)) 2 (2.4.35)
where

X i = (X i (1), • • • , X i (T )) and the distribution of X(0) is normal, namely g i (X(0)) ∼ N µ i ; σ 2 i 1-ρ 2 i .

The Marginal Likelihood

The state space models can be considered as a particular case of incomplete data models. The observations Y i (t) can be interpreted as incomplete due to missing state variable X i (t). The likelihood function known as Marginal likelihood can be written :

p(θ i ; Y) = n i=1 • • • p(Y i | X i ; θ)g i (X i )d(X i ). (2.4.36) = n i=1 • • • T t=1 q k=1 π ik (t) y ik (t) (1 -π ik (t)) 1-y ik (t) g i (X i )d(X i ). = n i=1 • • • T t=1 q k=1 [1 -π ik (t)] π ik (t) 1 -π ik (t) y ik (t) g i (X i )d(X i ). = n i=1 • • • T t=1 [1 -π ik (t)] q π ik (t) 1 -π ik (t) r i (t) g i (X i )d(X i ). = n i=1 • • • exp{ T t=1 r i (t).X t -q k=1 β k .r k } T t=1 q k=1 [1 + exp{X t -β k }] .g i (X i )d(X i ), (2.4.37)
where θ i = (β, σ 2 i ) are the parameters for first-order Markov process, or θ i = (β, ρ i , σ 2 i ) for the AR(1) latent process model. Since r i (t) = q k=1 y ik (t) is the score for the iindividual at the t-th occasion, and g i (X) is the joint distribution for first-order Markov chain latent process is described in equation (2.4.32) or joint distribution for AR(1) latent process is described in (2.4.35).

The EM algorithm

The Expectation-Maximization (EM) algorithm was proposed by [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. In fact, EM algorithm is a good approach to the iterative computation of maximum likelihood estimation (MLE). In a variety of incomplete-data problems, the EM algorithm is better than algorithms such as the Newton-Raphson method which may turn out to be more complicated. On each iteration of the EM algorithm, there are two steps-called the Expectation step or the E-step and the Maximization step or the M-step. The cases where the EM algorithm can be performed include not only evidently incomplete-data situations, where there are missing data, truncated distributions, or censored or grouped observations. The term "incomplete data" means there are two samples Y the observed data and X the state variables are not observed, or there are missing data. The EM algorithm aims at computing the MLE of the marginal likelihood by iteratively applying two steps. More precisely, at the step p + 1, to compute θ (p+1) using the value θ (p) which is obtained at the p step. The (p + 1)-th cycle of the EM algorithm consists of the following two steps for p = 0, 1, • • • : Expectation -step: This step computes the conditional mean of the complete loglikelihood knowing the current values of the estimators

Q(θ | θ (p) ) = E log [f (Y, X; θ)] | Y, θ (p) , = n i=1 • • • log{g i (X, θ)} + log{p(Y i , X i ; θ | Y, θ (p) )} ×p(X i | Y i , θ (p) )dX i , (2.4.38)
Bousseboua and Mesbah (2010) calculated this step for Rasch latent Markov processes as follows

Q(θ | θ (p) ) = - T 2 n i=1 log(2πσ 2 i ) + H 1 + H 2 , (2.4.39)
where

H 1 = - n i=1 1 2σ 2 i • • • X i (1) 2 + T t=2 (X i (t) -X i (t -1)) 2 ×p(X i | Y i , θ (p) )d(X i ). H 2 = + n i=1 q k=1 T t=1 • • • log exp{Y ik (t)(X i (t) -β k } 1 + exp{X i (t) -β k } ×p(X i | Y i , θ (p) )d(X i ).
Here, P (X i | Y i , θ (p) ) is the conditional density of the latent vector

X i = (X i (1), • • • , X i (T ))
given the observation vector Y i .

Bousseboua and Mesbah (2010) approximated this density by:

p(X i | Y i , θ (p) ) ∝ exp T t=1 X i (t).Y ik (t) -q k=1 β p k .Y ik (t) -C T t=1 q k=1 [1 + exp(X i (t) -β p k ] , (2.4.40)
where

C = 1 2σ 2(p) i X 2 i (1) + T t=2 (X i (t) -X i (t -1)) 2 .
The maximization with respect to σ 2 i s yields :

σ 2(p+1) i = 1 T • • • X 2 i (1) + T t=2 (X i (t) -X i (t -1)) 2 ×p(X i | Y i , θ (p) )d(X i ). (2.4.41)
Maximizing with respect to β includes only the double sum over i and t of the last terms of (2.4.39). That requires solving the q following equations,

for i = 1, • • • , n n i=1 T t=1 Y ik (t) = n i=1 T t=1 • • • [1 + exp(X i (t) -β k )] -1 ×p(X i | Y i , θ (p) )d(X i ), k = 1, • • • , q. (2.4.42)
On other hand, the estimation of Rasch latent AR(1) processes is as follows. Denote the parameters of the model by θ

= (β, ρ, σ 2 ) with β = (β 1 , • • • , β k ) , ρ = (ρ 1 , • • • , ρ n ) and σ 2 = (σ 2 1 , • • • , σ 2 n ) , ∀i, σ 2 i = 0.
Then, E-step is given by:

Q(θ | θ (p) ) = - T 2 n i=1 log(2πσ 2 i ) + 1 2 n i=1 log(1 -ρ 2 i ) + G 1 + G 2 , (2.4.43)
where

G 1 = - n i=1 1 2σ 2 i • • • (1 -ρ 2 i )X 2 1 + T t=1 (X t -ρ i X i (t -1)) 2 ×p(X i | Y i , θ (p) )d(X i ). G 2 = + n i=1 q k=1 T t=1 • • • log exp{Y ik (t)(X i (t) -β k } 1 + exp{X i (t) -β k } ×p(X i | Y i , θ (p) )d(X i ). and p(X i | Y i , θ (p)
) is the conditional density of latent vector

X i = (X i (1), • • • , X i (T )) given Y i .
Bousseboua and Mesbah (2010) approximated this density by:

p(X i | Y i , θ (p) ) ∝ exp T t=1 X i (t).Y ik (t) -q k=1 β p k .Y ik (t)-C T t=1 q k=1 [1 + exp(X i (t) -β p k )]
(2.4.44)

C = 1 2σ 2(p) i (1 -ρ 2 (p) i )X 2 i (1) + T t=2 (X i (t) -ρ (p) i X i (t -1)) 2 .
Here, the maximization with respect to σ 2 i also includes only the first and third terms in (2.4.43), and maximizing with respect to ρ i includes the second and third term of this expression. One obtains :

σ 2(p+1) i = 1 T n i=1 • • • (1 -ρ 2 i )X 2 i (1) + T t=2 (X i (t) -ρ i X i (t -1) 2 ×p(X i | Y i , θ (p) )d(X i ). i = 1, • • • , n, (2.4.45)
ρ i 1 -ρ 2 i = 1 σ 2 i • • • ρ i X 2 i (1) + T t=2 X i (t -1)(X i (t) -ρ i X i (t -1)) ×p(X i | Y i , θ (p) )d(X i ). (2.4.46)
and maximizing with respect to β includes only the double sum over i and t of the last term of the expression (2.4.43). For i = 1, • • • , n, that leads to solving the q following equations :

n i=1 T t=1 Y ik (t) = n i=1 T t=1 • • • [1 + exp(X i (t) -β k )] -1 (2.4.47) ×p(X i | Y i , θ (p) )d(X i ). k = 1, • • • , q.
Maximization-step: The second step aims to find the value θ (p+1) that achieves the maximum of the quantity Q(θ | θ (p) ) :

Q(θ | θ (p+1) ) = arg max θ Q(θ | θ (p) ) (2.4.48)
The iterative of the E and M stage are concluded if

| Q(θ (p+1) ) -Q(θ (p) ) |< ,
where is a prior fixed quantity. Bousseboua and Mesbah (2010) applied the Newton Raphson method to calculate the value θ (p+1) , and the integral are approximated numerically using the Gauss Hermite quadrature formulas.

Summary of the convergence properties of the EM algorithm

The convergence properties of the EM algorithm was studied by [START_REF] Wu | On the convergence properties of the EM algorithm[END_REF]. He proved the EM sequence converges to the unique maximum likelihood estimate. He specified the following two special cases:

(a) The unobserved complete-data density function can be described by a curved exponential family distribution with parameters having compact space. The likelihood is unimodal when it has a unique mode.

(b) The likelihood function is unimodal and can undoubtedly satisfy a differentiability condition.

These properties of the EM algorithm help us to prove the asymptotic normality of the parameters. The properties of the EM algorithm can be briefly presented as follows:

(i) The likelihood and L(θ (p) ) increases with any EM sequence {θ (p) }. If it is bounded, converges to some L * , where L * is a stationary value of L.

(ii) The continuity of Q achieved for every density function has a curved exponential

family if Q(θ | θ (p) ) is continuous in both θ (p) and θ. If θ (p) converges to some point θ * , θ * is a stationary point under the continuity condition of DQ(θ | θ (p) ) in θ (p) and θ. DQ(θ | θ (p)
) is first-order derivative with respect to the θ (p) , that is

DQ(θ | θ (p) ) = ∂ ∂θ Q(θ | θ (p) ) | θ=θ (p) .
(iii) If, in addition to (ii), Q is not trapped at any a stationary point θ 0 , but not a local maximum of L, i.e.

sup θ∈Θ Q(θ | θ 0 ) > Q(θ 0 | θ 0 ),
consequently, L * is also a local maximum of L.

(iv) If, in addition to (ii) or (iii), θ (p+1) -θ (p) → 0,1 as p → ∞ and the stationary point (local maxima) set with a L took a discrete value, thereafter θ (p) converges to a stationary point (local maximum).

(v) If, in addition to (ii) or (iii), unacceptable, there exist two different stationary points (local maxima) relates to the same L value, then θ (p) converges to a stationary point (local maximum).

(vi) θ (p) converges to the unique maximizer θ * of L(θ), if L(θ) achieves that it is unimodal in Θ and has only stationary point and D(θ | θ (p) ) is continuous in θ and θ (p) ,

Consistency and Asymptotic normality of the maximum likelihood estimator (MLE)

The inferential methods for any model must achieve the asymptotic properties of ML estimators. Under regularity assumptions which are determined by the proposed to the model, the following asymptotic properties of the estimators θ for the generalized state space model

Asymptotic existence and uniqueness: The estimators exist and are (locally) unique as a sample size n → +∞.

Consistency: If θ refer to the "true" parameter value, then as n → ∞, then θ → P θ in probability (weak consistency) or with probability p > 0.5 (strong consistency).

Asymptotic normality: The distribution of the MLE is normal for n → ∞, namely, θ ∼ a.s N (θ, I(θ))

and n 1/2 ( θ -θ) ∼ N (0, I -1 (θ)).

i.e., θ approximately is normal with approximate (or "asymptotic") covariance matrix cov( θ) = I(θ),

where I -1 (θ) is the inverse of the Fisher matrix and it can be calculated. Furthermore, the MLE is asymptotically efficient according to a broad class of other estimators.

Oakes' identity to find the information matrix via the EM algorithm

Oakes (1999) assumed that the observed data Y with likelihood L(φ, Y ) depending on the parameter vector φ may be specified as a many-to-one parameter-free map 2 of full data X with log-likelihood L(φ, Y ). The objective is to maximize L(φ, Y ) in φ.

Q(φ p | φ) = E X|Y ;φ L 0 (φ p , x).
The fundamental identity

L(φ p , X) = Q(φ p | φ) -E X|Y ;φ log{p(X | Y ; φ p )}, (2.5.1)
φ p denoted the result from the iteration p of the algorithm. The EM algorithm yields iteratively, at each step selecting φ p to maximize (φ p | φ) in its first argument φ p with its second argument φ held fixed. The procedure will be useful when as is often the case, the functional form of L(φ p , x) is appreciably simpler than that of L(φ, y). If L(φ p , X) is of exponential family form the maximization can be divided into a separate "E-step" -calculate the conditional expectation of the sufficient statistics in X given the observed data Y , using the current parameter estimate φ-and M-step -maximize the log-likelihood L(φ p , x) with these computed values for the sufficient statistics and choose the new parameter value φ p .

2.5.1 Derivative of Q(φ p | φ) and L(φ)

Oakes (1999) assumed that the procedures of expectation usually interchange with respect to X and differentiation in φ hold for log{p(X | Y ; φ p )}. Thus

E X|Y ;φ ∂[log{p(X | Y ; φ p )}] ∂φ = 0 (2.5.2)
and

E X|Y ;φ ∂ 2 [log{p(X | Y ; φ p )} ∂φ 2 = -E X|Y ;φ ∂[log{p(X | Y ; φ p )}] ∂φ ∂[log{p(X | Y ; φ p )}] T ∂φ . (2.5.3)
Differentiation of equation (2.5.1) in φ p gives

∂L ∂φ p = ∂Q(φ p | φ) ∂φ p -E X|Y ;φ ∂[log{p(X | Y ; φ p )}] ∂φ p (2.5.4)
The substitution φ = φ p makes the last term disappear by equation (2.5.2). So the score statistic for the observed data yields

∂L ∂φ = ∂Q(φ p | φ) ∂φ p φ p =φ . (2.5.5)
Differentiation of equation (2.5.4) in φ p and φ gives respectively

∂ 2 L ∂φ p 2 = ∂ 2 Q(φ p | φ) ∂φ p 2 -E X|Y ;φ ∂ 2 [log{p(X | Y ; φ p )}] ∂φ p 2 , 0 = ∂ 2 Q(φ p | φ) ∂φ∂φ p -E X|Y ;φ ∂[log{p(X | Y ; φ p )}] ∂φ p ∂[log{p(X | Y ; φ p )}] ∂φ .
Substituting φ = φ p , adding the two equations and using equation (2.5.3) yields

∂ 2 L ∂φ 2 = ∂ 2 Q(φ p | φ) ∂φ p 2 + ∂ 2 Q(φ p | φ) ∂φ p ∂φ φ p =φ , (2.5.6)
which is valid for all φ. The second term in equation (2.5.6) is naturally called the "missing information",i.e. due to the fact that only Y and not X is observed.

Chapter 3

Longitudinal multicategorical processes : Generalized state space models with Gaussian noises

Introduction

In this chapter, we introduce a new class of longitudinal multivariate processes which are multicategorical. The data are from a longitudinal study, where the patients participate in an interview. The interview aims at measuring the patients' health at regular intervals, the dates of which are determined before the study. It typically involves filling out a questionnaire in which they are asked multiple choice questions. These questionnaires are used in the quality of life studies in health scope, constructed in order to measure the patient's perceived health at the time of the visit. It is recognized that health is a latent multidimensional concept. We are interested to latent variables X i (t) produced by an individual i,

(i = 1, • • • , n), at time t, (t = 1, • • • , T ).
The X i (t) s may be the patients health, a latent trait, etc.. We only observe Y i (t) instead of X i (t). The Y i (t) s are the responses of the individuals to the questionnaire where each item in the questionnaire has multiple options, the individual selects one of them.

In the past studies, the latent variables X i (t) s were described by first-order autoregressive model AR(1) with gaussian noise. In this chapter they are described by a firstorder conditional heteroscedastic non-linear (CHARN) models with gaussian noise. In the state space models whether linear or non-linear if the observations and the latent variable have a gaussian distribution, the posterior distribution also has a gaussian distribution [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF], [START_REF] Kitagawa | Introduction to time series modeling[END_REF]). Consequently, the Kalman filtering is used to estimate the posterior mean and variance-covariance matrix. In the present context, the observations are from a multinomial distribution and the latent variables from a gaussian distribution. Since the posterior distribution is not symmetric we con-sider the posterior mode as an estimator of the latent variable. An alternative recording might be: 1 (awake) , 2 (active sleep), 3 (indeterminate sleep), 4 (quiet sleep) However, the diverse records leads to different results. In order to reduce the risk of multiple interpretations, the categorical time series is rewritten (regardless of the measurement scale) as the vector

U t = (U t1 , • • • , U ts ) of length s = c-1, with elements U tj = 1,
if the category j is observed at time t 0, with otherwise. The first recording of the example is written as follows:

-U t = (1, 0, 0) indicate "quiet sleep" -U t = (0, 1, 0) indicate "indeterminate sleep " -U t = (0, 0, 1)
indicate " active sleep" -U t = (0, 0, 0) indicate " awake".

After this explanation, we present our model which assumes that n individuals are investigated. The model consists of a finite trajectory of a multivariate process:

{(Y ik (t), X i (t)) : 1 ≤ i ≤ n, 1 ≤ k ≤ q, 1 ≤ t ≤ T } Figure 3.1: The Path of longitudinal multicategorical processes {Y 1 11 (t), • • • , Y c 1 11 (t)} Y 11 (t) , • • • , {Y 1 1k (t), • • • , Y c k 1k (t)} Y 1k (t) , • • • , {Y 1 1q (t), • • • , Y cq 1q (t)} Y 1q (t) • • • • • • {Y 1 n1 (t), • • • , Y c 1 n1 (t)} Y n1 (t) , • • • , {Y 1 nk (t), • • • , Y c k nk (t)} Y nk (t) , • • • ,{Y 1 nq (t), • • • , Y cq nq (t)} Ynq (t)     Item 1     Item k     Item q     Item 1     Item k     Item q Y 1 (t) • • • • • • Yn (t) ↑ ↑ X 1 (t) X n (t) • • • • • • • • • • • • • • • • • •
where for all i = 1,

• • • , n, k = 1, • • • , q (Y ik (t) = (Y 1 ik (t), • • • , Y c k ik (t))
) is the vector of longitudinal observations and, for all i, (X i (t)) is an unobservable (latent) process. The variable Y s ik (t), s = 1, • • • , c k , denotes the response to category s at instant t, of the person i to the item k. These model designs the q × (c

1 + c 2 + • • • + c q ) vector responses (Y i1 (t), Y i2 (t), • • • , Y iq (t)
) at every occasion t, the responses of individual i are categorical for q-items. In other words, every item has a vector of categories responses. This questionnaire is administered on different occasions to the same individuals. A characteristic of the models considered is that each response vector Y ik (t) depends only on the corresponding latent trait X i (t), as shown in figure (3.1). In this context , the model can be represented at any time t by the processes of observations {Y s ik (t)} which have values 0 or 1. The value 1 indicates that the category s has been observed, and 0 if it is not, s = 1, • • • , c k , with c k denoting the categories number of item k. For each individual at instant t one has q(c 1 + • • • + c q ) vectors of responses with values 0 and 1.

Modeling

In this section, we construct the observation and state equations for our model. The marginal likelihood is written. We compute the parameters estimation of the model, and study their consistency and asymptotic normality. One denotes the conditional probability π by

Pr[Y ik (t) = (y 1 ik (t), • • • , y c k ik (t)) | X i (t) = x i (t)] = π((y 1 ik (t), • • • , y c k ik (t)) | X i (t)), i = 1, • • • , n, k = 1, • • • , q, t = 1, • • • , T,

The observation equation

We make the following assumptions:

1. The conditional probability of Y ik (t) given X i (t) is a multinomial distribution.

More explicitly, for all

i = 1, • • • , n, k = 1, • • • , q, t = 1, • • • , T, Pr[Y ik (t) = (y 1 ik (t), • • • y c k ik (t)) | X i (t) = x i (t)] = c k s=1 [π s ik (t)] y s ik (t) , (3.3.1)
where

π s ik (t) = exp[η s ik (t)] 1 + c k j=1 exp[η j ik (t)] if s < c k , π c k ik (t) = 1 1 + c k j=1 exp[η j ik (t)] , (3.3.2)
and c k s=1 y s ik (t) = 1, c k s=1 π s ik (t) = 1. The link function η s ik (t) is defined with the logit function as follows

η s ik (t) = logit(π s ik (t)) = log π s ik (t) π c k ik (t) = log π s ik (t) 1 -c k -1 j=1 π j ik (t) = u i (t)β s k + X i (t),
and

η c k ik (t) = log π c k ik (t) π c k ik (t) = log(1) = 0. u i (t) = (u i1 (t), • • • , u ir (t)
) is the vector of independent covariate variables and r is their number. For k = 1, • • • , q, the β s k = (β s k1 , • • • , β s kr ) are vectors of unknown regression parameters.

The vectors Y

ik (1), • • • , Y ik (T ) are conditionally independent given the vector of latent variable X i (1) = x i (1), • • • X i (T ) = x i (T ) : Pr[Y ik (1) = y ik (1), • • • , Y ik (T ) = y ik (T ) | X i (1) = x i (1), • • • , X i (T ) = x i (T )] = T t=1 Pr[Y ik (t) = y ik (t) | X i (t) = x i (t)] 3. The vectors Y i1 (t), • • • , Y ik (t) are conditionally independent given the latent vari- able X i (t) : Pr[Y i1 (t) = y i1 (t), • • • , Y ik (t) = y ik (t) | X i (t) = x i (t)] = q k=1 Pr[Y ik (t) = y ik (t) | X i (t) = x i (t)]

The state equation

We first give some examples of usual nonlinear time series models.

The nonlinear time series models

Tjøstheim (1994) and Turkman et al.(2014) defined nonlinear time series models that broadly are classified into the following categories:

1. Parametric models for the conditional mean In these models appears the conditional mean function of the process X t as a nonlinear function of the past observations. In other hand, the conditional variance is constant. The general model is written as :

X t = F (X t-1 , θ) + ε t ,
where F is a known nonlinear function, (ε t ) is the noise process and θ is the unknown parameter vector.

E[X i (t) | X i (t -1) = x] = F (x, θ).
Here are some examples of these models :

(a) Exponential AR:

F (x, θ) = {θ 1 + θ 2 ) exp(-θ 3 x 2 )}x, θ = (θ 1 , θ 2 , θ 3 ) (b) Logistic AR: F (x, θ) = θ 1 x+θ 2 x{[1+exp(-θ 3 (x-θ 4 ))] -1 -1/2}, θ 3 > 0, θ = (θ 1 , θ 2 , θ 3 , θ 4 )

Parametric models for the conditional variance

In these models, the conditional variance function of the process X t is a nonlinear function of the past observations. Whereas the conditional mean is constant, the general model is written as :

X t = H(X t-1 , θ)ε t
This model is widely used in financial applications as ARCH and GARCH models. ARCH(1) model has the form

X t = (θ 1 + θ 2 X 2 t-1 ) 1/2 ε t , where Var[X t | X t-1 = x] = (θ 1 + θ 2 x 2 t-1 )σ 2 ε , with Var(ε t ) = σ 2 ε .

Mixed parametric models for the conditional mean and variance

In these models both the conditional mean and variance functions of the process X t appear as nonlinear functions of the past observations. A general model containing both a conditional mean and conditional variance component is the conditional heteroscedastic autoregressive nonlinear (CHARN) model, written on the form

X t = F (X t-1 , θ 1 ) + H(X t-1 , θ 2 )ε t .
This model contains special bilinear models of the form:

X t = θ 1 + θ 2 X t-1 + θ 3 X t-1 ε t + ε t .
The above example includes nonlinear dynamics both in the mean and the variance, since E

(X t | X t-1 ) = θ 1 + θ 2 X t-1 and Var[X t | X t-1 ] = (1 + θ 3 X t-1 ) 2 σ 2 ε .
Some properties of such models are studied in Tjøstheim (1986), Tjøstheim (1994), Brockwell and Davis (1991), Brockwell and Davis( 2013), Shumway and Stoffer(2001), and [START_REF] Tong | Non-linear Time series. A Dynamic System Approach[END_REF]. The parameter estimation is discussed in Ngatchou-Wandji (2008) among others.

First-order CHARN latent process

We describe the state equation by first-order CHARN model :

X i (t) = F [X i (t -1), u i (t), γ] + H[X i (t -1), u i (t), δ]ε i (t), (3.3.3)
where:

• γ, δ are the model parameters.

• (ε i (t)) is the gaussian noise process for the state process :

ε i (t) ∼ N (0; R t ), R t > 0.
• (u i (t)) is the covariate variable process, u i (t) ∈ R r , we recall r is the number of covariate variable process.

• F (., ., .) : R × R r × R l -→ R, is a non-linear function.

• H(., ., .) : R × R r × R l -→ R, is a non-linear function.

Remark: From simple computations, if one assumes that the ε i 's have common density function f, then the conditional distribution of X i (t) given X i (t -1) = x is given by

f (i) X i (t-1)=x (z) = 1 H(x, u i (t), δ) f z -F (x, u i (t), γ) H(x, u i (t), δ) = 1 2πR 2 t H(x, u i (t), δ) exp -(z -F (x, u i (t), γ)) 2 2R t H 2 (x, u i (t), δ) , (3.3.4)
From the above remark, it's clear that (X i (t) :

1 ≤ t ≤ T ) satisfies p(X i (t) | X i (t -1), X i (t -2), • • • , X i (1)) = p(X i (t) | X i (t -1)) ∼ N (µ i (t), V i (t)),
where for all i = 1, • • • , n, t = 1, • • • , T, the conditional mean µ i (t) and the conditional variance V i (t) are defined by

µ i (t) = E[X i (t) | X i (t -1), u i (t)] = F [X i (t -1), u i (t), γ] V i (t) = Var[X i (t) | X i (t -1), u i (t)] = H 2 [X i (t -1), u i (t), δ]R t
Then, the joint law

g i of X i = (X i (0), X i (1), • • • , X i (T )
) is obtained easily by successive conditionings:

g i (X i ) = T t=1 p(X i (t) | X i (t -1)) × p(X i (0)) = T t=0 N (µ i (t), V i (t)), = 1 (2π) T T t=0 V i (t) exp - T t=0 [X i (t) -µ i (t)] 2 2V i (t) . (3.3.5)

The marginal likelihood

For all i = 1, • • • , n, let Y i = (Y i (0), Y i (1), • • • , Y i (T )) , Y i (t) = (Y i1 (t) • • • , Y iq (t)) , t = 0, 1, 2, • • • , T. with Y ik (t) = (Y (1) ik (t), • • • , Y (c k ) ik (t)) , t = 0, 1, 2, • • • , T, k = 1, • • • , q, and
X i = (X i (0), X i (1), • • • , X i (T )) , d(X i ) = (d(X i (0)), d(X i (1)), • • • , d(X i (T ))) .
Let θ = (β, γ, δ) with β a q-dimensional vectors β = (β 1 , • • • , β q ) , with β k s, being c k × r matrices, where c k is categories number of item k, k = 1, • • • , q, and r denotes the number of covariates. Let y i = (y i (0),

y i (1), • • • , y i (T )) with y i (t) = (y i1 (t) • • • , y iq (t)) , define Pr(Y 1 = y 1 , Y 2 = y 2 , • • • , Y n = y n ) = n i=1 • • • p(y i | X i ; θ)g i (X i )d(X i ) = n i=1 • • • T t=0 q k=1 Pr(Y ik (t) = (y 1 ik (t), • • • , y c k ik (t)) | X i (t))g i (X i )d(X i ) = n i=1 • • • T t=0 q k=1 c k s=1 [π s ik (t)] y s ik (t) × g i (X i )d(X i ). = n i=1 • • • T t=0 q k=1 c k s=1 exp[η s ik (t)] 1 + c k j=1 exp[η j ik (t)] y s ik (t) × g i (X i )d(X i ). (3.3.6)
where g i (X i ) is the joint density function of the latent variables X i defined by equation (3.3.5). Then the likelihood is given by

p(Y 1 , Y 2 , • • • , Y n ) = n i=1 • • • T t=0 q k=1 c k s=1 exp[u i (t)β s k + X i (t))] 1 + c k j=1 exp[u i (t)β j k + X i (t)] Y s ik (t) ×g i (X i )d(X i ), (3.3.7)

The EM algorithm

The EM algorithm is a general iterative method to obtain maximum likelihood estimators in incomplete data cases. If θ (0) denotes a starting value for θ, the (p + 1)-th cycle of the EM algorithm consists of the following two steps for p = 0, 1, • • • :

Expectation -step: compute the expectation Q(θ | θ (p) ) as follows Q(θ | θ (p) ) = E{log[f (Y, X; θ)] | Y, θ (p) } = n i=1 • • • [log{g i (X i , θ i ) + log{p(Y i | X i )}] ×p(X i | Y i , θ (p) )d(X i ).
The E-step for our model yields :

Q(θ | θ (p) ) = - nT 2 log(2π) + G 1 + G 2 (3.3.8)
where

G 1 = - 1 2 n i=1 • • • T t=0 log V i (t) + [X i (t) -µ i (t)] 2 V i (t) ×p(X i | Y i , θ (p) )d(X i ). (3.3.9) G 2 = n i=1 q k=1 T t=0 c k s=1 Y s ik (t) • • • log[π s ik (t)] ×p(X i | Y i , θ (p) )d(X i ), (3.3.10) and p(X i | Y i , θ (p)
) is the conditional density of the latent vectors X i given the observation vectors Y i which is later computed by the particle filtering algorithm.

Maximizing step:

θ (p+1) = arg max θ Q(θ | θ (p) )

Estimation of first-order CHARN latent processes

We recall that the parameter vector is θ

= (β 1 , • • • , β q , γ , δ ) , k = 1, • • • , q, and β k = (β 1 k , • • • , β c k k )
. By applying the E-M algorithm MLE can be found as follows : maximize with respect to the β, only the part G 2 , β s(p+1) k is the solution of the following equation:

n i=1 T t=0 • • • u i (t)D ik (X t )Σ -1 ik (X t )[Y s ik (t) -π s ik (t)] × p(X i | Y i , θ (p) )d(X i ) = 0,
where

π s ik (t) = exp[u i (t)β s k + X i (t)] 1 + c k j=1 exp[u i (t)β j k + X i (t)]
.

The derivation of β s k is given in Appendix (B.1). Where

D s ik (X t ) = ∂π s ik (t) ∂η s ik = ∂{exp[η s ik (t)]/[1 + c k j=1 exp[η j ik (t)]} ∂η s ik (t) = (1 + c k j=1 exp[η j ik (t)]).(exp[η s ik (t)]) -(exp[η s ik (t)]) 2 (1 + c k j=1 exp[η j ik (t)]) 2 = exp[η s ik (t)]{(1 + c k j=1 exp[η j ik (t)]) -(exp[η s ik (t)])} (1 + c k j=1 exp[η j ik (t)]) 2 = exp[η s ik (t)] 1 + c k j=1 exp[η j ik (t)] × (1 + c k j=1 exp[η j ik (t)]) -exp[η s ik (t)] 1 + c k j=1 exp[η j ik (t)] = π s ik (t) × 1 + c k j=1 exp[η j ik (t)] 1 + c k j=1 exp[η j ik (t)] - exp[η s ik (t)] 1 + c k j=1 exp[η j ik (t)] = π s ik (t) × 1 - exp[η s ik (t)] 1 + c k j=1 exp[η j ik (t)] = π s ik (t)[1 -π s ik (t)] Σ ik (X t ) = cov(Y ik (t)) has generic elements σ sm ik (t) = π s ik (t)[1 -π s ik (t)], if s = m -π s ik (t)π m ik (t) if s = m
Maximizing the term G 1 with respect to γ, γ (p+1) is the solution of the following equation:

n i=1 • • • T t=0 [X i (t) -µ i (t)] V i (t) ∂µ i (t) ∂γ × p(X i | Y i , θ (p) )d(X i ) = 0.
Maximizing the term G 1 with respect to δ, δ (p+1) is the solution of the following equation:

1 2 n i=1 • • • T t=0 -1 V i (t) + [X i (t) -µ i (t)] 2 V 2 i (t) ∂V i (t) ∂δ × p(X i | Y i , θ (p) )d(X i ) = 0. (3.3.11)

Information matrix via the EM algorithm

In this subsection, we compute the observed information matrix of our model by implementing Oakes' identity (1999) presented in chapter 2. Oakes (1999) derived the following identity that involves two components:

I(θ) = - ∂ 2 Q(θ | θ (p) ) ∂θ 2 | θ (p) =θ + ∂ 2 Q(θ | θ (p) ) ∂θ (p) ∂θ | θ (p) =θ . (3.3.12)
The first component is directly computed by the EM algorithm. This component is the second-order derivative of the conditional expected value of the complete data loglikelihood given the observed data. The first component in (3.3.12) the block-diagonal matrix1 is defined as follows:

∂ 2 Q(θ | θ (p) ) ∂θ 2 = diag ∂ 2 Q(β | β (p) ) ∂β 2 , ∂ 2 Q(γ | γ (p) ) ∂γ 2 , ∂ 2 Q(δ | δ (p) ) ∂δ 2
The second component in (3.3.12) is the first-order derivative of the score, for the same expected log-likelihood, with respect to the current value of the parameters. This component known as the "missing information", is defined as :

∂ 2 Q(θ | θ (p) ) ∂θ (p) ∂θ = ∂ 2 Q(β | β (p) ) ∂θ (p) ∂β , ∂ 2 Q(γ | γ (p) ) ∂θ (p) ∂γ , ∂ 2 Q(δ | δ (p) ) ∂θ (p) ∂δ
First component : Second-order derivative with respect to β s k yields :

∂ 2 Q(β | β (p) ) ∂β s2 k = - n i=1 T t=0 • • • u i (t)D ik (X t )Σ -1 ik (X t )D ik (X t )u i (t) ×p(X i | Y i , θ (p) )d(X i ). (3.3.13)
second-order derivative with respect to γ and δ as follows:

∂ 2 Q(γ | γ (p) ) ∂γ 2 = n i=1 • • • ∂ 2 µ i (t)/∂γ 2 V i (t) × p(X i | Y i , θ (p) )d(X i ). = n i=1 • • • ∂ 2 ∂γ 2 (F [X i (t -1), u i (t), γ]) H 2 (x t-1 , u i (t), δ)R 2 ×p(X i | Y i , θ (p) )d(X i ). (3.3.14) ∂ 2 Q(δ | δ (p) ) ∂δ 2 = 1 2 n i=1 • • • T t=0 1 V 2 i (t) -2 [X i (t) -µ i (t)] 2 V 3 i (t) ∂ 2 V i (t) ∂δ 2 ×p(X i | Y i , θ (p) )d(X i ). ( 3
∂ 2 Q(β | β (p) ) ∂θ (p) ∂β s k = n i=1 T t=0 • • • u i (t)D ik (X t )Σ -1 ik (X t )[Y s ik (t) -π s ik (t)] × ∂p(X i | Y i , θ (p) )θ (p) ) ∂θ (p) d(X i ). (3.3.16) ∂ 2 Q(γ | γ (p) ) ∂θ (p) ∂γ = n i=1 • • • T t=0 [X i (t) -µ i (t)] V i (t) ∂µ i (t) ∂γ × ∂p(X i | Y i , θ (p) )θ (p) ) ∂θ (p) d(X i ). (3.3.17)
and

∂ 2 Q(δ | δ (p) ) ∂θ (p) ∂δ = 1 2 n i=1 • • • T t=0 -1 V i (t) + [X i (t) -µ i (t)] 2 V 2 i (t) ∂V i (t) ∂δ × ∂p(X i | Y i , θ (p) )θ (p) ) ∂θ (p) d(X i ). (3.3.18)
The derivation of the posterior density is obtained by using a numerical methods to approximate the derivative via the particle filter proposed by [START_REF] Poyiadjis | Maximum likelihood parameter estimation in general state-space models using particle methods[END_REF].

The formula of derivation of the posterior density is presented later in the section of the posterior distribution.

Fisher information :

The Fisher information matrix for each parameter can be written as follows :

• The Fisher information matrix of parameter β as follows :

I(β) = - ∂ 2 Q(β | β (p) ) ∂β s2 k + ∂ 2 Q(β | β (p) ) ∂θ (p) ∂β s k , (3.3.19)
where

∂ 2 Q(β|β (p) ) ∂β s2 k and ∂ 2 Q(β|β (p) )
∂θ (p) ∂β s k are obtained by equations ( 3.3.13 ) and (3.3.16), respectively.

• The Fisher information matrix of parameter γ is :

I(γ) = - ∂ 2 Q(γ | γ (p) ) ∂γ 2 + ∂ 2 Q(γ | γ (p) ) ∂θ (p) ∂γ , (3.3.20)
where

∂ 2 Q(γ|γ (p) ) ∂γ 2 and ∂ 2 Q(γ|γ (p) )
∂θ (p) ∂γ are computed by equations ( 3.3.14) and (3.3.17), respectively.

• The Fisher information matrix of parameter δ is given by :

I(δ) = - ∂ 2 Q(δ | δ (p) ) ∂δ 2 + ∂ 2 Q(δ | δ (p) ) ∂θ (p) ∂δ , (3.3.21)
where

∂ 2 Q(δ|δ (p) ) ∂δ 2 and ∂ 2 Q(δ|δ (p) ) ∂θ (p) ∂δ
are given by equations ( 3.3.15) and (3.3.18), respectively.

Assumptions

In the literature to have the consistency of MLE's, the underlying density (likelihood function) must satisfy certain "regularity conditions" about the sample Y i and the parameters θ. For more details, see Lehmann and Casella (1998, section 6.3). Under the following assumptions, the asymptotic properties of MLE's of our model can be established:

(A.1) For all i = 1, • • • , n, the observations Y i = (Y i (0), Y i (1), • • • , Y i (T )) , are
such that for all t = 1, • • • , T and k = 1, • • • , q, the conditional probability Pr(Y ik (t) | X i ) is described by a multinomial distribution.

(A.2) For all i = 1, • • • , n and k = 1, • • • , q, the vectors Y ik (1), • • • , Y ik (T ) are conditionally independent given the vector of latent variable

X i (1) = x i (1), • • • X i (T ) = x i (T ) : Pr(Y ik (1) = y ik (1), • • • , Y ik (T ) = y ik (T ) | X i (1) = x i (1), • • • , X i (T ) = x i (T )) = T t=1 Pr(Y ik (t) = y ik (t) | X i (t) = x i (t))
where Pr(Y ik (1) = y ik (1),

• • • , Y ik (T ) = y ik (T ) | X i (1) = x i (1), • • • , X i (T ) = x i (T )) is the joint mass function of Y i conditional to X i . (A.3) For all i = 1, • • • , n and t = 1, • • • , T, the vectors Y i1 (t), • • • , Y iq (t)
are conditionally independent knowing the latent variable X i (t) :

Pr(Y i1 (t) = y i1 (t), • • • , Y ik (t) = y ik (t) | X i (t) = x i (t)) = q k=1 Pr(Y ik (t) = y ik (t) | X i (t) = x i (t))
where Pr(Y i1 (t) = y i1 (t),

• • • , Y ik (t) = y ik (t) | X i (t) = x i (t)
) is the joint mass function conditional to X i (t).

(A.4) The process X i (t) follows equation (4.2.3) with noise process having Gaussian density function.

(A.5) The distribution of the observation has common support. This assumption holds since Y ik (t) = (y 1 ik (t), • • • , y c k ik (t)) where y s ik (t) has a value 0 or 1

, i = 1, • • • , n, k = 1, • • • , q, t = 1, • • • , T, s = 1, • • • , c k . (A.6) The parameters are identifiable; that is, if θ = θ , then the marginal likelihood satisfies p(Y 1 , • • • , Y n , θ) = p(Y 1 , • • • , Y n , θ ).
(A.7) The parameter space Θ contains an open ball H such that the true parameter value θ 0 is an interior point.

For the vector of parameters θ = (β, γ, δ), as n → +∞, and T → +∞, we give a short outline of standard (N ) 1/2 -asymptotic where N = (n × T ). For this case typical "regularity assumptions " (A.1-A.7) are weak conditions, in particular convergence of I(θ N )/N = cov( θN )/N, say as N → +∞,

I(θ N )/N → I(θ). (3.3.22)
We have under regularity assumptions (A.1-A.7)

N -1/2 L (θ) → N (0, I(θ)), (3.3.23)
where L (θ) = s(θ) is first-order derivative of the likelihood function with respect to θ, and is called the score function. Then the following Theorem hold under regularity assumptions (A.1-A.7) and the convergence properties to the EM algorithm

Theorem Let Y i = (Y i (0), Y i (1), • • • , Y i (T )) , i = 1, • • • , n, and N = T × n.
Let θ be the MLE's of θ. Under the regularity assumptions (A.1-A.7) and convergence properties of EM algorithm θ is asymptotically normal :

N 1/2 ( θ -θ) ∼ N (0; I -1 (θ)). ( 3

.3.24)

Proof: The likelihood function L(θ; Y ) depends on θ. Denote first-order derivative of the log likelihood (with respect to θ) by L , and second-order, and third-order derivative by L and L . Now, by using the Taylor expansion for L ( θ) in a neighbourhood of θ.

For some θ between θ and θ, one has

L ( θ) = L (θ) + ( θ -θ)L (θ) + 1 2 ( θ -θ) 2 L ( θ) = L (θ) + ( θ -θ)L (θ) + O(| θ -θ | 2 ). ( 3

.3.25)

O is used for describing the limiting behaviour of sequences for which the term O(| θ -θ |2 ) is asymptotically negligible. As we know L ( θ) is 0, rearranging and multiply through by √ N , one obtains:

√ N ( θ -θ) = √ N L (θ)(-L (θ)) -1 . (3.3.26)
Therefore, the first factor in (3.3.26) can be written as

√ N L (θ) = √ N 1 N n i=1 T t=1 L (θ) -0 = √ N 1 N n i=1 T t=1 L (θ) -E[L (θ)] → N [0; I(θ)],
where the score or efficient score 2 is the gradient (the vector of partial derivatives) with respect to some parameter θ. Namely L (θ) = s(θ) is the score function with E[L (θ)] = 0 and the Fisher information for θ is I(θ) with blok-diagonal matrix as follows I(θ) = diag(I(β), I(γ), I(δ)) . By the law of large numbers -L (θ) → I(θ)

Then √ N ( θ -θ) = √ N L (θ)(-L (θ)) -1 → N 0; I(θ)[I(θ)I (θ)] -1 → N (0; I -1 (θ))

The Posterior distribution

The estimation formulas for parameters uses the conditional distribution for the state space given the observations (posterior distribution) p(X i | Y i ). This density can be computed by a Bayesian approach. In the linear-gaussian state-space model the posterior distribution certainly is gaussian. Consequently, the conditional mean and variance-covariance matrix of the state vector are computed by Kalman filter and smoother recursions. With the general non-gaussian state-space model, the distribution of the state space X i (t) is generally non-gaussian. Actually, the model in equations (3.3.1)and (3.3.3) is a non-gaussian state space model. Therefore, it is necessary to use the particles methods to find the approximation to the posterior distribution.

Here, we compute this distribution by the auxiliary iterated extended Kalman particle filter (AIEKPF) method, proposed by [START_REF] Xi | The auxiliary iterated extended Kalman particle filter[END_REF]. We develop this algorithm with our model where we find the derivation of the equations of posterior mode and posterior covariance to the extended Kalman filtering method, presented in Appendix A.2.

The Auxiliary Iterated Extended Kalman Particle Filter (AIEKPF)

As we mentioned in chapter 2, the main idea of the (AIEKPF) algorithm is that the importance density function is generated by the Iterated Extended Kalman Particle Filter (IEKF) method within an Auxiliary Particle filter (APF) method. The support of state variable (X i (t)) for each individual of state and observation equa-

tions (3.3.1, 3.3.3) is given as particle form {X m i (t), m = 1, 2, • • • , N } with associated weights {w m i (t), m = 1, 2, • • • , N } and X i (t) = (X i (0), X i (1), • • • , X i (T )) is the set of the state variables for one individual at time t = 1, 2, • • • , T. Recall the observation to an individual Y i (t) = (Y i (1), • • • , Y i (T )) where i = 1, • • • , n, t = 1, • • • , T.
In chapter 2, we noted that one can approximate the posterior distribution by a set of N particles as follows

p(X t | Y t-1 ) ≈ N m=1 w (m) t|t-1 δ(X t -X (m) t ).
The posterior mean is approximated as

X t|t-1 ≈ N m=1 w (m) t|t-1 X (m) t .
Then, the posterior probability density function can be approximated for an individual as

p(X i (t) | Y i (t)) ≈ N m=1 w m i (t)δ(X i (t) -X m i (t)), (3.4.1)
where δ(.) refers to the Dirac delta function, and the weights are normalized as

N m=1 w m i (t) = 1.
The weight w m i (t) of an individual is selected using the formula of importance sampling, which is written as

w m i (t) ∝ p[X m i (t) | Y i (t)] q[X m i (t) | Y i (t)] , (3.4.2)
where q[X m i (t) | Y i (t)] defines the importance density. Actually, APF methodology depends on the idea that introduces an auxiliary variable, ς i , i = 1, • • • , n which plays an important role of index of the mixture component, the augmented joint distribution p(X i (t), ς i | Y i (t)) with this extra auxiliary variable is updated for an individual i as follows :

p(X i (t), ς i = m | Y i (t)) ∝ p(Y i (t) | X i (t))p(X i (t), ς i = m | Y i (t)) = p(Y i (t) | X i (t))p(X i (t) | ς i = m, Y i (t)) ×p(m | Y i (t)). = p(Y i (t) | X i (t))p(X i (t) | X m i (t -1))w m i (t -1) = q k=1 p(Y ik (t) | X i (t)) × p(X i (t) | X m i (t -1))w m i (t -1) = q k=1 M(π ik (t)) × N (E[X m i (t)], Var[X m i (t)]) ×w m i (t -1), (3.4.3)
Previously, we assume that

Y i (t) = (Y i1 (t), • • • , Y iq (t)) , and
p(Y i (t) | X i (t)) = q k=1 p(Y ik (t) | X i (t)),
where the symbol M refers to multinomial distribution and symbol N refers to Normal distribution with mean and variance given respectively by

E[X m i (t)] = F (x m i (t -1), u m i (t), γ m ), Var[X m i (t)] = H 2 (x m i (t -1), u m i (t), δ m )R t . (3.4.4)
The APF generates a sample {X m i (t), ς m i } N m=1 , from the joint distribution p(X i (t), ς i | Y i (t)), where ς m i (ς m i ≡ {ς i = m}) refers to the index of particle m's parent to an individual. The importance density used in APF has the property that

q(X i (t), ς i | Y i (t)) ∝ p(Y i (t) | µ m i (t))p(X i (t) | X m i (t -1))w m i (t -1) ∝ q k=1 p(Y ik (t) | X i (t)) × p(X i (t) | X m i (t -1))w m i (t -1) ∝ q k=1 M [π ik (u m i (t), µ m i (t))] × N (E[X m i (t)]; Var[X m i (t)]) × w m i (t -1) (3.4.5)
and µ m i (t) is some description of X i (t) given x m i (t -1) like the mean, mode or another statistics. If it is the mean then

µ m i (t) = E[X i (t) | X m i (t -1)] or a sample µ m i (t) ∼ P (X i (t) | X m i (t -1)
). Then, the auxiliary variable is omitted, the sample {X m i (t)} N m=1 is obtained. The importance density can be rewritten as

q(X i (t), ς i | Y i (t)) = q(ς i | Y i (t))q(X i (t) | ς i , Y i (t)). (3.4.6) Then q(X i (t) | ς m i , Y i (t)) p(X i (t) | X m i (t -1)) = N (E[X m i (t)]; Var[X m i (t)]), (3.4.7)
where E[X m i (t)] and Var[X m i (t)] are as in equations (3.4.4). By substituting Eqs. (3.4.6) and (3.4.7) into Eq. (3.4.5), one obtains for an individual i :

q(ς m i | Y i (t)) p(Y i (t) | µ m i (t))w m i (t -1) = q k=1 M π ik (u m i (t), µ m i (t)) w m i (t -1). (3.4.8)
The importance weights of an individual i are recursively updated as

w m i (t) ∝ p[X m i (t), ς m i | Y i (t)] q[X m i (t), ς m i | Y i (t)] . ( 3
w m i (t) ∝ p(Y i (t) | X m i (t))p(X i (t) | X m i (t -1))w m i (t -1) p(Y i (t) | µ ς m i i (t))p(X i (t) | X m i (t -1))w m i (t -1) = p[Y i (t) | X m i (t)] p[Y i (t) | µ ς m i i (t)] = q k=1 p[Y ik (t) | X m i (t)] q k=1 p[Y ik (t) | µ ς m i i (t)] = q k=1 M [π ik (X m i (t), u m i (t))] M π ik (µ ς m i i (t), u m i (t))
.

(3.4.10)

The IEKF recursions computes the following recursive approximation to the true posterior filtering density of an individual

p(X i (t) | Y i (t)) = N X ij (t); P ij (t) (3.4.11) where i = 1, • • • , n , t = 0, • • • , T, j = 1, • • • , c
and c is the number of the Iterated Extended Kalman Filter (IEKF), X ij (t) is the iterative value of the X i (t) on the jth iteration, and P ij (t) is the covariance matrix of X ij (t).

According to the IEKF update equation , X ij (t) and P ij (t) can be updated for an individual as follows :

X ij (t | t) = X ij (t | t -1) + K ij (t)[ Y i (t) -π i (t)] (3.4.12) 1 ≤ i ≤ n, 1 ≤ k ≤ q, 1 ≤ t ≤ T, 1 ≤ s ≤ c k , where Y i (t) = D -1 i (t){Y i (t) -π i (t)} + η i (t)
and

π i (t) = ( π i1 (t), • • • , π ik (t), • • • π iq (t)) , π ik (t) = ( π 1 ik (t), • • • , π ck ik (t)) since π s ik (t) = exp[u i (t)β s k + X i (t | t -1)] 1 + c k r=1 exp[u i (t)β s k + X i (t | t -1)] , s = 1, • • • , c k . K ij (t) = P ij (t | t -1)B ij (t)[B ij (t)P ij (t | t -1)B ij (t) + Σ -1 it ] -1 P ij (t | t) = (I -K ij (t)B ij (t)P ij (t | t -1), (3.4.13)
where K ij (t) the Kalman gain matrix for individual i at instant t in the jth iteration.

The importance density for each particle to one individual i is generated by using IEKF method within the APF method as follows

q(X m i (t), ς i | Y i (t)) = N X ς m i ij (t); P ς m i ij (t) , (3.4.14)
where m = 1, 2, • • • N, and N is the number of the particles.

In brief, AIEKPF method uses IEKF method to update the equations with the new observations to compute the mean and covariance of the importance density for each particle at time t -1 to one individual. Subsequently, the mth particle is sampled from the distribution.

Auxiliary Iterated Extended Kalman Filter (AIEKPF) Algorithm

The AIEKPF algorithm can be performed for an individual by the following steps:

1. Initialization (t = 0) : For m = 1, • • • , N , generate the states (particles) X m i (0) from the prior p(X i (0)) = N ( X m i (0); P m i (0)), and set

X m i (0) = E[X m i (0)] = F (x m i (0), u m i (0), δ m )
and

P m i (0) = E[(X m i (0) -X m i (0))(X m i (0) -X m i (0)) ] = Var(X m i (0)) = H 2 (x m i (0), u m i (0), δ m )R 0 For t = 1, • • • , T, repeat the following steps: 2. For m = 1, • • • , N , generate µ m i (t) ∼ N (E[X m i (t)]; Var[X m i (t)]).
3. For m = 1, • • • , N , using the IEKF algorithm to update the particles 3-1. Calculate the Jacobians A m i (t), C m i (t) of the process model

A m i (t) = ∂F (x, u i (t), γ) ∂x | x=x m i (t-1|t-1) C m i (t) = H(x m i (t -1 | t -1), u i (t), δ) 3-2.
Predict the particle with the IEKF:

X m i (t | t -1) ≈ F (x m i (t -1 | t -1), u i (t), γ) P m i (t | t -1) = A m i (t)P m i (t -1 | t -1)A m i (t) + C m i (t)R t C m i (t) 3-3. For j = 1, • • • c (c is the iteration number of the IEKF) a-Calculate the Jacobians , B m ij (t) B m ij (t) = ∂π it (u i (t), x) ∂x | x=x m ij (t|t-1)
b-The state estimation error covariance P m ij (t) is updated with equation 3.4.13 c-The state estimate is updated X m ij (t) with equation 3.4.12.

For

m = 1, • • • N, compute w m i (t) = q(m | Y i (t)) ∝ q k=1 M [π ik (u m i (t), µ m i (t))] w m i (t -1), 5
. Apply the Resample method to obtain the index ς m i of particle m's parent. (3.4.15) where j = c, where the j = 1, • • • , c the last iteration of IEKF before this step is equal to c.

Importance sampling : for

m = 1, • • • , N , 6-1. Draw samples X m i (t) ∼ q(X i (t), ς m i | Y i (t)) = N ( X ς m i ij (t), P ς m i ij (t)),

6-2.

Compute the importance weights of particles by using

w m i (t) = p[Y i (t) | X m i (t)] p[Y i (t) | µ ς m i i (t)] = q k=1 p[Y ik (t) | X m i (t)] q k=1 p[Y ik (t) | µ ς m i i (t)] = q k=1 M [π ik (X m i (t), u m i (t))] M π ik (µ ς m i i (t), u m i (t))
.

(3.4.16) 6-3. Normalize the weights

w m i (t) = w m i (t) N m=1 w m i (t)
.

Output: a set of weighted particles (samples) to an individual [{X

m i (t), w m i (t)} N m=1 ], i = 1, • • • , n.
m i (t), w m i (t)} N m=1 ] as N → ∞, for all i = 1, • • • , n and t = 1, • • • , T, √ N ( ϕ i (t) -φi (t)) ⇒ N [0; σ 2 AP F (ϕ i (t))], (3.4.17)
where ϕ i (t) is the function of trajectories X i (t) and φi (t) is the expectation of ϕ i (t) with respect to the filtering distribution.

φi (t) = ϕ i (t)p(X i (t) | Y i (t))dX i (t)
and

ϕ i (t) = N m=1
W m i (t)ϕ m i (t), (3.4.18) where

W m i (t) = w m i (t)[ N m=1 w m i (t)]
w m i (t) = p[Y i (t) | X m i (t)] p[Y i (t) | µ ς m i i (t)] = q k=1 p[Y ik (t) | X m i (t)] q k=1 p[Y ik (t) | µ ς m i i (t)] = q k=1 M [π ik (X m i (t), u m i (t))] M π ik (µ ς m i i (t), u m i (t))
.

At time t = 1 for an individual, the formula of the variance yields

σ 2 AP F (ϕ i (1)) = (p(X i (1) | Y i (1)) 2 q 1 (X i (1)) (ϕ i (1) -φi (1)) 2 dX i (1), (3.4.20)
and for t > 1,

σ 2 AP F (ϕ i (t)) = [p(X i (1) | Y i (t))] 2 q 1 (X i (1)) ϕ (1,t) i (X i (1)) 2 dX i (1) + t-1 r=2 [p(X i (1 : r) | Y i (t))] 2 p(X i (r -1) | Y i (r))q r (X i (r) | X i (r -1)) × ϕ (r,t) i (X i (r)) 2 dX i (r) + [p(X i (t) | Y i (t))] 2 p(X i (t -1) | Y i (t))q t (X i (t) | X i (t -1)) ×(ϕ i (t) -φi (t)) 2 dX i (t). (3.4.21)
where

ϕ (r,t) i (X i (r)) = ϕ i (X i (t))p(X i (r + 1 : t) | Y i (r + 1 : t))dX i (r + 1 : t) -φi (t)
and

q 1 (X i (1)) = N (E(X m i (1)), Var(X m i (1))) q t (X i (t) | X i (t -1)) = N (E(X m i (t)), Var(X m i (t))) p(X i (t) | Y i (t)) ≈ N m=1 w m i (t)δ(X i (t) -X m i (t)).
where δ(.) denotes the Dirac delta function.

Derivation of the posterior density function

We calculated the information matrix via the EM algorithm, in the second component the derivative of posterior density with respect to θ (p) is used, we recall θ (p) is the vector of parameters at iterative p. By using the approach proposed by [START_REF] Poyiadjis | Maximum likelihood parameter estimation in general state-space models using particle methods[END_REF]. This approach approximates the filter derivative based on the sequence of marginal distributions p(X i (t)

| Y i (t)). Let p(X i (t) | Y i (t), θ (p) ) = ξ[X i (t), Y i (t), θ (p) ] ξ[X i (t), Y i (t), θ (p) ]dX i (t) , (3.4.22)
where

ξ(X i (t), Y i (t), θ (p) ) = p(Y i (t) | X i (t), θ (p) )p(X i (t) | Y i (t -1), θ (p) ) = q k=1 M[π ik (t)] × p(X i (t) | Y i (t -1), θ (p) ),
with M denoting to the multinomial distribution. Then the gradient of (3.4.22 ) is given by

∇p(X i (t) | Y i (t), θ (p) ) = ∇ξ[X i (t), Y i (t), θ (p) ] ξ[X i (t), Y i (t), θ (p) ]dX i (t) -p(X i (t) | Y i (t)) ∇ξ[X i (t), Y i (t), θ (p) ]dX i (t) ξ[X i (t), Y i (t), θ (p) ]dX i (t) , (3.4.23)
where ∇ is the first-order derivatives of the density with respect to θ (p) , and

∇ξ(X i (t), Y i (t), θ (p) ] = p(Y i (t) | X i (t)) p(X i (t) | X i (t -1))p(X i (t -1) | Y i (t -1) ×[∇ log p(Y i (t) | X i (t)) + ∇ log p(X i (t) | X i (t -1))]dX i (t -1) + p(Y i (t) | X i (t)) p(X i (t) | X i (t -1)) ×∇p(X i (t -1) | Y i (t -1))dX i (t -1).
The filter derivative approximations by the weighted particles can be obtained as follows:

∇p(X i (t) | Y i (t), θ (p) ) = N m=1 w m i (t)α m i (t)δ(X i (t) -X m i (t)), (3.4.24)
where δ(.) denotes the Dirac delta function, and for all i = 1,

• • • , n, t ≥ 1,
w m i (t) = w m i (t) N m=1 w m i (t)
, with

w m i (t) = q k=1 M [π ik (X m i (t), u m i (t))] M π ik (µ ς m i i (t), u m i (t)) . (3.4.25)
Recall that w m i (t) is the weighted particles calculated by the Auxiliary Particle Filter (APF) algorithm or Auxiliary Iterated Extended Kalman Particle Filter(AIEKPF) Algorithm, and

w m i (t)α m i (t) = ρ m i (t) N m=1 w m i (t) -w m i (t) N m=1 ρ m i (t) N m=1 w m i (t) , (3.4.26)
where α m i (t) corresponds to an approximation of the so-called score ∇p(X m i (t) | Y i (t), θ (p) ) :

α m i (t) = [w m i (t)] -1 ∇p(X m i (t) | Y i (t), θ (p) ) q[X m i (t) | Y i (t), θ (p) ] = [w m i (t)] -1 q k=1 ∇M [π ik (X m i (t), u m i (t))] M π ik (µ ς m i i (t), u m i (t)) , (3.4.27)
and

ρ m i (t) = ∇ξ[X m i (t), Y i (t), θ (p) ] p[X m i (t) | Y i (t)] (3.4.28) since ∇ξ[X m i (t), Y i (t), θ (p) ] = N m=1 w m i (t)p(Y i (t) | X i (t))p(X i (t) | X m i (t -1)) [∇ log p(Y i (t) | X i (t)) + ∇ log p(X i (t) | X m i (t -1)) + α m i (t -1)] = N m=1 w m i [ q k=1 M(π ik (t)] × N (µ m i (t), V m i (t)) ∂l(β) ∂β s k + ∂l(γ) ∂γ + ∂l(δ) ∂δ + α m i (t -1)
where

∂l(β) ∂β s k = n i=1 T t=1 [y s ik (t) -π s ik (t)]u i (t)D ik (X t )Σ -1 ik (X t ) ∂l(γ) ∂γ = n i=1 T t=1 [X i (t) -µ i (t)] V i (t) × ∂µ i (t) ∂γ ∂l(δ) ∂δ = n i=1 T t=1 1 V i (t) - [X i (t) -µ i (t)] 2 V 2 i (t) × ∂V i (t) ∂δ .

Posterior mode estimation

In nonlinear non gaussian state space model, there are two approaches to estimate the state variable. First, the Bayesian approach which uses A maximum a posteriori (MAP). Regrettably, this approach leads to numerical integration problems, often difficult for more complicated models. Several methods are proposed to solving this problem. One of these Recent Markov Chain Monte Carlo methods solving this problem by repeated sampling from approximative posterior distributions. But there are still open questions about achieving the convergence of sampling schemes. Second, penalized likelihood estimation approach was proposed by Fahrmeir and Kaufmann (1991), [START_REF] Fahrmeir | Posterior mode estimation by extended Kalman filtering for multivariate dynamic generalized linear models[END_REF] and [START_REF] Fahrmeir | Penalized likelihood estimation and iterative Kalman smoothing for non-Gaussian dynamic regression models[END_REF]. This approach can also be interpreted as nonparametric method for the state space models. They estimated the posterior mode by using Fisher scoring method by iterative Kalman filtering and smoothing. They used the numerical method to maximize the penalized likelihood called "Working Kalman filtering and smoother (WKFS)" with the exponential family distribution.

In this chapter, we used the two approaches to find the estimation of state variable. We develop the two approaches with our model to find the posterior mode as follows :

1. Penalized likelihood estimation : a -Gauss-Newton and Fisher-scoring Filtering and smoothing algorithms.

b -Working extended Kalman filter and smoother algorithms.

Maximum a posteriori sequence estimation (MAP):

a -The auxiliary iterated extended Kalman filter particle filter-MAP.

Penalized likelihood estimation

The posterior mode smoother is defined as follows

a ≡ a (0 | T ), a (1 | T ), • • • , a (T | T ) ∈ R m ,
with m = (T + 1)n, where

a(0 | T ) = ((a 1 (0 | T ), • • • , a n (0 | T )) , . . . . . . . . . . . . a(T | T ) = (a 1 (T | T ), • • • , a n (T | T )) .
The posterior likelihood of X by Bayes' theorem

p(X | Y) = 1 p(Y) n i=1 q k=1 T t=1 p(Y ik (t) | X i (t)) × n i=1 T t=1 g i (X i (t)) n i=1 g i (X i (0)) (3.5.1) p(Y) does not depend on X p(X | Y) ∝ n i=1 q k=1 T t=1 p(Y ik (t) | X i (t)) × n i=1 T t=1 g i (X i (t)) n i=1 g i (X i (0)) (3.5.2)
where the prior likelihood g i (X i (0))

g i (X i (0)) = 1 2πV i (0) exp -[X i (0) -µ i (0)] 2 2V i (0) . (3.5.3)
Taking logarithms and inserting the densities of X in equation (3.3.5 ) and X i (0) in equation (3.5.3), the panelized log-likelihood function yields

P L : R m -→ R, m = (T + 1)n P L(X) = n i=1 q k=1 T t=1 {log p(Y ik (t) | X i (t))} (3.5.4) + n i=1 T t=1 logg i (X i (t)) + n i=1 log g i (X i (0)) where log p(Y ik (t) | X i (t)) = log c k s=1 [π s ik (t)] Y s ik (t) = c k s=1 Y s ik (t) log π s ik (t). Then P L(X) = n i=1 q k=1 T t=1 c k s=1 {Y s ik (t) log π s ik (t)} + G 1 + G 2 (3.5.5)
where

G 1 = - 1 2 n i=1 (X i (0) -µ i (0))V -1 i (0)(X i (0) -µ i (0)) G 2 = - 1 2 n i=1 T t=1 (X i (t) -µ i (t)) V -1 i (t)(X i (0) -µ i (t)),
The posterior mode smoother yields

a ≡ (a (0 | T ), a (1 | T ), • • • , a (T | T )) = arg max X {P L(X)}, (3.5.6)
The maximization of p(X | Y) is equivalent to maximization of the penalized loglikelihood (3.5.4). Numerical maximization of the penalized log-likelihood can be performed by various algorithms. In this chapter we use two algorithms to find the posterior mode. First, the iterative forward-backward Gauss-Newton (Fisher-scoring) algorithms. Second, the Working Extended Kalman Filter and Smoother (WEKFS) which gives the same numerical results computed by Gauss-Newton (Fisher-scoring) algorithms.

Gauss-Newton and Fisher-scoring Filtering and smoothing

A maximization of the penalized log-likelihood P L(X) with generally best performance to approximation quality can be computed by Gauss-Newton or Fisher-scoring iterations. In other words, this can be also achieved by applying extended Kalman filtering and smoothing to a "working model" in each Fisher-scoring iteration. Therefore, the penalized log-likelihood criterion (3.5.4) can be presented in compact matrix notation as:

P L(X) = l(X) - 1 2 X KX, (3.5.7)
where

l(X) = n i=1 q k=1 T t=0 c k s=1 {Y s ik (t) log(π s ik (t)} .
The penalty matrix for one individual K i is symmetric and block-tridiagonal, with blocks that can be easily computed from (3.5.4):

K i =         K i,00 K i,01 • • • • • • 0 K i,10 K i,11 K i,12 • • • . . . . . . K i,21 . . . . . . . . . . . . . . . . . . K i,T -1,T 0 • • • • • • K i,T,T -1 K i,T,T         with µ i (t) = F (X i (t-1), u i (t), γ)
, the conditional mean is rewritten by matrix notation as follows

F (X i (t -1), u i (t), γ) = Ψ it X i (t), then K i,t-1,t = K i,t,t-1 , i = 1, • • • , n, t = 1, • • • , T K i,00 = Ψ i,1 R -1 t Ψ i,1 , K i,tt = R -1 t + Ψ i,t+1 R -1 t Ψ i,t+1 , i = 1, • • • , n, t = 1, • • • , T Ψ i,T +1 = 0, K i,t-1,t = -Ψ i,t R -1 t , i = 1, • • • , n, t = 1, • • • , T,
We now describe the steps of Fisher scoring in vector notation. Recall that for all i = 1, • • • , n,

Y i = (Y i (0), Y i (1), • • • , Y i (T )) , Y i (t) = (Y i1 (t) • • • , Y iq (t)) , t = 0, 1, 2, • • • , T.
Fahrmeir and Wagenpfeil (1997) assumed Y i (0) = a i (0) and the vectors of expectations are defined by

Π i (X) = (π i0 (X 0 ), π i1 (X 1 ), • • • , π iT (X T ))
Fahrmeir and Wagenpfeil (1997) assumed π i (0) = X i (0). We recall the conditional mean and variance of one individual

E(Y ik (t) | X i (t)) = π s ik (t), s = 1, • • • , c k Var(Y ik (t) | X i (t)) = Σ ik (X t ), where Y ik (t) = (Y 1 ik (t), • • • , Y c k ik (t))
, and Σ ik (X t ) has generic elements

σ sm ik (t) = π s ik (t)[1 -π s ik (t)], if s = m -π s ik (t)π m ik (t) if s = m (3.5.8)
The diagonal covariance matrix of an individual i at time t is

Σ i (X) = diag(V i (0), Σ i1 (X 1 ), • • • , Σ iT (X T )),
where

V i (0) = H 2 [X i (0), u i (0)
, δ]R 0 , and the diagonal matrix

D i (X) = diag(1, D i1 (X 1 ), • • • , D iT (X T )).
Since D it (X t ) is first-order derivative of the conditional probability π i (t) evaluated at η i (t), the score function of l(X) in (3.5.7) is given by

S i (X) = ( S i0 (X 0 ), S i1 (X 1 ), • • • , S iT (X T )) ,
where for all i = 1,

• • • , n, S i (X) := D i (X)Σ -1 i (X) {Y i (t) -Π i (X)} , (3.5.9)
with components

S i (X 0 ) = V -1 i (0)(a i (0) -X i (0)) (3.5.10) S it (X t ) = D it (X t )Σ -1 it (X t ) {Y i (t) -π it (X t )} , t = 1, • • • , T, (3.5.11)
The first-order derivatives of P L(X) in (3.5.7) are

M (X) = ∂P L(X)/∂X = S(X) -KX, (3.5.12)
and the expected information matrix is given by

I i (X) = (I i0 (X 0 ), I i1 (X 1 ) • • • , I iT (X T ))
where for all i = 1, • • • , n

I i (X) = D i (X)Σ -1 i (X)D i (X) (3.5.13)
with diagonal blocks

I i0 (X 0 ) = V -1 i (0) (3.5.14) I it (X t ) = D it (X t )Σ -1 it (X t )D it (X t ), t = 1, • • • , T. (3.5.15)
Taylor expansion of the score function about X 0 yields

M ( X | Y) ≈ M (X 0 ) -I(X 0 ) × X 1 -X 0 since M ( X | Y) = 0, a single Fisher scoring to the next iterate X 1 ∈ R m , with m = (T + 1)n is as follows I(X 0 ) × X 1 -X 0 = M (X 0 ).
This can be rewritten as

X 1 = I(X 0 ) + K -1 I(X 0 ) Y, (3.5.16)
with "working" observation

Y = ( Y 1 , • • • , Y n ) and Y i = ( Y i (0), Y i (1), • • • , Y i (T )))
, then can compute as

Y i := D -1 i ((X)) {Y i -Π i (X))} + η i (X), (3.5.17)
with components

Y i (0) = a i (0) Y i (t) = D -1 it (X t ) {Y i (t) -π it (X t )} + η it (X t ), t = 1, • • • , T, where η i ((X)) = (η i1 (X 1 ), • • • , η iT (X T ))
, is the vector of link function for i individual.

Working Extended Kalman Filter and Smoother (WEKFS)

In the following algorithm, a t|t-1 , , a t|t , a t|T are numerical approximations to predicted, filtered, smoothed values of posterior mode to X and the corresponding P t|t-1 , P t|t , P t|T are numerical approximations to predicted, filtered, smoothed values of error covariance matrices. We expand the non-linear functions F and H in the state equation and the conditional probability π s ik (t) in Taylor series expansion about the estimation value of X. The derivation of Extended Kalman filter and smoother is presented in Appendix (A.2) Initialization:

a i (0 | 0) = a i (0), V i (0 | 0) = V i (0).
(3.5.18)

Prediction: for t = 1, • • • , T , a i (t | t -1) = F (a i (t -1 | t -1), u i (t), γ) P i (t | t -1) = A i (t)P i (t -1 | t -1)A i (t) + C i (t)R t C i (t), (3.5.19)
where

A i (t) = ∂F (x, u i (t), γ) ∂x | x=a i (t-1|t-1)
.

Filtering :

for t = 1, • • • , T , a i (t | t) = a i (t) + K i (t)( Y i (t) -a i (t | t -1)) K i (t) = P i (t -1 | t -1)B i (t)(B i (t)P i (t -1 | t -1)B i (t) + I -1 i (t)) -1 P i (t | t) = (I -K i (t)B i (t))P i (t -1 | t -1), (3.5.20)
where

B i (t) = ∂π it ∂x | x=a i (t|t-1) .

Iteratively extended Kalman filter and smoother (IWEKFS)

Initialization: Calculate X 0 = (a 0 0|t , a 0 1|t , • • • , a 0 T |t ) with General Kalman filter and smoother (GKFS) recursions. Set iteration index k = 0. STEP 1: Starting with X k , compute X k+1 by performing Working extended Kalman filter and smoother (WEKFS) recursions.

STEP 2:

If the convergence criterion is achieved : STOP, else set k = k + 1 and go to STEP 1.

The maximum a posteriori (MAP) with particle filtering and smoothing algorithm

The second approach to finding the posterior mode is the MAP. Firstly, the extension of the MAP via the particle filtering is proposed by [START_REF] Saha | On the Monte Carlo marginal MAP estimator for general state space models[END_REF] with our model. Secondly, we developed the MAP via the Auxiliary Iterated Extended Kalman Filter (AIEKPF) Algorithm. The filtering and smoothing step with our model is as follows:

Filtering step

As we mentioned in chapter 2, the particle based filter MAP estimator (pf-MAP) proposed by [START_REF] Saha | On the Monte Carlo marginal MAP estimator for general state space models[END_REF] has the advantage that the posterior density can be approximated not only by particles shape at clouds but at any point.

The MAP estimate of the filtering density at time t to an individual is defined as follows

X M AP i (t) = arg max X i (t) p(X i (t) | Y i (t)). (3.5.21)
By using the Baye's rule , the filtering density yields

p(X i (t) | Y i (t)) = p(X i (t) | Y i (t), Y i (t -1)) = p(Y i (t) | X i (t), Y i (t -1))p(X i (t) | Y i (t -1)) p(Y i (t) | Y i (t -1)) = p(Y i (t) | X i (t))p(X i (t) | Y i (t -1)) p(Y i (t) | Y i (t -1)) (3.5.22)
Therefore, the MAP estimate is as follows

X M AP i (t) = arg max X i (t) p(Y i (t) | X i (t))p(X i (t) | Y i (t -1)) p(Y i (t) | Y i (t -1)) . (3.5.23) Since the denominator p(Y i (t) | Y i (t -1)) is independent of X i (t)
, the MAP estimate can be expressed as

X M AP i (t) = arg max X i (t) [p(Y i (t) | X i (t))p(X i (t) | Y i (t -1))] = arg max X i (t) q k=1 p(Y ik (t) | X i (t))p(X i (t) | Y i (t -1)
= arg max

X i (t) q k=1 M(π ik (t))p(X i (t) | Y i (t -1) . (3.5.24)
The main principle for evaluating MAP is the computation of the predictive density

P (X i (t) | Y i (t - 1
)) (the second terms in equation (3.5.24) ) which is not available in known formula. The prediction density formula in Bayesian approach is as follows

p(X i (t) | Y i (t -1)) = p(X i (t), X i (t -1) | Y i (t -1))dX t-1 = p(X i (t) | X i (t -1), Y i (t -1))p(X i (t -1) | Y i (t -1))dX i (t -1) = p(X i (t) | X i (t -1)) × p(X i (t -1) | Y i (t -1))dX i (t -1) = N (E[X m i (t)]; Var[X m i (t)]) × p(X i (t -1) | Y i (t -1))dX i (t -1). (3.5.25)
The posterior distribution can be approximated in particle form by a set of N weighted particles as:

p(X i (t) | Y i (T )) ≈ N m=1 w m i (t)δ(X i (t) -x m i (t)) (3.5.26)
Now, use equation (3.5.26)

to approximate p(X i (t) | Y i (t -1)) as p(X i (t) | Y i (t -1)) ≈ N j p(X i (t) | X j i (t -1))w (j) i (t -1). ≈ N j N (µ j i (t); V j i (t))w (j) i (t -1). (3.5.27)
Substituting (3.5.27) into (3.5.24), the MAP estimation is obtained by finding the global maxima of the posterior distribution which is approximated by particle as

X M AP i (t) = arg max X (m) i (t) p(Y i (t) | X (m) i (t)) × j p(X (m) i (t) | X (j) i (t -1))w (j) i (t -1) = arg max X (m) i (t) q k=1 M(π ik (t)) × j N (µ j i (t), V j i (t))w (j) i (t -1)
,

where i = 1, • • • , n, t = 1, • • • , T, m = 1, • • • , N.
X M AP i (t | T ) = arg max X i (t) p(X i (t) | Y i (T )) = arg max X (m) i (t) p(X (m) i (t) | Y i (t)) w (m) i (t | T ) w (m) i (t) . (3.5.28)
The filtered density p(X i (t) | Y i (t)) at the particle cloud {X

} N m=1 can be approximated during the forward filtering step as

p(X i (t) | Y i (t)) ≈ p(Y i (t) | X (m) i (t)) j p(X (m) i (t) | X (j) i (t -1))w (j) i (t -1) p(Y i (t) | Y i (t -1)) (3.5.29) Since p(Y i (t) | Y i (t -1)) in equation (3.5.29) is independent of X (m) i (t), X M AP i (t | T ) can be obtained by replacing p(X (m) i (t) | Y i (t)) by the filtered density q(X (m) i (t) | Y i (t)) = p(Y i (t) | X (m) i (t)) j p(X (m) i (t) | X (j) i (t -1))w (j) i (t -1). = q k=1 M(π m ik (t)) j N (µ j i (t), V j i (t))w (j) i (t -1) (3.5.30)
The derivation of the Forward-Backward smoothing is given in Appendix (A.4).

Auxiliary Iterated Extended Kalman Particle Filter -MAP(AIEKPF-MAP) Algorithm

the MAP estimation can be found via the Auxiliary Iterated Kalman Particle Filter (AIEKPF) with the weights to individual i at time t as follows :

w m i (t) = p[Y i (t) | X m i (t)] p[Y i (t) | µ ς m i i (t)] = q k=1 M[π ik (X m i (t), u i (t))] M[π ik (µ ς m i i (t), u i (t))]
The AIEKPF algorithm can be performed with the following steps:

1. Initialization (t = 0) : For m = 1, • • • , N , draw the states (particles) X m i (0) from the prior p(X i (0)) = N ( X m i (0), P i (0) m ), and set

X m i (0) = E[X m i (0)] = F (X m i (0), u m i (0), δ m )
and

P m i (0) = E[(X m i (0) -X m i (0))(X m i (0) -X m i (0)) ] = Var(X m i (0)) = H 2 (X m i (0), u m i (0), δ m )R 0 For t = 1, • • • , T, repeat the following steps: 2. For m = 1, • • • , N , generate µ m i (t) ∼ N (E(X m i (t)); Var(X m i (t))).
3. For m = 1, • • • , N , update the particles via the IEKF algorithm 3-1. Compute the Jacobians A m i (t), C m i (t) of the model as follows

A m i (t) = ∂F (x, u i (t), γ) ∂x | x=X m i (t-1|t-1) C m i (t) = H(X m i (t -1 | t -1), u i (t), δ) 3-2.
Predict the particle with the IEKF:

X m i (t | t -1) ≈ F (X m i (t -1 | t -1), u i (t), γ) P m i (t | t -1) = A m i (t)P m i (t -1 | t -1)A m i (t) + C m i (t)R t C m i (t) 3-3 For j = 1, • • • c (c is the iteration number of the IEKF) a-Compute the Jacobians , B m ij (t) B m ij (t) = ∂π it (u i (t), x) ∂x | x=x m ij (t|t-1)
b-Update the state estimation error covariance P m ij (t) with Eq.( 3.4.13) c-Update the state estimate X m ij (t) with Eq.(3.4.12).

For

m = 1, • • • N, calculate w m i (t) = q(m | Y i (t)) ∝ q k=1 M(π ik (u i (t), µ m i (t))w m i (t -1), 5.
Resample to obtain the index ς m i of particle m's parent.

Importance sampling : for

m = 1, • • • , N , 6-1. Draw samples X m i (t) ∼ q(X i (t), ς m i | Y i (t)) = N ( X ς m i ij (t); P ς m i ij (t))
, where j = c 6-2. Calculate importance weights of particles by using

w m i (t) = p[Y i (t) | X m i (t)] p[Y i (t) | µ ς m i i (t)] = q k=1 M[π ik (x m i (t), u i (t))] M[π ik (µ ς m i t (t), u i (t))] 6-3. Normalize the weights w m i (t) = w m i (t) N m=1 w m i (t)
.

Forward Filtering step

X M AP i (t) = arg max X (m) i (t) p(Y i (t)(t) | X (m) i (t)) j p(X (j) i (t) | X (j) i (t -1))w (j) i (t -1) = arg max X (m) i (t) q k=1 M(π m ik (t)) j N (µ j i (t); V j i (t)).w (j) i (t -1) .
Backward smoothing step

• Set w (m) i (T | T ) = w (m) i (t) • For t = T -1, • • • , 1 compute the smoother importance weights as w (m) i (t | T ) = w (m) i (t) N j=1 w (j) i (t + 1 | T ) p(X (j) i (t + 1) | X (m) i (t)) N r=1 p(X (j) i (t + 1) | X (r) i (t))w (r) i (t) = w (m)(t) i N j=1 w (j) i (t + 1 | T ) N (µ j i (t); V j i (t)) N r=1 N (µ r i (t); V r i (t))w (r) i (t)
.

• Compute the approximate smoother MAP as

X M AP i (t | T ) = arg max X (m) i (t) q(X (m) i (t) | Y i (t)) w (m) i (t | T ) w (m) i (t) . = arg max X (m) i (t) q k=1 M(π m ik (t)) j N (µ j i (t); V j i (t))w (j) i (t -1) w (m) i (t | T ) w (m) i (t)
.

Chapter 4

Longitudinal multicategorical processes : Generalized state space models with exponential family state noise

Introduction

In this chapter, we present longitudinal multivariate processes which are multicategorical

Y ik (t) = (Y 1 ik (t), • • • , Y c k ik (t)) , i = 1, • • • , n, t = 1, • • • , T and k = 1, • • • , q, .
As in chapter 3, the latent processes X i (t) are described by a CHARN model but with the state noise from exponential family distribution. Practically, many of the most common distributions belong to the exponential families, as Normal, Exponential, Gamma, Chi-squared, Beta, Dirichlet, Bernoulli, Multinomial, Poisson, Wishart, Inverse Wishart and many others distributions. In this chapter, we generalize our model with the state noises assumed to be an exponential family distribution. In the state space models if the observations and the latent variable have gaussian distribution, the posterior distribution is also a gaussian distribution. Therefore, the classical Kalman filtering can be used to estimate the posterior mean and variance matrix [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF], [START_REF] Kitagawa | Introduction to time series modeling[END_REF]). Here, the observations are from a multinomial distribution and the latent variables are from a CHARN model with noise from an exponential family distribution. Consequently, the posterior distribution is non symmetric. That is why we estimate the latent variables by the posterior mode instead of the posterior mean. As in chapter 3, we find the posterior mode by two approaches. First, via the Working Extended Kalman Filtering recursions (WEKF). Second, by using the Auxiliary Iterated Extended Kalman Particle Filtering (AIEKPF-MAP). The models parameters are estimated through the Maximum Likelihood (MLE) method via the Expectation-Maximization (EM) algorithm. Their consistency and asymptotic normality are estab-lished. The posterior distribution p(X i | Y i ) is approximated by the Auxiliary Iterated Extended Kalman Particle filter (AIEKPF).

The state space models

In this section we define as in the preceding chapter the observation and state equations. The observation equation does not change. In other words, the observations have multinomial distribution. The state equation CHARN model but the noise process has an exponential family distribution.

The observation equation

1. We assume that the conditional probability of

Y ik (t) given X i (t) is a multinomial distribution. for all i = 1, • • • , n, k = 1, • • • , q, t = 1, • • • , T, the conditional probability can clearly written as follows Pr[Y ik (t) = (y 1 ik (t), • • • y c k ik (t)) | X i (t) = x i (t)] = c k s=1 [π s ik (t)] y s ik (t) , (4.2.1)
where

π s ik (t) = exp[η s ik (t)] 1 + c k j=1 exp[η j ik (t)] if s < c k , π c k ik (t) = 1 1 + c k j=1 exp[η j ik (t)] , (4.2.2)
and c k s=1 y s ik (t) = 1, and c k s=1 π s ik (t) = 1. We recall that the link function η s ik (t) is defined with the logit function as follows

η s ik (t) = logit(π s ik (t)) = log π s ik (t) π c k ik (t) = log π s ik (t) 1 -c k -1 j=1 π j ik (t) = u i (t)β s k + X i (t), and
η c k ik (t) = log π c k ik (t) π c k ik (t) = log(1) = 0. u i (t) = (u i1 (t), • • • , u ir (t)
) is the vector of independent covariate variables and r is their number. For k = 1, • • • , q, the β s k = (β s k1 , • • • , β s kr ) are vectors of unknown regression parameters.

The vectors Y

ik (1), • • • , Y ik (T ) are conditionally independent given the vector of latent variable X i (1) = x i (1), • • • X i (T ) = x i (T ) : Pr[Y ik (1) = y ik (1), • • • , Y ik (T ) = y ik (T ) | X i (1) = x i (1), • • • , X i (T ) = x i (T )] = T t=1 Pr[Y ik (t) = y ik (t) | X i (t) = x i (t)] 3. The vectors Y i1 (t), • • • , Y ik (t) are conditionally independent given the latent vari- able X i (t) : Pr[Y i1 (t) = y i1 (t), • • • , Y ik (t) = y ik (t) | X i (t) = x i (t)] = q k=1 Pr[Y ik (t) = y ik (t) | X i (t) = x i (t)]

The state equation

The state equation is the following

X i (t) = F [X i (t -1), u i (t), γ] + H[X i (t -1), u i (t), δ]ε i (t), (4.2.3)
where

ε i (t) ∼ expf (ν i (t), φ i (t))
, that is the density functions of the ε i (t) s are :

f ε i (t) (z) = exp zν i (t) -b[ν i (t)] φ i (t) + c[z, φ i (t)] . (4.2.4)
with γ, δ the model parameters.

It results that the conditional distribution of x i (t) given x i (t -1) = x has density function :

f i X i (t-1)=x (z) = 1 H[x, u i (t), δ] f z z -F [x, u i (t), γ] H[x, u i (t), δ] (4.2.5)
• (ε i (t)) is the noise process for the state process.

• F (., ., .) : R × R r × R l -→ R is a non-linear function.

• H(., ., .) : R × R r × R l -→ R is a non-linear function.

• ν i (t) is a canonical parameter.

• φ i (t) is the dispersion or a scale parameter.

• b[ν i (t)], c[z, φ i (t)] are functions taking different forms depending on the distribution of the ε i (t) s.

The common distributions from exponential family with their functions and parameters are presented in Appendix (C.1)).

The latent process (X i (t) : 1 ≤ t ≤ T ) satisfies:

p(X i (t) | X i (t -1), X i (t -2), • • • , X i (1)) = p(X i (t) | X i (t -1)).
Recall that

µ i (t) = E[X i (t) | X i (t -1), u i (t)] = F [X i (t -1), u i (t), γ] V i (t) = Var[X i (t) | X i (t -1), u i (t)] = H 2 [X i (t -1), u i (t), δ]R t
Then the joint law of variables

X i = (X i (0), X i (1), • • • , X i (T )
) is deduced easily by conditioning

g i (X i ) = T t=1 p(X i (t) | X i (t -1)) × p(X i (0)) = 1 T t=0 H[X i (t -1), u i (t), δ] exp T t=0 Z i (t)ν i (t) -b[ν i (t)] φ i (t) + c[Z i (t), φ i (t)] . (4.2.6)
where

Z i (t) = X i (t) -F [X i (t -1), u i (t), γ] H[X i (t -1), u i (t), δ] (4.2.7)

The marginal likelihood

In this section, as chapter 3, we recall that for all i = 1, • • • , n,

Y i = (Y i (0), Y i (1), • • • , Y i (T )) , Y i (t) = (Y i1 (t) • • • , Y iq (t)) , t = 0, 1, 2, • • • , T. with Y ik (t) = (Y (1)
ik (t), • • • , Y (c k ) ik (t)) , t = 0, 1, 2, • • • , T, k = 1, • • • , q, and
X i = (X i (0), X i (1), • • • , X i (T )) , d(X i ) = (d(X i (0)), d(X i (1)), • • • , d(X i (T ))) .
The parameters of the model are θ = (β, γ, δ). β q-dimensional vectors β = (β 1 , • • • , β q ) , with β k s, are c k × r matrices, where c k is categories of item k, k = 1, • • • , q. r denotes the number of covariates.

We recall y i = (y i (0), y i (1), • • • , y i (T )) with y i (t) = (y i1 (t) • • • , y iq (t)) , the density function can be written as follows

Pr(Y 1 = y 1 , Y 2 = y 2 , • • • , Y n = y n ) = n i=1 • • • p(y i | X i ; θ)g i (X i )d(X i ) = n i=1 • • • T t=0 q k=1 Pr(Y ik (t) = y 1 ik (t), • • • , y c k ik (t) | X i (t))g i (X i )d(X i ) = n i=1 • • • T t=0 q k=1 c k s=1 [π s ik (t)] y s ik (t) × g i (X i )d(X i ). = n i=1 • • • T t=0 q k=1 c k s=1 exp[η s ik (t)] 1 + c k j=1 exp[η j ik (t)] y s ik (t) × g i (X i )d(X i ). (4.2.8)
where g i (X i ) is the joint density function of the latent variables X i defined by equation (4.2.6). Thus, the likelihood is

p(Y 1 , Y 2 , • • • , Y n ) = n i=1 • • • T t=0 q k=1 c k s=1 exp[u i (t)β s k + X i (t))] 1 + c k j=1 exp[u i (t)β j k + X i (t)] Y s ik (t) ×g i (X i )d(X i ), (4.2.9)

The EM algorithm

As in chapter 3 the Expectation and Maximizing step are:

Expectation -step: compute the expectation Q(θ | θ (p) ) given by

Q(θ | θ (p) ) = E{log[f (Y, X; θ)] | Y, θ (p) } = n i=1 • • • [log{g i (X i , θ i ) + log{p(Y i | X i )}] ×p(X i | Y i , θ (p) )d(X i )
The E-step for first-order Markov latent process from Exponential family distribution is :

Q(θ | θ (p) ) = G 1 + G 2 + G3 (4.2.10)
where

G 1 = n i=1 • • • T t=0 Z i (t)ν i (t) -b[ν i (t)] φ i (t) + c[Z i (t), φ i (t)] ×p(X i | Y i , θ (p) )d(X i ). (4.2.11) G 2 = n i=1 q k=1 T t=0 c k s=1 y s ik (t) • • • log[π s ik (t)] ×p(X i | Y i , θ (p) )d(X i ), (4.2.12) G3 = -log(H[X i (t -1), u i (t), δ]) (4.2.13) and p(X i | Y i , θ (p)
) is the conditional density of latent vector X i given the observation vector Y i (we further use the particle filtering algorithm to find it ).

Maximizing step:

θ (p+1) = arg max Q(θ | θ (p) )

Estimation of first-order CHARN latent processes

We recall that the parameter vector is θ

= (β 1 , • • • , β q , γ , δ ) , k = 1, • • • , q, and β k = (β 1 k , • • • , β c k k )
. We apply the E-M step for finding the MLE as follows maximize with respect to the

β s k , k = 1, • • • , q, s = 1, • • • , c k , only the part G 2 , β s(p+1) k
is the solution of the following equation:

n i=1 T t=0 • • • u i (t)D ik (X t )Σ -1 ik (X t )[Y s ik (t) -π s ik (t)] × p(X i | Y i , θ (p) )d(X i ) = 0.(4.2.14)
The derivation of β s k is given in Appendix (B.1). D ik (X t ) is a matrix with generic elements

D s ik (X t ) = π s ik (t)[1 -π s ik (t)] (4.2.15)
and Σ ik (X t ) = cov(Y ik (t)) has generic elements p+1) is the value for which the derivative of G 1 with respect to γ is nil. That is the

σ sm ik (t) = π s ik (t)[1 -π s ik (t)], if s = m -π s ik (t)π m ik (t) if s = m γ (
solution of ∂ ∂γ (G 1 ) = 0. Since ∂ ∂γ (G 1 ) = ∂ ∂γ n i=1 • • • T t=0 Z i (t)ν i (t) -b[ν i (t)] φ i (t) + c[Z i (t), φ i (t)] ×p(X i | Y i , θ (p) )d(X i ), = n i=1 • • • T t=0 ν i (t) φ i (t) + c [Z i (t), φ i (t)] × ∂Z i (t) ∂γ ×p(X i | Y i , θ (p) )d(X i ) = 0, (4.2.16)
where Z i (t) defined by equation (4.2.7), then, γ (p+1) is the solution of equation (4.2.16)= 0. Also, δ (p+1) is the value for which the derivative of G 1 + G 3 with respect to δ is nil. That is the solution of:

∂ ∂δ {G 2 + G 3 } = n i=1 • • • T t=0 ν i (t) φ i (t) + c [Z i (t), φ i (t)] ∂Z i (t) ∂δ × p(X i | Y i , θ (p) )d(X i ) - H [X i (t -1), u i (t), δ] H[X i (t -1), u i (t), δ] = 0.

Information matrix via the EM algorithm

In this section, as in chapter 3 we compute the Fisher information matrix via the Oakes' identity

First component: The first component is the second-order derivative with respect to β s k :

∂ 2 Q(β | β (p) ) ∂β s2 k = - n i=1 T t=1 • • • u i (t)D ik (p) (X t )Σ -1(p) ik (X t )D (p) ik (X t )u i (t) ×p(X i | Y i , θ (p) )d(X i ). (4.2.17)
The second-order derivatives with respect to γ and δ are given by :

∂ 2 Q(γ | γ (p) ) ∂γ 2 = n i=1 • • • T t=0 c [Z i (t), φ i (t)] × ∂ 2 Z i (t) ∂γ 2 ×p(X i | Y i , θ (p) )d(X i ), (4.2.18) ∂ 2 Q(δ | δ (p) ) ∂δ 2 = n i=1 • • • T t=0 c [Z i (t), φ i (t)] ∂ 2 Z i (t) ∂δ 2 × p(X i | Y i , θ (p) )d(X i ) - H [X i (t -1), u i (t), δ] × H[X i (t -1), u i (t), δ] -2H [X i (t -1), u i (t), δ] H 2 [X i (t -1), u i (t), δ] .(4.2.19)
Second component: The second component is computed as follows

∂ 2 Q(β | β (p) ) ∂θ (p) ∂β s k = n i=1 T t=1 • • • u i (t)D ik (X t )Σ -1 ik (X t )[Y s ik (t) -π s ik (t)] × ∂p(X i | Y i , θ (p) ) ∂θ (p) d(X i ). (4.2.20) ∂ 2 Q(γ | γ (p) ) ∂θ (p) ∂γ = n i=1 • • • T t=0 ν i (t) φ i (t) + c [Z i (t), φ i (t)] × ∂Z i (t) ∂γ × ∂p(X i | Y i , θ (p) ) ∂θ (p) d(X i ). (4.2.21)
and

∂ 2 Q(δ | δ (p) ) ∂θ (p) ∂δ = n i=1 • • • T t=0 ν i (t) φ i (t) + c [Z i (t), φ i (t)] ∂Z i (t) ∂δ ∂p(X i | Y i , θ (p) ) ∂θ (p) d(X i ). - H [X i (t -1), u i (t), δ] H[X i (t -1), u i (t), δ] . (4.2.22)
The derivative of the posterior density is done as in chapter 3.

Fisher information :

The information matrix (Fisher information) of parameter β can be written as

I(β) = - ∂ 2 Q(β | β (p) ) ∂β s2 k + ∂ 2 Q(β | β (p) ) ∂θ (p) ∂β s k , (4.2.23)
where

∂ 2 Q(β|β (p) ) ∂β s2 k and ∂ 2 Q(β|β (p) )
∂θ (p) ∂β s k are given by equations ( 4.2.17 ) and (4.2.20), respectively.

The information matrix of parameter γ is :

I(γ) = - ∂ 2 Q(γ | γ (p) ) ∂γ 2 + ∂ 2 Q(γ | γ (p) ) ∂θ (p) ∂γ , (4.2.24)
where

∂ 2 Q(γ|γ (p) ) ∂γ 2 and ∂ 2 Q(γ|γ (p) )
∂θ (p) ∂γ are given by equations ( 4.2.18) and (4.2.21), respectively.

The information matrix of parameter δ is given by :

I(δ) = - ∂ 2 Q(δ | δ (p) ) ∂δ 2 + ∂ 2 Q(δ | δ (p) ) ∂θ (p) ∂δ , ( 4

Consistency and Asymptotic normality of the maximum likelihood estimator(MLE)

We recall the parameters of our model is θ = (β , γ , δ ) , where

β = (β 1 , • • • , β q ) and β k = (β 1 k , • • • , β c k k ) , k = 1, • • • , q.
As in Chapter 3, we discuss the consistency and asymptotic normality to the model parameters.

Assumptions

Our assumptions remain the same as in chapter 3. Only (A4) changes to (A.4) The state space X i (t) follows equation (4.2.3) with state noise from exponential family distribution.

As in chapter 3, under the regularity assumptions (A.1-A.7) and the convergence properties to the EM algorithm (which had been explained in chapter 3) we conclude that

N 1/2 ( θ -θ) ∼ d N (0, I -1 (θ)). (4.2.26)

The Posterior distribution

As in chapter 3, the estimation formulas for parameters uses the the posterior distribution p(X i | Y i ). This density can be computed by the Auxiliary Iterated Extended Kalman Particle Filter (AIEKPF) method previously is presented in chapter 3. Therefore, its algorithm is recalled here.

Auxiliary Iterated Extended Kalman Filter (AIEKPF) Algorithm

The AIEKPF algorithm can be performed for the individual i in the following steps:

1. Initialization (t = 0) : For m = 1, • • • , N , draw the states (particles) X m i (0) from the prior p(X i (0)) ∼ f i X i (0) , and set

X m i (0) = E[X m i (0)] = F (x m i (0), u m i (0), δ m )
and 1) . 3. For m = 1, • • • , N , update the particles using the IEKF algorithm 3-1. Compute the Jacobians A m i (t), C m i (t) of the process model

P m i (0) = E[(X m i (0) -X m i (0))(X m i (0) -X m i (0)) ] = Var(X m i (0)) = H 2 (x m i (0), u m i (0), δ m )R 0 For t = 1, • • • , T, repeat the following steps: 2. For m = 1, • • • , N , generate µ m i (t) ∼ f i X i (t-
A m i (t) = ∂F (x, u i (t), γ) ∂x | x=x m i (t-1|t-1) C m i (t) = H(x m i (t -1 | t -1), u i (t), δ)
3-2. Predict the particle with the IEKF:

X m i (t | t -1) ≈ F (x m i (t -1 | t -1), u i (t), γ) P m i (t | t -1) = A m i (t)P m i (t -1 | t -1)A m i (t) + C m i (t)RC m i (t) 3-3. For j = 1, • • • c (c is the iteration number of the IEKF) a-Compute the Jacobians , B m ij (t) B m ij (t) = ∂π it (u i (t), x) ∂x | x=x m ij (t|t-1)
b-Update the state estimation error covariance P m ij (t) :

P ij (t | t) = (I -K ij (t)B ij (t)P ij (t | t -1) (4.3.1) K ij (t) = P ij (t | t -1)B ij (t)[B ij (t)P ij (t | t -1)B ij (t) + Σ -1 i (t)] -1
c-Update the state estimate X m ij (t) :

X ij (t | t) = X ij (t | t -1) + K ij (t)[Y i (t) -π i (t)] (4.3.2) 4. For m = 1, • • • N, calculate w m i (t) = q(m | Y i (t)) ∝ q k=1 M [π ik (u m i (t), µ m i (t))] w m i (t -1).
5. Resample to obtain the index ς m i of particle m's parent.

Importance sampling : for

m = 1, • • • , N , 6-1. Draw samples X m i (t) ∼ q(X i (t), ς m i | Y i (t)) = N ( X ς m i ij (t); P ς m i ij (t)), (4.3.3)
6-2. Calculate importance weights of particles by using

w m i (t) = p[Y i (t) | X m i (t)] p[Y i (t) | µ ς m i i (t)] (4.3.4) = q k=1 p[Y ik (t) | X m i (t)] q k=1 p[Y ik (t) | µ ς m i i (t)] = q k=1 M[π ik (x m i (t), u m i (t)] M[π ik (µ ς m i i (t), u m i (t
ϕ i (t) = N m=1 W m i (t)ϕ m i (t), (4.3.7)
where

W m i (t) = w m i (t)[ N m=1 w m i (t)] -1 , and
w m i (t) = p(Y i (t) | X m i (t)) p(Y i (t) | µ ς m i i ) = q k=1 M(π ik (u m i (t), X m i (t))) M(π ik (u m i (t), µ ς m i i ))
. σ 2 AP F (ϕ i (t)) as equations (3.4.20), (3.4.21) in chapter 3, but with q t (X i (t)

| X i (t -1)) = f i X i (t-1)

Posterior mode estimation

In this chapter, the observations are from a multinomial distribution and the latent variables are described by a CHARN model with noise from an exponential family distributions. In this chapter we extend the numerical solutions used in chapter 3 to find the posterior mode. We used the following two approaches :

1. Penalized likelihood estimation : a -Gauss-Newton and Fisher-scoring Filtering and smoothing algorithms. b -Working extended Kalman filter and smoother algorithms.

2. Maximum a posteriori sequence estimation (MAP) via the auxiliary iterated extended Kalman filter particle filter-MAP.

Penalized likelihood estimation

The posterior mode smoother is given by

a ≡ a (0 | T ), a (1 | T ), • • • , a (T | T ) ∈ R m ,
with m = (T + 1)n, where

a(0 | T ) = ((a 1 (0 | T ), • • • , a n (0 | T )) , . . . . . . . . . . . . a(T | T ) = (a 1 (T | T ), • • • , a n (T | T )) .
The posterior distribution of X by Bayes' theorem

p(X | Y) = 1 p(Y) n i=1 q k=1 T t=1 p(Y ik (t) | X i (t)) × n i=1 T t=1 g i (X i (t)) n i=1 g i (X i (0)) (4.4.1) p(Y) does not depend on X p(X | Y) ∝ n i=1 q k=1 T t=1 p(Y ik (t) | X i (t)) × n i=1 T t=1 g i (X i (t)) n i=1 g i (X i (0)) (4.4.2)
where

g i (X i (0)) = 1 H[X i (0), u i (0), δ] exp Z i (0)υ i (0) -b[υ i (0)] φ i (0) + c[Z i (0), φ i (0)] . (4.4.3)
Taking logarithms and inserting the densities of X i (t) in equation (4.2.6 ) and X i (0) in equation (4.4.3), the penalized log-likelihood function is given by

P L : R m -→ R, m = (T + 1)n P L(X) = n i=1 q k=1 T t=1 {log p(Y ik (t) | X i (t))} (4.4.4) + n i=1 T t=1 logg i (X i (t)) + n i=1 log g i (X i (0)) where log p(Y ik (t) | X i (t)) = log c k s=1 [π s ik (t)] Y s ik (t) = c k s=1 Y s ik (t) log π s ik (t). Then P L(X) = n i=1 q k=1 T t=1 c k s=1 {Y s ik (t) log π s ik (t)} + G 1 + G 2 (4.4.5)
where

G 1 = n i=1 Z i (0)υ i (0) -b[υ i (0)] φ i (0) + c[Z i (0), φ i (0)] -log(H[X i (0), u i (0), δ]) G 2 = n i=1 T t=1 Z i (t)υ i (t) -b[υ i (t)] φ i (t) + c[Z i (t), φ i (t)] -log(H[X i (t -1), u i (t), δ]).
Numerical maximization of the penalized log-likelihood can be achieved by various algorithms. As in chapter 3 we use two methods. First, iterative forward-backward Gauss-Newton (Fisher-scoring) algorithms. Second, Working Extended Kalman Filter and Smoother(WEKFS).

Gauss-Newton and Fisher-scoring Filtering and smoothing

The penalized log-likelihood criterion (4.4.4) can be presented in compact matrix notation as: P L(X) = l 1 (X) -l 2 (X), (4.4.6) where

l 1 (X) = n i=1 q k=1 T t=0 c k s=1 {Y s ik (t) log(π s ik (t)} , and l 2 (X) = X Aυ -1 Ab(υ) + c(X, φ)
where

• X is a (n × T ) matrix of latent variables.

• A = diag( 1 φ ) with size (n × n), with φ = 0.
• υ is the (n × T ) matrix of link function.

• 1 is a (n × T ) matrix of ones.

• b(υ), c(X, φ) are a (n × T ) matrices taking a different forms depending on the distribution of the X i (t).

For the description of Fisher scoring steps in matrix notation, the tables of observations for all i = 1, • • • , n, can be written as follows

Y i = (Y i (0), Y i (1), • • • , Y i (T )) , Y i (t) = (Y i1 (t) • • • , Y iq (t)) , t = 1, 2, • • • , T.
Fahrmeir and Wagenpfeil (1997) assumed Y i (0) = a i (0). Correspondingly we define the tables of expectations by

Π i (X) = (π i0 (X 0 ), π i1 (X 1 ), • • • , π iT (X T )) .
They assumed π i (0) = X i (0). We recall the conditional mean and variance of individual i

E(Y ik (t) | X i (t)) = π s ik (t), s = 1, • • • , c k Var(Y ik (t) | X i (t)) = Σ ik (X t ),
where Σ ik (X t ) has generic elements

σ sm ik (t) = π s ik (t)[1 -π s ik (t)], if s = m -π s ik (t)π m ik (t) if s = m (4.4.7)
The diagonal covariance matrix of an individual i at time t is

Σ i (X) = diag(V i (0), Σ i1 (X 1 ), • • • , Σ iT (X T )),
where V i (0) = H 2 [X i (0), u i (0), δ]R 0 , and the diagonal matrix

D i (X) = diag(1, D i1 (X 1 ), • • • , D iT (X T ))
Since D it (X t ) is first-order derivative of the conditional probability π i (t) evaluated at η i (t).

The score function of l(X) in (4.4.6) is given by S i (X) = ( S i0 (X 0 ), S i1 (X 1 ),

• • • , S iT (X T )) ,
where for all i = 1,

• • • , n, S i (X) := D i (X)Σ -1 i (X) {Y i (t) -Π i (X)} , (4.4.8)
with components

S i (X 0 ) = V -1 i (0)(a i (0) -X i (0)) (4.4.9) S it (X t ) = D it (X t )Σ -1 it (X t ) {Y i (t) -π it (X t )} , t = 1, • • • , T, (4.4.10)
the first-order derivatives of P L(X) in (4.4.6) are M (X) = ∂P L(X)/∂X = S(X) -S(υ), (4.4.11) where S(υ) = X A -1Ab (υ).

The expected information matrix is given by

I i (X) = (I i0 (X 0 ), I i1 (X 1 ) • • • , I iT (X T )),
where for all i = 1, • • • , n,

I i (X) = D i (X)Σ -1 i (X)D i (X) (4.4.12)
with diagonal blocks

I i0 (X 0 ) = V -1 i (0) (4.4.13) I it (X t ) = D it (X t )Σ -1 it (X t )D it (X t ), t = 1, • • • , T. (4.4.14)
As in chapter 3, the Taylor expansion of the score function about X 0 yields

M ( X | Y) ≈ M (X 0 ) -I(X 0 ) × X 1 -X 0
since M ( X | Y) = 0, a single Fisher scoring to the next iterate X 1 ∈ R m , with m = (T + 1)n is as follows

I(X 0 ) + J (ν 0 ) × X 1 -X 0 = M (X 0 ).
This can be rewritten as

X 1 = I(X 0 ) + J (ν 0 ) -1 I(X 0 ) Y, (4.4.15)
where

J (υ 0 ) = -1Ab (υ 0 ), with "working" observation Y = ( Y 1 , • • • , Y n ) and Y i = ( Y i (0), Y i (1), • • • , Y i (T )))
, then can compute as

Y i := D -1 i ((X)) {Y i -Π i (X))} + η i (X), (4.4.16)
with components

Y i (0) = a i (0) Y i (t) = D -1 it (X t ) {Y i (t) -π it (X t )} + η it (X t ), t = 1, • • • , T, where η i ((X)) = (η i1 (X 1 ), • • • , η iT (X T ))
, is the vector of link function for i individual.

Working Extended Kalman Filter and Smoother (WEKFS)

In the following algorithms, a t|t-1 , , a t|t , a t|T are numerical approximations to predicted, filtered, smoothed values of posterior mode to X and the corresponding P t|t-1 , P t|t , P t|T are numerical approximations to predicted, filtered, smoothed values of error covariance matrices.

Initialization:

a i (0 | 0) = a i (0), V i (0 | 0) = V i (0). ( 4
.4.17)

Prediction for t = 1, • • • , T a i (t | t -1) = F (a i (t -1 | t -1), u i (t), γ) P i (t | t -1) = A i (t)P i (t -1 | t -1)A i (t) + C i (t)R t C i (t). ( 4

.4.18)

where

A i (t) = ∂F (x, u i (t), γ) ∂x | x=a i (t-1|t-1)
.

Filtering for t = 1, • • • , T , a i (t | t) = a i (t) + K i (t)( Y i (t) -a i (t | t -1)) K i (t) = P i (t -1 | t -1)B i (t)(B i (t)P i (t -1 | t -1)B i (t) + I -1 (t)) -1 P i (t | t) = (I -K i (t)B i (t))P i (t -1 | t -1) (4.4.19)
where

B i (t) = ∂π i (t) ∂x | x=a i (t|t-1) .

The maximum a posteriori (MAP) with particle filtering and smoothing algorithm

As in chapter 3, we find the maximum a posteriori (MAP) via the Auxiliary Iterated Extended Kalman Particle Filter (AIEKPF) method. We present the algorithm with the state noise processes from exponential family distribution.

Auxiliary Iterated Extended Kalman Particle Filter -MAP(AIEKPF-MAP) Algorithm

We find the MAP estimation by using the Auxiliary Iterated Extended Kalman Particle Filter (AIEKPF), where:

w m i (t) = p[Y i (t) | X m i (t)] p[Y i (t) | µ ς m i i (t)] (4.4.20) = q k=1 p[Y ik (t) | X m i (t)] q k=1 p[Y ik (t) | µ ς m i i (t)] = q k=1 M[π ik (X m i (t), u i (t))] M[π ik (µ ς m i i (t), u i (t))]
The AIEKPF algorithm can be implemented by the following steps:

1. Initialization (t = 0) : For m = 1, • • • , N , draw the states (particles) X m i (0) from the prior p(X i (0)) = f i X i (0)
, and set

X m i (0) = E[X m i (0)] = F (X m i (0), u m i (0), δ m ) and P m i (0) = E[(X m i (0) -X m i (0))(X m i (0) -X m i (0)) ] = Var(X m i (0)) = H 2 (X m i (0), u m i (0), δ m )R 0 For t = 1, • • • , T, repeat the following steps: 2. For m = 1, • • • , N , generate µ m i (t) ∼ f i X i (t-1) .
3. For m = 1, • • • , N , Update the particles using the IEKF algorithm 3-1. Compute the Jacobians A m i (t), C m i (t) of the process model

A m i (t) = ∂F (x, u i (t), γ) ∂x | x=X m i (t-1|t-1) C m i (t) = H(X m i (t -1 | t -1), u i (t), δ)
3-2. Predict the particle with the IEKF:

X m i (t | t -1) ≈ F (X m i (t -1 | t -1), u i (t), γ) P m i (t | t -1) = A m i (t)P m i (t -1 | t -1)A m i (t) + C m i (t)R t C m i (t) 3-3. For j = 1, • • • c (c is the iteration number of the IEKF) a-Compute the Jacobians , B m ij (t) B m ij (t) = ∂π it (u i (t), x) ∂x | x=x m ij (t|t-1)
b-Update the state estimation error covariance P m ij (t) :

P ij (t | t) = (I -K ij (t)B ij (t)P ij (t | t -1) (4.4.21) K ij (t) = P ij (t | t -1)B ij (t)[B ij (t)P ij (t | t -1)B ij (t) + Σ -1 i (t)] -1 c-Update the state estimate X m ij (t) : X ij (t | t) = X ij (t | t -1) + K ij (t)[Y i (t) -π i (t)] (4.4.22) 4. For m = 1, • • • N, calculate w m i (t) = q(m | Y i (t)) ∝ q k=1 M(π ik (u i (t), µ m i (t))w m i (t -1), 5.
Resample to obtain the index ς m i of particle m's parent.

Importance sampling : for

m = 1, • • • , N , 6-1. Draw samples X m i (t) ∼ q(X i (t), ς m i | Y i (t)) = N ( X ς m i ij (t), P ς m i ij (t))
, where j = c 6-2. Calculate importance weights of particles by using

w m i (t) = p[Y i (t) | X m i (t)] p[Y i (t) | µ ς m i i (t)] = q k=1 M[π ik (x m i (t), u i (t))] M[π ik (µ ς m i t (t), u i (t))] 6-3. Normalize the weights w m i (t) = w m i (t) N m=1 w m i (t) .
Forward Filtering step

X M AP i (t) = arg max X (m) i (t) p(Y i (t) | X (m) i (t)) j p(X (j) i (t) | X (j) i (t -1))w (j) i (t -1) = arg max X (m) i (t) q k=1 p(Y ik (t) | X (m) i (t)) j p(X (j) i (t) | X (j) i (t -1))w (j) i (t -1) = arg max X (m) i (t) q k=1 M(π m ik (t)) j f i X i (t-1) .w (j) i (t -1) . Backward smoothing step • Set w (m) i (T | T ) = w (m) i (t)
• For t = T -1, • • • , 1, evaluate the smoother importance weights as

w (m) i (t | T ) = w (m) i (t) N j=1 w (j) i (t + 1 | T ) p(X (j) i (t + 1) | X (m) i (t)) N r=1 p(X (j) i (t + 1) | X (r) i (t))w (r) i (t) = w (m)(t) i N j=1 w (j) i (t + 1 | T ) f j X i (t-1) N r=1 f r X i (t-1) )w (r) i (t)
.

• Evaluate the approximate smoother MAP as

x M AP i (t | T ) = arg max X (m) i (t) q(X (m) i (t) | Y i (t)) w (m) i (t | T ) w (m) i (t) . = arg max X (m) i (t) q k=1 M(π m ik (t)) j f j X i (t-1) .w (j) i (t -1) w (m) i (t | T ) w (m) i (t)
.

Chapter 5

Practical considerations

Introduction

In this chapter, the numerical results and their discussion are presented. There are two parts. In the first part, the numerical simulation where the longitudinal multicategorical data are generated with latent variables described by a CHARN model and the state noise is generated from exponential family distributions. In the second part, we apply our results to real data in the quality of life of patients surged for cancer.

In the two parts, we created R codes from the methods described in chapters 3 and 4. Our objective in this chapter is to estimate the latent variables by posterior mode via the working extended Kalman filtering recursions.

Simulation experiments

In this part, we produce data from the observation equation described by a multinomial distribution defined in equation (3.3.1), and the state equation described by a CHARN model defined in equation (3.3.3) with state noise from exponential families. We consider gaussian and exponential distributions. There are two parts. The first aims at testing the efficiency of the working extended Kalman filter recursions. Here, the parameters of the models are assumed to be known. The second part uses the EM algorithm for estimates the parameters of the model, before applying the working extended Kalman filter recursions.

Simulation experiments I

In order to investigate the efficiency of the working extended Kalman filter recursions (WEKF), we consider that an individual fills out a questionnaire constituted of multiple choice questions administered at t occasions. The outline of simulation experiments I is as follows

The outline of simulation I

In this outline we have one individual, so we omit the subscript i in the equations.

1. Generate the multicategorical longitudinal data as follows:

(a) Suppose one individual in the longitudinal study to whom is asked a set of (q = 5) items, for each item k, (c k = 6) categories administered at t occasions.

(

π s k (t) = exp[η s k (t)] 1 + c k j=1 exp[η j k (t)] where η s k (t) = u (t)β s k + X(t). (

Simulation experiments II

In this part, the longitudinal data are generate, and it is considered as real data. In other words, after we generate the data, the state variables and the parameters of models are considered as unknown. In the following, we outline how we estimate the latent variables.

Recall that

Y i = (Y i (0), Y i (1), • • • , Y i (T )) , X i = (X i (0), X i (1), • • • , X i (T )) .
Calculate the posterior distribution p(X i | Y i ) via the Auxiliary Iterated Extended Kalman Particle Filter (AIEKPF) algorithm.

3. Set iteration p = 0, and apply the classical Kalman Filtering Recursions to calculate the initial value a 0 i (t) of posterior mode.

4. Starting with a 0 i (t), compute the model's parameters β p+1 , γ p+1 , δ p+1 via EM algorithm.

ε i (t) -E(ε i (t)) Var(ε i (t))
In steps [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF][START_REF] Bartholomew | Latent variable models and factor analysis: A unified approach[END_REF][START_REF] Bartolucci | Modeling Longitudinal Data by Latent Markov Models with Application to Educational and Psychological Measurement[END_REF][START_REF] Bartolucci | Longitudinal analysis of selfreported health status by mixture latent auto-regressive models[END_REF], the initial values are supposed λ (0) = 0.6, and β (0) as the Gaussian distribution.

To apply the working Kalman filtering recursions, we supposed ρ (p) = 1/λ (p) and R (p) t = 1/λ (p)2 . As the previous section, there are two parts. First, (n ≤ 15) with occasions (T > 75). Second (n > 75) with occasions (T ≤ 10).

• n small with T large As the gaussian state noise we consider n = 1, 2, 3, 4, 10, and T = 75, 100, 200, 300. As seen in Figures 5.17

• n large with T small

In this part,we consider n = 75, 100, 200, 300, and T = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. As the gaussian state noise at instant T = 1, 2, 3, 4, the red circle refers to the latent variables, the blue circle refers to the posterior mode via the filtering step is , and the posterior mode via the prediction step is the green circle. show at instant T = 1, 2, 3, 4, the results are excellent, but at T = 5, 10 there are some of the spaces in the curves, but the results are acceptable.

Discussion of the results

A disadvantage of EM algorithm is its slow convergence as we have seen when applying the simulation II. [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] proved that converges only linearly to the maximum likelihood estimate. Then we used the numerical analysis through the Fisher scoring method. The Fisher scoring method is an iterative algorithm for obtaining the maximum of likelihood. It started with initial values for the parameters. The iterations continue until the maximum likelihood estimates are achieved.

The key for performing the working Kalman filter recursions is using the parameters which are calculated in the EM algorithm as initial values to estimate of the posterior mode. In some cases, we observed that the calculated values of the parameters achieved the condition of convergence for likelihood in equation ( 5.2.2). But it is not optimal initial values for the working extended Kalman filtering. So this explains the existence of some unsatisfactory results. We conclude that the optimal choice for the initial values of the parameters leads to the best estimators for the state variables. The patients filled the questionnaire several times, the patients' questionnaires were completed at their clinic visits, or they were given a postage-paid envelope to return the questionnaires. The questionnaires considered personality traits which was completed before the surgery. The LOT "life Orientation Test" questionnaire and the trait section of the STAI-B " State-Trait Anxiety Inventory" instrument were used, the Table 5.1 shows that the French version (STAI-B) consists of twenty items. The patient responses to each item are classified into 4 categories (almost never, sometimes, often, and almost always). Here, the latent variable is the patient fatigue after surgery. This variable is assumed to be quantitative and varying over time around a mean value assumed to be nil in our work. Ten covariates are determined for the study : age, marital status, family situation, number of children and children's age, education, employment status, the group chemotherapy, the step of treatment, and the distance between patient's home and hospital are collected at the baseline assessment.

Application to real data

Data analysis

We obtained data from the above mentioned centres. Some of the patients' informations are missed (Not registered). So these patients' informations are removed from the data, and it remains 435 patients with complete information. Each patient is asked to fill the questionnaire at 10 instants. The questionnaire has 20 items. We selected 3 covariates: the marital status, the step of treatment, and the group of chemotherapy. These are said by the expert to be related to the latent variable studied that is the fatigue of the patient. They are: marital status are scored : 1= "single", 2= "cohabitation", 3= "bride", 4="widow", 5="divorced", and 6= "bride/cohabitation". The scores of the step of treatment are : 1= "step I", 2= "step II", and 3= "step III". The scores of the group chemotherapy are : 1= " the group without chemotherapy ", and 2= "the group with chemotherapy ". From Figure 5.32 we can see that the individual 1 is tired at T = 1, and he feels rested at T = 2, 3, 4. Then, he returns tired at T = 5, while he feels rested at T = 6, 7, 8, 9, 10. The individual 6 is tired at T = 1, and he feels rested at T = 2, 3, 4. Then, he is tired at T = 5, but he feels rested at T = 6, 7, whereas at T = 8 he is tired. Finally, at T = 9, 10 he feels rested.

CHARN(1,1) model

As in the preceding subsection, the outline of the analysis of data is performed with CHARN(1,1) as in equation (5.2.5). We set the initial values for β 0 and R 0 t as in the AR(1) model, and set ρ 0 1 = 0.6, ρ 0 2 = 0.4. Figures 5.37 -5.41 display the filter and predictive steps for the first ten individuals of every hundred individuals. Figures 5. 37-5.41 show that the curves of predictive and filter steps for each individual with CHARN(1,1) model are the same as the curves for AR(1) model. In other words, the effect of CHARN(1,1) is similar to that the AR [START_REF] Agresti | Categorical Data Analysis[END_REF]. We can explain that the predictive and the filter values depend on the function F = ρX i (t -1), which is supposed that it is the same in two model. The variance-covariance matrix for predictive, filter P t|t-1 , P t|t respectively also have values less than 0.001. On Figures 5.37 -5.41 too, positive values the rested patient and negative values the tired patient. Figure 5.37 shows the curves of prediction and filter steps for the 1st individual to the fifth individual are identical. We can also see that the individual 1 is tired at T = 1, and he feels rested at T = 2, 3, 4, and he returns tired at T = 5. Finally, he feels rested at T = 6, 7, 8, 9, 10. The curves of prediction and filter steps for the sixth individual to the tenth individual are similar. We also see that the individual 6 is tired at T = 1, and he feels rested at T = 2, 3, 4. Then, he is tired at T = 5, but he feels rested at T = 6, 7, while at T = 8 he is tired. Eventually, at T = 9, 10 he feels rested. This thesis focuses on the generalization of Kalman filter and smoother recursions with generalized state space models. We took a practical model, the observations are multi-categorical and longitudinal, and the state equation is described by a CHARN models with noise from exponential family distributions. This chapter presents the contributions of the study, conclusions and some recommendation for future studies.

Conclusion

In this thesis, we have investigated and discussed the generalization of Kalman filter and smoother recursions with generalized state space models. The main purpose of this study is to estimate the state by generalized Kalman filter recursions. This thesis also discusses the Kalman filter recursion with linear and nonlinear state space model, EM algorithm for estimating the model parameters. We propose two approaches in this study to estimating the state variables. The first, Working Extended Kalman Filtering Recursion (WEKFS). The second, the maximum a posteriori (MAP) via auxiliary Iterated Kalman particle filter. To achieve this objective, we apply the working extended Kalman filtering recursions with various simulation studies have been designed. The simulation studies are designed with different characteristics such as size of samples (number of individuals and length of occasions), the form of the distribution of state noise (gaussian or exponential), and the type of state equation (AR(1) model, CHARN [START_REF] Agresti | Categorical Data Analysis[END_REF][START_REF] Agresti | Categorical Data Analysis[END_REF], and CHARN(0,1)). In addition to the simulation, an application to real data in the quality of life of patients surged for breast cancer. In order to evaluate our approach in this study, we present the figures that appear the state variables which are simulated with their predictive and filter values. Another criteria we use the P t|t-1 , P t|t are the variance-covariance matrix to predictive and filter respectively, where each matrix have values less than 0.001 refer to the best estimator. The first objective is to propose a new multicategorical longitudinal in quality of life, where we use the state space model approach. The observation equation is a conditional probability of Y ik (t) given the state variables X i (t) described by a multinomial distribution. The state equation is described by conditional heteroskedastic autoregressive nonlinear (CHARN) model. The second objective of this study is to propose an alternative and generalisation of existing methods for estimating the latent trait in quality of life.

Areas of Future Studies

In this section, some future recommendations would be mentioned based on our current study, which can play an important role in the future research. We can summarize the suggested topic for future studies as follows.

• Our approach succeeds for estimating a latent variable in the health field.

Consequently, our approach can be applied in the following fields: the economic field for estimating the business confidence or morale of customers, in the industrial field for estimating the level of anxiety due to the machines or robots on workers in factories.

• In chapter 5, we apply the first approach the working Kalman filter recursion through created the R-packages, we can create the R-package for the maximum a posteriori (MAP) via the Auxiliary Iterated Kalman Particle Filter (AIEKPF) method. One can compare between two approaches which are proposed.

• In real data analysis, we found the data Incomplete, where missed the covariates or the responses for the individuals. We found challenging to analyse this data. Therefore we ignored the individuals have missing data. One can develop our approach to estimating the latent trait for the Incomplete data.

• In our study, we proposed the link function is linear predictor η s ik (t) = u i (t)β s k + X i (t). Hastie and [START_REF] Tibshirani | Generalized additive models[END_REF] introduced the class of generalized additive models which replaces the linear form by a sum of smooth function j s j (u j ), where s j (.) s are unspecified functions that are estimated using a scatterplot smoother. By the similar way one can be developed our model with η s ik (t) = s i (u i (t)) + X i (t).

where, p(X t | Y t ) and p(X t+1 | Y t ) are the filtering density and one step ahead predictive density respectively, at time t. Thus, starting with p(X T | Y T ), p(X t | Y T ) can be recursively obtained from p(X t+1 | Y T ). Using the above recursion, the marginal smoothing distribution can now be approximated by the weighted particle cloud. Here, one starts with the forward filtering pass for computing the filtered distribution at each time step using particle filter as eq (3.4.1), then one performs the backward smoothing pass as given by (A. T . It is important to note that the forward-backward smoother keeps the same particle support as used in filtering step and re-weights the particles to obtain the approximated particle based smoothed distribution. To obtain the marginal MAP smoother, one needs the posterior density p(X t | Y T ) from the above cloud representation. Here, we proceed as follows. Using the Bayes' rule, one can write the one step ahead predictive density in equation (A. [START_REF] Bartolucci | Modeling Longitudinal Data by Latent Markov Models with Application to Educational and Psychological Measurement[END_REF] The likelihood is rewritten by x k = 1, = exp q-1 k=1

l(β) = q k=1 T t=1 c k -
Dev = 2 n i=1 {[x i νi -x i νi ] -[b(ν i ) -b(ν i )]} = 2 n i=1 {[x i log x i 1 -x i -x i log πi 1 -πi ] -[log 1 1 -x i -log 1 1 -πi ]} = 2 n i=1 {x i log x i πi -x i log 1 -x i 1 -πi + log 1 -x i 1 -πi } = 2
x k log( π k π q ) + log(π q ) .

where q k=1 π k = 1,

ν = log π 1 π q , log π 2 π q , • • • , log π q π q , φ = 1, b(ν) = -log (π q ) , c(x, φ) = 0 Dev = 2 n i=1 {[x i νi -x i νi ] -[b(ν i ) -b(ν i )]} = 2 n i=1 { q k=1
x ik log x ik -q k=1

x ik log πik } setting µ i (t) = X i (t) gives the maximum likelihood achievable log likelihood, namely l(X i (t), V i (t), X i (t)) = -T 2 log(2πV i (t)), so that the scaled deviance function is

Dev * (X i ; µ i ) = 2{l(X i (t), V i (t), X i (t))-l(X i (t), µ i (t), V i (t)} = T i=1

(X i (t) -µ i (t)) 2 V i (t) , or Dev * (X i ; µ i ) = 2 T t=1 X i (t)[ν i (t) -νi (t)] -b[ν i (t)] + b[ν i (t)] φ i (t) = 2 T t=1 X i (t)[X i (t) -μi (t)] -1 2 X 2 i (t) + 1 2 μ2 i (t) V i (t) = 2 T t=1 1 2 X 2 i (t) -X i (t)μ i (t) + 1 2 μ2 i (t) V i (t) = T t=1 X 2 i (t) -2X i (t)μ i (t) + μ2 i (t) V i (t) = T t=1 [X i (t) -μi (t)] 2 V i (t)
and

Dev(X i , µ i ) = V i (t) × T t=1 [X i (t) -μi (t)] 2 V i (t) (C.2.3) = T t=1 [X i (t) -μi (t)] 2 Then φd = Vi (t) = T t=1 [X i (t) -μi (t)] 2 T -2
Example 2: we have for an individual

X i (t) ∼ Gamma(α, λ α ),
where ν i (t) = -1 λ i (t) , φ i (t) = 1 α i (t) , b[ν i (t)] = -log[ν i (t)], µ i (t) = α i (t) λ i (t) , V i (t) =

α i (t) λ 2 i (t) ,

4 . 4 Definition 3 . 1 (

 4431 Optimality of the Auxiliary Iterated Extended Kalman Particle Filter(AIEKPF) Algorithm Johansen and Doucet (2008) proposed a novel interpretation of the APF, which allows them to present the first convergence results for APF algorithm. They defined the consistency and asymptotic normality of a weighted particles (samples) [{X m t , w m t } N m=1 ]. Asymptotic normality). Under the regularity condition given in [Chopin (2004), Theorm 1] or [Del Moral, 2004, Section 9.4, pp.300-306], one has the convergence in distribution

Smoothing step:

  [START_REF] Saha | On the Monte Carlo marginal MAP estimator for general state space models[END_REF] extended the MAP estimation concept to the marginal smoothing

Figure 2 . 2 :

 22 Figure 2.2: The Path of longitudinal Rasch processes

.4. 9 )

 9 Substituting Eqs.(3.4.3 ) and (3.4.5) into Eq.(3.4.9), one obtains

3. 4 . 2 Definition 3 . 1 (

 4231 Optimality of the Auxiliary Iterated Extended Kalman Particle Filter(AIEKPF) Algorithm Under the following regularity condition C.1, C.2 and C.3 given in [Chopin (2004) , Douc and Moulines (2008) ] one can prove the consistency and asymptotic normality of a weighted particles (samples) [{x m t , w m t } N m=1 ]. C.1 The initial sample to an individual [{X m i (0), w m i (0)} N m=1 ] is consistent and asymptotically normal. C.2 The selection step consists in multinomial resampling. C.3 φ and σ 2 AP F (ϕ) are finite quantities, where φ is the expectation of a function ϕ and σ 2 AP F (ϕ) the variance . Asymptotic normality). Under the conditions (C.1-C.3) a sample for one individual [{X

3. 5 . 4 . 2

 542 Smoothing step Saha et al. (2008) extended the MAP estimation concept to the marginal smoothing

Figure 5 .

 5 Figure 5.1 shows one individual with different type of models and the state noise distribution. There the red colour refers to the latent variable, while the blue colour refers to the postrior mode via filtering step and green colour refers to the posterior mode via the prediction step. It is clear that the results are very good. The working Kalman filter recursion succeeds in producing the posterior mode with different type of state space model. Moreover, the values via two steps are equal to the actual value of the latent variables. There is no curve for the posterior mode via the predictive step in Figures 5.1(e) and (f) because the prediction step does not exist for the model considered. Indeed, it depends on the function F of the state equation, this function does not exist for the CHARN(0,1) model studied.

  Figure 5.1 shows one individual with different type of models and the state noise distribution. There the red colour refers to the latent variable, while the blue colour refers to the postrior mode via filtering step and green colour refers to the posterior mode via the prediction step. It is clear that the results are very good. The working Kalman filter recursion succeeds in producing the posterior mode with different type of state space model. Moreover, the values via two steps are equal to the actual value of the latent variables. There is no curve for the posterior mode via the predictive step in Figures 5.1(e) and (f) because the prediction step does not exist for the model considered. Indeed, it depends on the function F of the state equation, this function does not exist for the CHARN(0,1) model studied.

Figure 5 . 1 :

 51 Figure 5.1: The graphs for n = 1individual and T = 100

  (a) AR(1) with noise gaussian (b) AR(1) with noise exponential (c) CHARN(1,1) with noise gaussian (d) CHARN(1,1) with noise exponential (e) CHARN(0,1) with noise gaussian (f) CHARN(0,1) with noise exponential

Figure 5 . 2 :

 52 Figure 5.2: The graphs for n = 1 individual (a) T = 75 (b) T = 100

Figure 5 . 3 :

 53 Figure 5.3: The graphs for n = 2 individuals (a) T = 75 (b) T = 100

Figure 5 . 4 :

 54 Figure 5.4: The graphs for n = 3 individuals (a) T = 75 (b) T = 100

Figure 5 . 5 :

 55 Figure 5.5: The graphs for n = 4 individuals (a) T = 75 (b) T = 100

Figure 5 . 6 :

 56 Figure 5.6: The graphs for n = 10 individuals (a) T = 75 (b) T = 100

i

  (t) and the initial values of parameters β (0) , γ (0) , δ (0) .

Figure 5 . 7 : 2 Figure 5 . 8 :

 57258 Figure 5.7: The graphs for n = 75 individuals with T = 1, 2 (a) T = 1 (b) T = 2

Figure 5 . 9 : 4 (Figure 5 . 10 : 10 Figure 5 . 11 : 2 Figure 5 . 12 :

 594510105112512 Figure 5.9: The graphs for n = 100 individuals with T = 1, 2, and 4 (a) T = 1 (b) T =

Figure 5 . 13 : 4 (= 4 Figure 5 . 14 :

 51344514 Figure 5.13: The graphs for n = 300 individuals with T = 1, 2, 3 and 4 (a) T = 1 (b) T = 2

Figure 5 . 15 :Figure 5 . 16 :

 515516 Figure 5.15: The graphs for n = 10 individuals (a) T = 75 (b) T = 100

  -5.18 some sharp jumps are appearing in the curve of latent variables, this leads to distancing values of posterior mode from the latent variables. Figures 5.19-5.21 show the results are good.

Figure 5 . 17 :

 517 Figure 5.17: The graphs for n = 1 individual (a) T = 75 (b) T = 100

Figure 5 . 18 :

 518 Figure 5.18: The graphs for n = 2 individuals (a) T = 75 (b) T = 100

Figure 5 . 19 :

 519 Figure 5.19: The graphs for n = 3 individuals (a) T = 75 (b) T = 100

Figure 5 . 20 :

 520 Figure 5.20: The graphs for n = 4 individuals (a) T = 75 (b) T = 100

Figure 5 . 21 :

 521 Figure 5.21: The graphs for n = 10 individuals (a) T = 75 (b) T = 100

Figure 5 . 22 : 4 (Figure 5 . 23 : 10 Figure 5 . 24 : 2 Figure 5 . 25 :

 5224523105242525 Figure 5.22: The graphs for n = 75 individuals with T = 1, 2, and 4 (a) T = 1 (b) T =

Figure 5 . 26 : 4 (= 4 Figure 5 . 27 : 10 Figure 5 . 28 : 2 Figure 5 . 29 :

 52644527105282529 Figure 5.26: The graphs for n = 200 individuals with T = 1, 2, 3 and 4 (a) T = 1 (b) T = 2

Figure 5 . 30 :

 530 Figure 5.30: The graphs for n = 2 individuals (a) T = 75 (b) T = 100

Figure 5 . 31 :

 531 Figure 5.31: The graphs for n = 3 individuals (a) T = 75 (b) T = 100

Rotonda

  et al. (2011) presented a study on quality of life. Their aim investigate factors correlated with cancer-related fatigue for women surgered for cancer. In this study, 502 patients were recruited from September 2008 until September 2010. Three French cancer centres received the patients, the Alexis Vautrin anti-cancer centre of Lorraine, the Georges-François Leclerc anti-cancer centre of Burgundy and the Paul Strauss anticancer centre of Alsaca, France.

Figure 5 . 32 :Figure 5 . 33 :Figure 5 . 34 :Figure 5 . 35 :Figure 5 . 36 :

 532533534535536 Figure 5.32: The graphs for the 1st individual to 10th individual with AR(1) model

Figure 5 . 37 :Figure 5 . 38 :Figure 5 . 39 :Figure 5 . 40 :Figure 5 . 41 :

 537538539540541 Figure 5.37: The graphs for the 1st individual to 10th individual with CHARN(1,1) model

4 . 1 ) 1 )

 411 to approximate the smoothing distribution p(X t | Y T ) = Where δ denoted the Dirac delta function. The smoothing weights are obtained through the following backward recursion:

C. 1 . 2 1 r- 1 .νC. 1 . 3

 121113 Negative Binomial distribution p(x) = x-1 r-1 π r (1 -π) x-r , x = {r, r + 1, • • • }, = exp x log(1 -π) + r log π 1 -π + log x-= log(1 -π), φ = 1, b(ν) = -r log π 1 -π , c(x, φ) = log( i νi -x i νi ] -[b(ν i ) -b(ν i )]} = 2 n i=1 {[x i log(1 -x i ) -x i log(1 -πi )] -[-r i log x i 1 -x i + r i log πi 1 Bernoulli distribution p(x) = π x (1 -π) 1-x , x ∈ {0, 1}, = exp x log( π 1 -π ) + log(1 -π) . ν = log π 1 -π , φ = 1, b(ν) = -log(1 -π), c(x, φ) = 0

n i=1 {x i log x i πi + ( 1 -

 1 x i ) log 1 -x i 1 -πi }

1 , 2 , 1 -(( 1 - 2 ) -x 2 /(2σ 2 ) 2 ) 2 ) -log x -log x 2 /(2σ 2 )x 2 i 2 Dev = 2 n i=1 {

 12112222222222i=1 • • • , n}, = exp {x log λ -λ -log x!} . ν = log λ, φ = 1, b(ν) = λ, c(x, φ) = -log(x!) i νi -x i νi ] -[b(ν i ) -b(ν i )i log x i -x i log λi ] -[x i -λi ]π k ) log x k + log[Γ( Γπ k )] . ν = (1 -π k ), φ = 1, b(ν) = log[Γ( Γπ k )], c(x, φ) = 0 Dev = 2 n i=1 {[x i νi -x i νi ] -[b(ν i ) -b(ν i )x ik ) log x ik -(1 -πik ) log x ik ] i νi -x i νi ] -[b(ν i ) -b(ν i )-log x -log x 2 2σ 2 ] . ν = µ, φ = σ 2 , b(ν) i νi -x i νi ] -[b(ν i ) -b(ν i )-2 log x i μi + μ2 i } = n i=1 (log x i -μi ) 2 log(2π) -3 log x 2 . ν = -1/µ 2 , φ = 2/λ, b(ν) = -2 µ , c(x, φ) = -λ 2x + log λ -log(2π) -3 log x [x i νi -x i νi ] -[b(ν i ) -b(ν i )α -1) log x + α log(α) -α log(λ) -log(Γ(α)) = exp -x( 1 λ ) -log(λ) 1 α + (α -1) log x + α log α -log(Γ(α)) . ν = -1 λ , φ = 1 α, b(ν) = log(λ), c(x, φ) = (α-1) log(x)+α log α-log(Γ(α))

  Xi et al. (2015) ont proposé un nouvel algorithme du filtre particulaire. Ils ont généré une densité d'importance par l'utilisation de l'algorithme du filtre de et Wagenpfiel (1997) et Fahrmer et Tutz (2013). Ceux-ci ont considéré que la série observée {Y t } est discrète et ont estimé, par maximum de vraisemblance pénalisé, la variable d'état par le mode a posteriori. Leur approche peut également être interprétée comme une méthode nonparamétrique pour le modèle espace-état. Ils ont estimé le mode a posteriori en utilisant la méthode du score de Fisher via le filtrage et le lissage itératif de Kalman. Deuxièmement, Bousseboua et Mesbah (2010) ont proposé une nouvelle classe de processus longitudinaux multivariés qui sont des résultats dichotomiques (Y ik (t)

	Kalman étendu itéré.
	L'estimation du Maximum a Posteriori (MAP) est utilisé pour obtenir une es-
	timation ponctuelle des variables latentes X t par maximum de vraisemblance
	:
	X M AP t	= arg max

Xt p(X t | Y t ).

Dans ce chapitre, nous présentons deux applications du modèle espace-état non-Gaussien. Nous présentons d'abord les modèles espace-état dont la loi du bruit appartient aux familles exponentielles (modèles linéaires généralisés). Ces modèles ont été étudiés par Fahrmeir

1 .

 1 Les simulations : Dans cette partie, nous produisons des données à partir de l'équation d'observation décrite par une distribution multinomiale définie par l'équation (3.3.1), et l'équation d'état décrite par un modèle CHARN défini par l'équation (3.3.3) avec bruit de loi appartenant à une famille exponentielle. Nous considérons les lois gaussienne et exponentielle. Il y a deux parties. La première partie vise à tester l'efficacité des récursivités du filtre de Kalman étendu, où nous considérons les paramètres du modèle connus. La deuxième partie utilise l'algorithme EM pour estimer les paramètres du modèle, avant d'appliquer des récursivités du filtre de Kalman étendu. Pour simuler, nous tenons compte de la

taille des échantillons (nombre d'individus et durée), la forme de la distribution du bruit d'état ( Gaussien ou exponentiel) et le type d'équation d'état (modèle AR (1), CHARN (1,1) et CHARN (0,1)). Afin d'évaluer notre approche dans cette étude, nous présentons les figures montrant la série des variables d'état simulées et celles de leurs valeurs obtenues par les algorithmes de prédiction et de filtrage.

  [START_REF] Bartolucci | Longitudinal analysis of selfreported health status by mixture latent auto-regressive models[END_REF] proposed a model for longitudinal categorical data. The latent process was modelled by a mixture of autoregressive AR[START_REF] Agresti | Categorical Data Analysis[END_REF] processes with different means and correlation coefficients, but with equal variance. They considered an application based on a longitudinal data set, concerning self-reported health status (SRHS) derived from the Health and Retirement Study (HRS) and conducted by the University of Michigan and supported by the US National Institute on Aging and the Social Security Administration.[START_REF] Bartolucci | Modeling Longitudinal Data by Latent Markov Models with Application to Educational and Psychological Measurement[END_REF] reviewed a class of models for longitudinal data. In this type of models, the unobserved individual characteristics of interest was represented by a sequence of discrete latent variables, which follows a Markov chain. He used two applications to illustrate his models. The first application is the evolution of the ability level in mathematics. There, Bartolucci considered testing the hypothesis of absence of tiring or learning-through-training phenomena during the administration of series of 12 items on Mathematics. The second application is the evolution of psychological traits in children. He considered data collected through a mental experiment based on tests administered at different occasions to pre-school children to measure two types of ability: inhibitory control and attentional flexibility.Fahrmeir and Wagenpfiel (1997) and[START_REF] Fahrmeir | Multivariate statistical modelling based on generalized linear models[END_REF] studied non-linear time series or discrete longitudinal observations. They created a method of inference based on the posterior mode or, equivalently, maximum penalized likelihood estimation. They acquired efficient smoothing by using the working Kalman filtering and smoothing. They analysed data taken from the IFO Institute for economic research in Munich which collects categorical monthly data of firms in various industrial branches.

Table 2 . 1 :

 21 The distributions and their approximations in the sequential Monte Carlo filter and smoother algorithm.

	Distributions Density function Approximations by particles
	Predictive distribution	p(X t | Y t-1 )	{p	(1)

  .[START_REF] Bousseboua | Processus de markov longitudinal latent rasch observable[END_REF] where S t the inverse function of G t in the Non-linear observation equation (2.2.3), and r t is the probability density function of the observation noise ξ t . Here α

	(j)

  .3.4.3 The Auxiliary Iterated Extended Kalman Particle Filter (AIEKPF)[START_REF] Xi | The auxiliary iterated extended Kalman particle filter[END_REF] proposed a new particle filter called the auxiliary iterated extended Kalman particle filter (AIEKPF). The basic idea of the AIEKPF algorithm is the use of the iterated extended Kalman filter (IEKF) algorithm within an APF algorithm to generate the importance density function. The proposed approach uses the iterated extended Kalman filter (IEKF) to merge the last measurements into state distribution. The robustness of the APF and the importance density generated by the IEKF leads to the new filter that can match the posterior distribution well. We adapted this algorithm to the non-linear non-gaussian models defined by (2.2.3,2.2.4) by deriving the equations of extended Kalman filtering. The derivation of equations for extended Kalman filtering recursions is presented in the Appendix (A.2). Let {X m t , m = 1, 2, • • • , N } be a set of support points with corresponding weights {w m

t , m = 1, 2, • • • , N }.

Remark: In the state space representation, the posterior densities {p(x 1:t | y 1:t ), t ≥ 1} are sequentially approximated by the Sequential Monte Carlo methods (SMC). The posterior distributions that are approximated by a set of N weighted random sample called particles, for t ≥ 1 leads to the approximation equation,

  Update the state estimation error covariance V m tj by Eq.( 2.3.33) c -Update the state estimate X m tj by Eq.(2.3.32).

  and Fahrmer and Tutz (2013). They considered the discrete time series observation {Y t } and estimated the state variable by posterior mode. Second, Bousseboua and Mesbah (2010) proposed a new class of longitudinal multivariate processes which are dichotomous outcomes (Y ik (t

  For example in binomial distribution h(π t ) = exp(η t )/{1 + exp(η t )} where η t = u i β, with (u i ) standing for the covariate variables and β the regression parameters. This function h is known as the logit model. In contrast, h(π t ) = Φ(π t ) is called the probit model.

.

4.4)

where h : R → R is a response function which is two-times continuously differentiable. Z t is a × k matrix, which may be expressed by covariates u t or also by past responses Y t-1 . The state process {X t , t = 1, 2, • • • }, X t ∈ R k , is the solution of the following equation, called state equation: X

  the variance-covariance matrix of the noise process ε t of the state equation (2.4.5). The description of Fisher scoring steps in matrix notation is as follows: rewrite the observations matrix with Y 0

  In this chapter, two approaches are presented for estimating the latent variable. First, the Working Extended Kalman Filtering recursions (WEKF). Second, Maximum A Posteriori via the Auxiliary Iterated Extended Kalman Particle Filtering (AIEKPF-MAP). The models parameters are estimated through the Maximum Likelihood Estimation method (MLE) via the Expectation-Maximization (EM) algorithm. The consistency and asymptotic normality for MLE's are established. The posterior distribution p[X i (t) | Y i (t)] is computed through the Bayesian approach, where we develop the Auxiliary Iterated Extended Kalman Particle filter(AIEKPF) with our model by deriving the equations of extended Kalman filter recursions.3.2 The path of the processIn the applications of life, like biometrics, economics, finance, psychometrics, public health, social sciences, the data are obtainable as longitudinal data or categorical time series. The categorical time series concept can be explained with an example. Let {U

t , t = 1, 2, • • • } be a time series that has c categories. The values of U t are 0, 1, 2, • • • , c-1. To clarify more, in health study for the children the sleep state score virtually is observed for an individual (newborn infants) the responses are (quiet sleep, indeterminate sleep, active sleep and awake), these scores can be recorded as : 1 ( quiet sleep), 2 (indeterminate sleep), 3 (Active sleep), 4 (awake).

  b) Draw the samples of latent variables X(t) by a state equation (CHARN model ) with state noise from exponential family distributions (gaussian or exponential distribution).

(c) For this individual produce 2 covariate variables u (t) : Age u 1 (t) and Sex u 2 (t), where u 1 (t) ∼ N (µ u 1 , σ 2 u 1 ) and u 2 (t) ∼ Bin(n, p). (d) For each item k, assume the β s k s are known and calculate the probabilities

  5.2.2.6 The exponential state noise with AR(1) modelIn this subsection, the simulation II outlined precedently is also Performed. We considered an AR(1) model for state equation, but with state noise from the exponential distribution ε i (t) ∼ exp(λ). In step (1) of the simulation II outline, set λ = 1, to generate the state noise. The state noise ε i (t) is centralized by the following

Table 5 .

 5 ne sais plus où j'en suis, je me sens perturbée 19-je me sens solide, posée 20-je me sens de bonne humeur, aimable steps of the first ten individuals in the ranking for every hundred individuals. We recall that the patient fatigue is varying over time around a mean value assumed to be nil in our work. Figure5.32 shows the curves of prediction and filter steps for the 1st individual to the fifth individual are identical , and the sixth individual to the tenth individual are similar. One observes the same thing with Figures 5.33-5.36. The reason is that the effects of the covariates are identical which make the values of link functions are also similar. Consequently, the predictive and filter values are almost equal for each group of individuals. The variance-covariance matrix for predictive, filter P t|t-1 , P t|t respectively have values less than 0.001. In general, on Figures 5.32-5.36, a positive value indicates that the patient is rested and a negative value indicates that he is tired.

	1: The STAI (State-Trait Anxiety Questionnaire)
	Non	Plutôt	Plutôt	Oui
		non	oui
	1-je me sens calme
	2-je me sens en sécurité
	3-je me sens tendue
	4-je me sens surmenée
	5-je me sens tranquille,bien dans ma peau
	6-je me sens bouleversée
	7-je me sens préoccupée par tout les malheurs possibles
	8-je me sens comblée
	9-je me sens effrayée
	10-je me sens à l'aise
	11-je me sens que j'ai confiance en moi
	12-je me sens nerveuse
	13-je suis affolée
	14-je me sens indécise
	15-je suis détendue
	16-je me sens satisfaite
	17-je suis préoccupée
	18-je

  .1) as p(X t+1 | Y t ) = p(X t+1 | Y t+1 )p(Y t+1 | Y t ) p(Y t+1 | X t+1 ) (A.4.2)Then equation (A.4.1) becomesp(X t | Y T ) = p(X t | Y t ) p(X t+1 | Y T )p(X t+1 | X t )p(Y t+1 | X t+1 ) p(X t+1 | Y t+1 )p(Y t+1 | Y t ) dX t+1 = p(X t | Y t ) p(Y t+1 | Y t ) p(X t+1 | X t )p(Y t+1 | X t+1 ) p(X t+1 | Y t+1 ) p(X t+1 | Y T )dX t+1 = p(X t | Y t ) p(Y t+1 | Y t ) p(X t+1 | X t )p(Y t+1 | X t+1 ) p(X t+1 | Y t+1 ) × P (X t+1 | Y T ).

	we putting the same coefficients together
	l(β) =	n i=1	q k=1	T t=1	{y 1 ik (t) log	π 1 ik (t) π c k ik (t)	+ y 2 ik (t) log	π 2 ik (t) π c k ik (t)	+ • • •
		+y c k -1 ik	(t) log	π c k -1 ik π c k ik (t) (t)	+ log π c k ik (t)}
	Where η s ik (t) = log	π s ik (t) π c k ik (t) and π c k ik (t) =	1+	1 j=1 exp[η j c k ik (t)]
										(A.4.3)

  We calculate the partial score, uses it only to about the double sum i and t:

	where									π ik =	exp[η s ik (t)] 1 + c k j=1 exp[η j ik (t)]
	∂π s ik (t) ∂η s ik (t)	=	(1 + c k j=1 exp[η j ik (t)]).(exp[η s ik (t)]) -(exp[η s ik (t)]) 2 (1 + c k j=1 exp[η j ik (t)]) 2
				=	exp[η s ik (t)]{(1 + c k j=1 exp[η j ik (t)]) -(exp[η s ik (t)])} (1 + c k j=1 exp[η j ik (t)]) 2
	1 s=1 exp[η s ik (t)] y s ik (t)η s ik (t) -log 1 + 1 + c k j=1 exp[η j ik (t)] (1 + c k j=1 exp[η j c k j=1 exp η j ik (t)]) -exp[η s ik (t)] × 1 + c k j=1 exp[η j ik (t)] ik (t) To calculate the partial score, recall the chain rule for multivariate functions to = = π s ik (t) × 1 + c k j=1 exp[η j ik (t)] 1 + c k j=1 exp[η j ik (t)] exp[η s ik (t)] -1 + c k j=1 exp[η j ik (t)]
	obtain			∂l(β) ∂β s k ik (t) × 1 -= ∂l(β) ∂η s ik (t) = π s 1 + c k ∂η s ik (t) ∂π s ik (t) exp[η s ik (t)] ∂π s ik (t) ∂η s ik (t) j=1 exp[η j ik (t)] = π s ik (t)(1 -π s ik (t))	∂η s ik (t) ∂β s k	,	(B.1.2)
	∂l(β) ∂η s ik (t) Since η s ik (t) = u i (t)β s = = D ik (X t ) n ∂ ∂η s ik (t) i=1 k + X i (t), we obtain q T c k -1 k=1 t=1 y s ik (t)η s ik (t) -log 1 + s=1 = n i=1 T t=1 y s ik (t) -1 1 + c k j=1 exp[η j × exp(η s c k j=1 ∂η s ik (t) ∂β s k = u i (t) ik (t)) exp[η j ik (t) (B.1.5) (B.1.6) ik (t)] Substitution of equations (B.1.3) ,(B.1.4), (B.1.5), (B.1.6) into (B.1.2) shows that
			= ∂l(β) k ∂β s	n i=1 =	T t=1 n i=1	y s ik (t) -T t=1 [y s ik (t) -π s 1 + c k exp(η s ik (t)) j=1 exp[η j ik (t)] ik (t)] u i (t)D ik (X t )Σ -1 ik (X t )	(B.1.7)
	= Then we using this result to find the score function of β s n T [y s ik (t) -π s ik (t)] k with EM algorithm: (B.1.3)
					i=1	t=1
	and	β	s(p) k	=		n i=1	T t=1	∂η s ik (t) ∂π s ik (t)	= Σ -1 ik (X t )	(B.1.4)
	and Σ ik (t) = cov(y s ik (t) with generic elements
	σ sm ik (t) ∼	π s ik (t)[1 -π s ik (t)], if s = m -π s ik (t)π m ik (t) if s = m
	To solve equations: ∂π s ik (t) ∂η s ik (t) we will using the quotient rule for differentiating of division two
									f g	(a) =	g(a).f (a) -f (a).g (a) [g(a)] 2

. denote the absolute value norm, if p, q are vectors, thenp -q 1 = n i=1 | p i -q i | .

many-to-one function : The parameters are a functions of full data (X, Y ). A function f which may (but does not necessarily) companion a given member of the range of f with more than one member of the domain of f . For example, trigonometric functions such as sin x are many-to-one where sinx = sin(2π + x) = sin(4π + x) = ...., .

Block matrix: In mathematics, a block matrix is a matrix that is interpreted as having been broken into sections called blocks or submatrices.

[START_REF] Cox | Theoretical statistics[END_REF], p.107

The aim to minimizing the expected value of the square of this vector E X t -xt|t

, this equivalent to minimizing the trace of the posterior covariance matrix V a t , which is minimized when its matrix derivative with respect to the gain matrix is zero.

Remerciements

. Bien sûr, atteindre ces objectifs n'aurait pas été possible sans l'aide des membres M. Antoine Lejay (Ferry, Monsieur Georges BILLANT, Madame Laurence QUIROT (responsable administrative

The latent variables are used to calculate the observations probabilities π ik (t) to generate the individuals' responses Y ik (t), after that, the latent variables are considered as unknown variables in the real data. Therefore, we apply the steps of outline to find the estimators for the model's parameters and the latent variables. Then the latent variables and their estimates are calculated via the working extended Kalman filtering are compared.

Numerical computation

The EM algorithm has the property of increasing the likelihood at each step, but the convergence is slow. As an alternative to convergence to a local maximum, we use Fisher scoring iteration method with an initial estimate θ(0) Fisher scoring iteration method is given by

where θ = (β , γ , δ ) and s( θ (p) ) the Fisher scoring of the parameters and I( θ (p) ) the Fisher information matrix of the parameters described in Chapters 3 and 4. As glm-R package the convergence occurs if

where, dev = -2 log(L), we take = 0.01

The gaussian state noise with AR(1) model

The outline of simulation experiments II is performed, where the latent variable considered follows an AR(1) model

.2.3)

To generate the individuals' responses, suppose ρ = 0.5, and R t = 1. The steps [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF][START_REF] Bartholomew | Latent variable models and factor analysis: A unified approach[END_REF][START_REF] Bartolucci | Modeling Longitudinal Data by Latent Markov Models with Application to Educational and Psychological Measurement[END_REF][START_REF] Bartolucci | Longitudinal analysis of selfreported health status by mixture latent auto-regressive models[END_REF] of outline are implemented with iterations p = 0, 1, • • • , 10, and in the initial iteration p = 0, suppose ρ (0) = 0.4, and R (0) t = 0.5. Recalling that , here η s ik (t) = β s 0k + β s 1k u 1 (t) + β s 21k u 21 (t) + β s 22k u 22 (t) (5.2.4) where u 1 (t) is the age variable and u 2 (t) is the sex variable. We divide it into two groups according to its value u 21 (t) and u 22 (t). We suppose the initial vectors of β 21q = 0.9 0.9 0.9 0.1 0.1 0.1

22q = 0.4 0.4 0.4 0.7 0.9 0.9 Two cases are considered. First, the number of individuals is small (n ≤ 15) with instant (T > 75). Second the number of individuals is large (n > 75) with instant (T ≤ 10).

• n small with T large

The graphs (5.2-5.6) show the latent variables with state noise from the gaussian distribution and their estimators (posterior mode) via the Working Extended Kalman Filter (WEKF). On the graphs, the red colour refers to the latent variables, while the blue colour refers to the posterior mode via filtering step and the green colour refers to the posterior mode via the prediction step. We took n = 1, 2, 3, 4, 10, and T = 75, 100, 200, 300. Figures 5.2-5.6 show the results. In Figure 5.2 the results are acceptable, but in Figures 5.2 (a) and (b), there are some jumps in the latent curve. Our approach has not adapted to it. In Figure 5. [START_REF] Bartholomew | Latent variable models and factor analysis: A unified approach[END_REF] The results are good except Figure 5.3 (a) show the posterior mode in two steps (prediction and filtering ) smaller than the latent variables. As seen in Figure 5.4, the results are also good and the series stable. Figure 5.5 shows that there are oscillations in the latent variables, while the posterior mod is stationary. Figures 5.6 show the results are good as the previous individuals Figures.

The gaussian state noise with CHARN(1,1)

The state space equation is described as follows

the simulation II outline is implemented. In step (1) set ρ 1 = 0.5, ρ 1 = 0.2 and R t = 1.

In steps [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF][START_REF] Bartholomew | Latent variable models and factor analysis: A unified approach[END_REF][START_REF] Bartolucci | Modeling Longitudinal Data by Latent Markov Models with Application to Educational and Psychological Measurement[END_REF][START_REF] Bartolucci | Longitudinal analysis of selfreported health status by mixture latent auto-regressive models[END_REF], set ρ

= 0.4, and R (0) t = 0.5. Figure 5.15 shows the results, where n = 10 at instants T = 75, 100, 200, 300. The working extended Kalman filter (WEKF) recursion succeeds in producing the best estimators for the state variables described by nonlinear state equation as the linear state equation with gaussian noise processes.

5.2.2.5

The gaussian state noise with CHARN(0,1)

The state space equation considered is as follows Here, the posterior mode in predictive step equal 0 due to the fact that we supposed the function F is equal 0 in the state equation. Hence the posterior mode in WEKF is a t|t-1 = F (a t|t , u i (t), γ), and it follows that a t|t-1 = 0. As seen in Figure 5.16 the results are good, this proves the efficiency of our approach.

The exponential state noise with CHARN(1,1)

The state space equation is described as follows

Here, also in step(1) of the simulation II outline suppose ρ 1 = 0.4, ρ 2 = 0.1 and λ = 1. Implement the steps (2-5) with ρ

2 = 0.3, and λ (0) = 0.6 . Figure 5.30 show the results are good as the case of a CHARN(1,1) model with gaussian state noise.

5.2.2.8

The exponential state noise with CHARN(0,1)

The state space equation is considered as

(5.2.8)

In step(1) of the simulation II outline suppose ρ 1 = 0.4, ρ 2 = 0.1 and λ = 1, the steps (2-5) are performed with ρ

2 = 0.3, and λ (0) = 0.6. Figure 5.31 display the results are good as the previous cases. As seen in Figure 5.31 the prediction equal 0 as we explained in gaussian state space.

The modeling

The observation equations is as follows

where

The link function η s ik (t) is defined to be the logit, namely

The • β s m1k , β s m2k are the parameters of the marital status classes.

• β s t1k , β s t2k , β s t3k are the parameters of the step of treatment classes.

• β s g1k , β s g2k are the parameters of the group of chemotherapy classes.

• X i (t) the fatigue of the individual i at instant t, is not observed (latent variable) and is assumed to have numerical values, and to a follow a CHARN model described by the equation:

The outline of data analysis

The following steps are performed to find the fatigue of the individuals 1. Read the data from excel. As we mentioned in preceding Chapters , the working extended Kalman filter recursions depends on F (., ., .) and H(., ., .) functions of the state equation. Then, two cases to the CHARN model are performed. First, the AR(1) model. Second, the CHARN(1,1) model.

AR(1) model

The steps of outline are implemented with AR(1) as in equation (5.3.4 ), set the initial values of parameters as: R 0 t = 0.5, ρ 0 = 0.6, and β s(0) rk is as follows

The Working Kalman filter recursions is performed. The fatigue of individuals is estimated by posterior mode, Figures 5.32 -5.36 display the filter and predictive

Appendix A Kalman Filter and Smother Recursions

A.

Derivation of the Kalman filter and smoother recursions

In this appendix, a brief derivation of the Kalman filter (Kalman (1960)) and the fixed interval smoothing algorithm [START_REF] Bierman | Numerical comparison of Kalman filter algorithms: Orbit determination case study[END_REF]) are presented.

A.1.1 One-step-Ahead prediction

At first we define

The prediction error of Y t is denoted by η t ; then from Y t = G t X t + ξ t , can be written as

Therefore, obtain

By using the facts that

where ⊕ refer to a direct sum, in the case of the normal distribution, the conditional expectation of X t given Y t is expressible by orthogonal projection as follows

Due to Proj(X t | η t ) is obtained by fall back X t on η t , from (A.1.4) and (A.1.5),

Therefore, we express

In contrast, from

we obtain

Then, the decomposition as follows

Consequently, we find that

Therefore, by putting

Here, considering that Z T produces Y T , , Obtains

Furthermore, using

Here, using

Substituting this into (A.1.18 ) yields

A.2 The derivation of the Extended Kalman filter Recursions

Let the observation process Y t : Y t ∈ R . Then the observation equation as follows:

The state process X t : X t ∈ R k , the state equation as follows:

Initially, the information certainly available is the mean, X a 0 , and the covariance, V 0 of the initial state, then the initial optimal estimate X a 0 and the covariance of error V a 0 as follows

Prediction suppose an optimal estimate is

The predictable part of X t is written as:

By using first-order Taylor expansion can be linearized the functions in the equation (A.2.2)

where

and (H.O.T) the higher order expressions are considered negligible.

(b) For a function that depends on two variables, x and y, first -order Taylor series about the point (a, b) as follows

Now we Expanding H(X t-1 , u t , δ)ε t , about two points (x a t-1 , εt ) by using (A.2.7), where

Then, obtains

where

and E[e t-1 | Y t-1 ] = 0. Consequently, the forecast value of X t yields

10)

Substituting (A.2.6) and (A.2.8) in the forecast error equation

The forecast error covariance is specified by

(A.2.12)

Filtering At time t have two value of input: the forecast value x f t with the covariance V f t and the observation Y t with conditional covariance Σ t . Our objective is to compute the best unbiased estimate. In the least squares concept, x a t of X t , one way is to suppose the estimate is a linear combination of both x f t and Y t as follows

From the unbiasedness definition

14)

Substituting (A.2.14) in (A.2.13):

As in model forecast step, expanding G(X t , u t , λ) in Taylor series about

where

Then the expectation on both sides of (A.2.16) given

where E[e f t | Y t ] = 0. Substituting (A.2.17) in (A.2.15), the state estimate (posterior mode) yields:

The error in the estimate x a t yields:

Then, the posterior covariance as follows

By using equation (A.2.12) obtains

Such as the standard Kalman filter, the posterior covariance formula for any K t can be computed by minimizing tr(V a t ) w.r.t

where the Kalman gain is:

Substituting this back in (A.2.20), obtains

A.3 The Monte Carlo Filter

This Appendix presented a derivation of the Monte Carlo filter algorithm . The readers are indicated to [START_REF] Kitagawa | Monte Carlo filter and smoother for non-Gaussian nonlinear state space models[END_REF]

A.3.1 One-step-Ahead Prediction

Suppose that m particles {f

t-1 } that can be considered as m independent realizations from the conditional density p(X t-1 | Y t-1 ) of the state X t-1 , and m particles {ε

t } from the system noise distribution p(ε) are known. Namely, assume the following

Then, the one-step-ahead predictive distribution of X t is written as

where the system noise ε t is independent of previous states and observations, the conditional distribution of the system noise yields

In contrast , since X t depends only on X t-1 and ε

where δ(.) indicates the Dirac delta function, obtains

Therefore, when realizations {ε (j) t } of p(ε t ) and {f

are given, which are m independent realizations of the distribution p(X t | Y t-1 ) , it is equivalent to approximation of the distribution p(X t | Y t-1 ) by the empirical distribution function

where I(X, a) = 0 for X < 0 and I(X, a) = 1 otherwise. In other words, that the predictive distribution p(X t | Y t-1 ) is approximated by the probability function

Then, the observation y t is given, the posterior distribution of X t is computed by

The cumulative distribution function yields

By using the formula (A.3.8), the approximation of the distribution of the filter can be obtained . We can facilely to re-approximate it by m particles f

with equal weights to perform the calculation of the prediction (A.3.1) in the next time step. This corresponds to appearing the distribution of (A.3.8) by the empirical distribution function

The m realizations {f

} can be obtained by re-sampling {p

t } with probabilities Pr(f

(A.3.10)

A.3.1.2 Smoothing

For smoothing, assume that Pr(X 1 = s

and ε (j)

t ∼ q(ε) and define (p

t|t-1 ) as follows:

Then, (p

t|t-1 ) can be considered as a realization from the joint distribution of (X 1 , • • • , X t ) when the observation Y t-1 is given. Next, given the observation y t , the distribution Pr(X 1 ≤ p

) can be updated as follows:

t|t-1 is the same as the particle p (j) t of the filter algorithm (A.3.5), the smoothing distribution p(X 1 , • • • , X t | Y t ) can be obtained by re-sampling m ndimensional vectors p

t|t-1 ) , j = 1, • • • , m with the same weights as the filter.

A.4

The derivative of the Forward-Backward smoothing

The marginal smoother can be obtained using forward-backward smoother as

Approximating the above integration by Monte Carlo integration method, one obtains

Further approximating the filtered density p(X t+1 | Y t+1 ) from the running particle filter as

we can rewrite equation (A.4.4) as

t+1|T .

The MAP estimate of the marginal smoothing density , p(X t | Y T ) can then be obtained by finding the location of its global maxima. This maximization can be performed using different optimization method. We maximize along the particles by using the equation (A.4.5), this leads to the approximate particle based MAP estimate as

Where m = 1, • • • , N and N is the number of particles in the cloud at each time step. By using equation (A.4.2), the estimator can be further simplified to

where the filtered density p(X t | Y t ) at the particle cloud {X (m) t

} N m=1 can be evaluated during the forward filtering step as

The first-order derivative of parameter β for the multinomial logistic regression model

The likelihood is a product of the multinomial probabilities

the log-likelihood is given by

By substituting

To calculate the second partial score, recall the chain rule for multivariate functions to obtain

and we can find the derivation for the probability to β s k , where

Substitution of equations (B.1.10) ,(B.1.11) into (B.1.9) shows that

Appendix C

The exponential family distributions C.1 Members of the exponential family distributions

Exponential distribution

A gamma distribution with shape parameter α = 1 and scale parameter λ is an exponential (λ) distribution.

Laplace distribution

The mean µ is known

The log likelihood for an individual based on a set of independent variables

where

,and φ i = (φ 1 (1), • • • , φ 1 (T )) . There are advantages in using as the goodness-of-fit criterion, not the log likelihood l(X i , µ i , φ i ) but a particular linear function, namely

where νi (t) = ν i [X i (t)] and νi (t) = ν[μ i (t)], which they called the scaled deviance. we present Dev * (X i , µ i ) in Appendix (F) for the exponential family distributions. Note that, for the exponential family distributions considered here, l(X i , φ i , X i ) is the maximum likelihood achievable for an exact fit in which the fitted values are equal to latent variables, they consider maximizing l(X i , µ i , φ i ) is equivalent to minimizing Dev * (X i , µ i ) with respect to µ i . Deviance, as a measure of goodness of fit, is defined as following Dev(X i , µ i ) = φ i .Dev * (X i , µ i ) (C.2.1)

where Dev(X i , µ i ) is known as the deviance for the current model and is a function of the data only. We can estimate φ i by the mean scaled Deviance

where N = (n × T ) is the sample size and r is the total number of unknown parameters. Dev * i = Dev i /φ i is approximately χ 2 (N -r) distributed with expectation (N -r). Thus φd is approximately unbiased.

Example 1: we have for an individual

C.2.2 Scaled Pearson estimation

The Pearson chi-squared statistics which take the form where for an individual Z i (t) 1 has an exponential family distribution