From 91d43c8da613fa83343e31d56bf8483f4f27384e Mon Sep 17 00:00:00 2001 From: "Bouchra R. Nasri" Date: Sat, 8 Jul 2023 09:50:10 +0000 Subject: [PATCH] version 1.1.1 --- DESCRIPTION | 8 +- MD5 | 20 +- NAMESPACE | 4 +- R/EstHMM1d.R | 136 +++----- R/EstRegime.R | 2 +- R/GofHMM1d.R | 14 +- R/em.step.R | 117 ------- man/EstHMM1d.Rd | 108 +++---- man/EstRegime.Rd | 2 +- man/GofHMM1d.Rd | 6 +- man/em.step.Rd | 38 --- src/HMM_est_1d_new.c | 609 ++++++++++++++++++++++++++++++++++++ src/registerDynamicSymbol.c | 19 ++ 13 files changed, 748 insertions(+), 335 deletions(-) delete mode 100644 R/em.step.R delete mode 100644 man/em.step.Rd create mode 100644 src/HMM_est_1d_new.c create mode 100644 src/registerDynamicSymbol.c diff --git a/DESCRIPTION b/DESCRIPTION index f245b06..70cd72c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,17 +1,17 @@ Package: GaussianHMM1d Title: Inference, Goodness-of-Fit and Forecast for Univariate Gaussian Hidden Markov Models -Version: 1.1.0 +Version: 1.1.1 Authors@R: c(person(c("Bouchra","R."), "Nasri", role = c("aut", "cre","cph"), email = "bouchra.nasri@umontreal.ca"),person(c("Bruno","N"), "Remillard", role = c("aut", "ctb","cph"), email = "bruno.remillard@hec.ca")) Description: Inference, goodness-of-fit test, and prediction densities and intervals for univariate Gaussian Hidden Markov Models (HMM). The goodness-of-fit is based on a Cramer-von Mises statistic and uses parametric bootstrap to estimate the p-value. The description of the methodology is taken from Chapter 10.2 of Remillard (2013) . Depends: R (>= 3.5.0), doParallel, parallel, foreach, stats License: GPL (>= 2) Encoding: UTF-8 RoxygenNote: 7.2.3 -NeedsCompilation: no -Packaged: 2023-06-23 18:43:54 UTC; Utilisateur +NeedsCompilation: yes +Packaged: 2023-06-26 22:37:30 UTC; Utilisateur Author: Bouchra R. Nasri [aut, cre, cph], Bruno N Remillard [aut, ctb, cph] Maintainer: Bouchra R. Nasri Repository: CRAN -Date/Publication: 2023-06-24 18:50:08 UTC +Date/Publication: 2023-07-08 10:50:10 UTC diff --git a/MD5 b/MD5 index 0495de8..c71465c 100644 --- a/MD5 +++ b/MD5 @@ -1,30 +1,30 @@ -916fbf849f9503d06d7287f1ff0f9abe *DESCRIPTION -592d12cff367047da859982091bc14b8 *NAMESPACE -2d39a1a6d4cfc5847d6035e41e63731b *R/EstHMM1d.R -03d43ff3b4b23104672986e2b2043aad *R/EstRegime.R +7c7c296c8bd3ba9dd80ff858d2fac811 *DESCRIPTION +27424a713c86437044a3eb29c591d127 *NAMESPACE +4ecada836e99f49bdd391bae46607f4b *R/EstHMM1d.R +69b7808db30393bf49798eb1d18d77a6 *R/EstRegime.R 1c83ac676483e28928098ca1845ddfb9 *R/ForecastHMMPdf.R 1dae3f8b23b474ca1a5b132ce25a7fd7 *R/ForecastHMMeta.R eaf9446716e911bc3fc9628d12f7ecba *R/GaussianMixtureCdf.R 084fb9c24fba98aa7819866ababdf63a *R/GaussianMixtureInv.R 1970f755a302644c1cbfbcdb6e3babf1 *R/GaussianMixturePdf.R -b58350964fcfbca70b34ff26cbf5cb0b *R/GofHMM1d.R +199c9f70228db2c4b04335d5b153162e *R/GofHMM1d.R 2b5b54455d7f3f903929b7636923d0c1 *R/Sim.HMM.Gaussian.1d.R 7ec2d8d5c39feae14f4b61dd3fc8b28a *R/Sim.Markov.Chain.R f63971bccfb643081a1bda8fa3403fa7 *R/SimHMMGaussianInv.R ffd7e50f029f133a6583725bbb37e665 *R/Sn.R 360a390cc4195e6fd8b3731c4700d296 *R/bootstrapfun.R -cea4c71f8d59658eca58f8dc3976a319 *R/em.step.R -cb1e3da8c25f58b00b89c30393c077fb *man/EstHMM1d.Rd -16299a226ca24e44de4924181eb8ccda *man/EstRegime.Rd +30bb603174892fc398ef57181b4ae379 *man/EstHMM1d.Rd +9146d9d5aabefec795992ba3f1a0c6f8 *man/EstRegime.Rd ed94c57f062999e12ab6ebe094306ead *man/ForecastHMMPdf.Rd a3c166de6107ba31e3ca35d63b451b6a *man/ForecastHMMeta.Rd bc5c2556ec071a4da5a02e75b65adc8d *man/GaussianMixtureCdf.Rd 10ed6a47dbd53800f3cdc6559b4d07cf *man/GaussianMixtureInv.Rd 3841bf18cc9855a94cf0c258490bab12 *man/GaussianMixturePdf.Rd -f489e418035168b8fe0f1b7d53c78be4 *man/GofHMM1d.Rd +9fc3a45898a9159aa21e31aad1cabe44 *man/GofHMM1d.Rd b6d82fed0e1f6a86189b25d2aa440ddc *man/Sim.HMM.Gaussian.1d.Rd 8396f64b05881fb05cc8d1b5054293cd *man/Sim.Markov.Chain.Rd 627f6f66e66c8668b4495f6fefe86dfc *man/SimHMMGaussianInv.Rd ea426af72e455f6866dc026d085ca69d *man/Sn.Rd 656e498dbac4d477ece81d33f67a19d1 *man/bootstrapfun.Rd -6ba231d877bf7dafc8b6980d89022f40 *man/em.step.Rd +505d24b5ccc9054a24b691e73f8542c3 *src/HMM_est_1d_new.c +a4c9f1d9379b0948c3b92b70cf52fe10 *src/registerDynamicSymbol.c diff --git a/NAMESPACE b/NAMESPACE index 3c2dc4e..b87c8a9 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,4 @@ -# Generated by roxygen2: do not edit by hand - +useDynLib(GaussianHMM1d,.registration = TRUE) export(EstHMM1d) export(EstRegime) export(ForecastHMMPdf) @@ -13,7 +12,6 @@ export(Sim.Markov.Chain) export(SimHMMGaussianInv) export(Sn) export(bootstrapfun) -export(em.step) import(foreach) importFrom(doParallel,registerDoParallel) importFrom(graphics,plot) diff --git a/R/EstHMM1d.R b/R/EstHMM1d.R index bd09d3a..a887a40 100644 --- a/R/EstHMM1d.R +++ b/R/EstHMM1d.R @@ -9,7 +9,7 @@ #'@param y (nx1) vector of data #'@param reg number of regimes #'@param max_iter maximum number of iterations of the EM algorithm; suggestion 10 000 -#'@param prec precision (stopping criteria); suggestion 0.0001. +#'@param eps precision (stopping criteria); suggestion 0.0001. #' #'@author Bouchra R Nasri and Bruno N Rémillard, January 31, 2019 #' @@ -19,7 +19,7 @@ #'@return \item{eta}{(n x reg) probabilities of being in regime k at time t given observations up to time t} #'@return \item{lambda}{(n x reg) probabilities of being in regime k at time t given all observations} #'@return \item{cvm}{Cramér-von Mises statistic for the goodness-of-fit test} -#'@return \item{W}{Pseudo-observations that should be uniformly distributed under the null hypothesis of a Gaussian HMM} +#'@return \item{U}{Pseudo-observations that should be uniformly distributed under the null hypothesis of a Gaussian HMM} #'@return \item{LL}{Log-likelihood} #' #'@references Chapter 10.2 of B. Rémillard (2013). Statistical Methods for Financial Engineering, @@ -27,103 +27,45 @@ #' #'@examples Q <- matrix(c(0.8, 0.3, 0.2, 0.7),2,2); mu <- c(-0.3 ,0.7) ; sigma <- c(0.15,0.05) #'data <- Sim.HMM.Gaussian.1d(mu,sigma,Q,eta0=1,100)$x -#'est <- EstHMM1d(data, 2, max_iter=10000, prec=0.0001) +#'est <- EstHMM1d(data, 2, max_iter=10000, eps=0.0001) #' #'@importFrom stats dnorm pnorm sd #'@export #' -EstHMM1d <- function(y,reg,max_iter=10000,prec=0.0001){ -n <- length(y) -r <- reg - - -if(r==1){ - mu <- mean(y) - sigma <- sd(y) - Q <- 1 - eta <- rep(1,n) - lambda <- eta - - Z = (y-mu)/sigma - L = dnorm(Z)/sigma - W = pnorm(Z); - - LL <- sum(log(L)) - cvm <- Sn(W) - nu <- 1} else{ - -n0 <- floor(n/r) - -x <- matrix(0,n0,r) - -ind0 <- (1:n0) - -for(j in 1:r){ - ind <- (j-1)*n0 + ind0 - x[,j] <- y[ind] -} - -mu0 <- apply(x,2,mean) -sigma0 <- apply(x,2,sd) -Q0 <- matrix(1,r,r)/r - -for(k in 1:100){ - curr.step <- em.step(y,mu0,sigma0,Q0) - #eta <- curr.step$eta - mu0 <- curr.step$mu - Q0 <- curr.step$Q - sigma0 <- curr.step$sigma -} - - -rprec= r*prec; - - -for(k in 1:max_iter){ - curr.step <- em.step(y,mu0,sigma0,Q0) - eta <- curr.step$eta - - sum1 <- sum(abs(mu0)) - sum2 <- sum(abs(curr.step$mu-mu0)) - - - if ( sum2 < (sum1*rprec) ){ - break - } - mu0 <- curr.step$mu - Q0 <- curr.step$Q - sigma0 <- curr.step$sigma - - -} - - -mu <- curr.step$mu -Q <- curr.step$Q -sigma <- curr.step$sigma -LL<-curr.step$LL -lambda<-curr.step$lambda -nu=curr.step$nu - -Z=matrix(0,n,r) - -for (i in 1:r){Z[,i]=(y-mu[i])/sigma[i]} - -U = pnorm(Z); - -eta00 = matrix(1,1,r)/r; - -w00 = rbind(eta00, eta) %*% Q; - -dd=dim(w00)[1] - -w = w00[-dd,]; - -W = as.vector(rowSums( w * U)); - -cvm = Sn(W); -} - -return(list(mu=mu,sigma=sigma,Q=Q,eta=eta,nu=nu,lambda=lambda,LL=LL,cvm=cvm, W=W)) +EstHMM1d <- function(y,reg,max_iter=10000,eps=0.0001){ + + + + n = length(y) + + out0=.C("est_hmm_1d_new", + as.double(y), + as.integer(reg), + as.integer(n), + as.double(eps), + as.double(max_iter), + mu=double(reg), + sigma = double(reg), + nu = double(reg), + Qvec = double(reg*reg), + etavec = double(n*reg), + lambdavec = double(n*reg), + U = double(n), + cvm = double(1), + L = double(n)) + +eta = matrix(out0$etavec,ncol=reg) +lambda = matrix(out0$lambdavec,ncol=reg) +Q= matrix(out0$Qvec,ncol=reg) + +mu <- out0$mu +sigma <- out0$sigma +nu <- out0$nu +U <- out0$U +cvm = out0$cvm +LL = sum(log(out0$L)) +if(is.na(sum(U))){out=NULL; +cat("Singular values. Number of regimes too large.")}else{ +out=list(mu=mu,sigma=sigma,Q=Q,eta=eta,lambda=lambda,nu=nu,U=U,cvm=cvm,LL=LL)} +return(out) } - diff --git a/R/EstRegime.R b/R/EstRegime.R index 512485e..17688af 100644 --- a/R/EstRegime.R +++ b/R/EstRegime.R @@ -22,7 +22,7 @@ #'@examples Q <- matrix(c(0.8, 0.3, 0.2, 0.7),2,2); mu <- c(-0.3 ,0.7) ; sigma <- c(0.15,0.05); #'data <- Sim.HMM.Gaussian.1d(mu,sigma,Q,eta0=1,100)$x #'t=c(1:100); -#'est <- EstHMM1d(data, 2, max_iter=10000, prec=0.0001) +#'est <- EstHMM1d(data, 2) #'EstRegime(t,data,est$lambda, est$eta) #' #'@importFrom graphics plot diff --git a/R/GofHMM1d.R b/R/GofHMM1d.R index 058b302..d60242c 100644 --- a/R/GofHMM1d.R +++ b/R/GofHMM1d.R @@ -7,7 +7,7 @@ #'@param y (n x 1) data vector #'@param reg number of regimes #'@param max_iter maxmimum number of iterations of the EM algorithm; suggestion 10 000 -#'@param prec precision (stopping criteria); suggestion 0.0001 +#'@param eps eps (stopping criteria); suggestion 0.0001 #'@param n_sample number of bootstrap samples; suggestion 1000 #'@param n_cores number of cores to use in the parallel computing #' @@ -31,7 +31,7 @@ #'\donttest{ #'Q <- matrix(c(0.8, 0.3, 0.2, 0.7),2,2); mu <- c(-0.3 ,0.7) ; sigma <- c(0.15,0.05) #'data <- Sim.HMM.Gaussian.1d(mu,sigma,Q,eta0=1,100)$x -#'gof <- GofHMM1d(data, 2, max_iter=10000, prec=0.0001, n_sample=100,n_cores=2) +#'gof <- GofHMM1d(data, 2, max_iter=10000, eps=0.0001, n_sample=100,n_cores=2) #'} #'@importFrom doParallel registerDoParallel #'@importFrom parallel makeCluster stopCluster @@ -42,14 +42,14 @@ #' -GofHMM1d <-function(y, reg, max_iter=10000, prec=0.0001, n_sample=1000,n_cores){ +GofHMM1d <-function(y, reg, max_iter=10000, eps=0.0001, n_sample=1000,n_cores){ cl <- makeCluster(n_cores) registerDoParallel(cl) #n_cores = detectCores()-2; - out0 = EstHMM1d(y,reg,max_iter, prec); + out0 = EstHMM1d(y,reg,max_iter, eps); print("End of estimation"); mu = out0$mu; @@ -63,10 +63,10 @@ GofHMM1d <-function(y, reg, max_iter=10000, prec=0.0001, n_sample=1000,n_cores){ LL = out0$LL; n= length(y); - fun = c('Sim.Markov.Chain','Sim.HMM.Gaussian.1d','EstHMM1d','Sn','bootstrapfun','em.step') + fun = c('Sim.Markov.Chain','Sim.HMM.Gaussian.1d','EstHMM1d','Sn','bootstrapfun') - # result <- foreach(i=1:n_sample, .packages='GaussianHMM1d') %dopar% bootstrapfun(mu,sigma,Q,max_iter,prec,n) - result <- foreach(i=1:n_sample, .export=fun) %dopar% bootstrapfun(mu,sigma,Q,max_iter,prec,n) + # result <- foreach(i=1:n_sample, .packages='GaussianHMM1d') %dopar% bootstrapfun(mu,sigma,Q,max_iter,eps,n) + result <- foreach(i=1:n_sample, .export=fun, .packages="GaussianHMM1d") %dopar% bootstrapfun(mu,sigma,Q,max_iter,eps,n) cvm_sim = rep(0,n_sample) for (i in 1:n_sample){ diff --git a/R/em.step.R b/R/em.step.R deleted file mode 100644 index 54cf1c9..0000000 --- a/R/em.step.R +++ /dev/null @@ -1,117 +0,0 @@ -#'@title Function to perform the E-M steps for the estimation of the paramters -#' -#'@description This function perform the E-M steps for the estimation of the parameters -#'of a univariate Gaussian HMM. -#' -#'@param y points at which the density function is comptuted (mx1); -#'@param mu vector of means for each regime (r x 1); -#'@param sigma vector of standard deviations for each regime (r x 1); -#'@param Q transition probality matrix (r x r); -#' -#'@author Bouchra R Nasri and Bruno N Rémillard, January 31, 2019 -#' -#'@return \item{f}{values of the density function at time n+k} -#'@return \item{w}{weights of the mixture} -#' -#'@references Chapter 10.2 of B. Rémillard (2013). Statistical Methods for Financial Engineering, -#'Chapman and Hall/CRC Financial Mathematics Series, Taylor & Francis. -#' -#'@examples mu <- c(-0.3 ,0.7) ; sigma <- c(0.15,0.05); Q <- matrix(c(0.8, 0.3, 0.2, 0.7),2,2) ; -#'data <- Sim.HMM.Gaussian.1d(mu,sigma,Q,eta0=1,100)$x -#'out <- em.step(data,mu,sigma,Q) -#'@importFrom stats dnorm -#' -#'@export -#' -em.step <- function(y,mu,sigma,Q){ - - n <- length(y) - r <- length(mu) - gamma <- matrix(0,n,r) - gammabar <- matrix(0,n,r) - lambda <- matrix(0,n,r) - f <- matrix(0,n,r) - Lambda <- array(0, dim=c(r,r,n)) - ## - for(j in 1:r){ - z <- (y-mu[j])/sigma[j] - f[,j] <- dnorm(z)/sigma[j] - } - - ## - gammabar[n,]<-1/r - - for(k in 1: (n-1)){ - i <- n-k - j=i+1 - v = ( gammabar[j,] * f[j,] ) %*% t(Q) - gammabar[i,] = v/sum(v) - #gammabar[i,] <- ( gammabar[i+1,] * f[i+1,] ) * t(Q) - } - - ## - L=matrix(0,1,n) - eta=matrix(0,n,r) - eta0 <- matrix(1,1,r)/r - - v <- ( eta0 %*% Q) * f[1,] - L[1]=sum(v) - eta[1,] <- v/sum(v) - - for(i in 2:n){ - v <- ( eta[i-1,] %*% Q) * f[i,] - L[i]=sum(v) - eta[i,] <- v/L[i] - } - - LL=sum(log(L)) - ## - v <- eta * gammabar - sv0=rowSums(v) - for(j in 1:(r)){lambda[,j]=v[,j]/sv0} - - #sv <- matrix(1,1,r) %*% apply(v,1,sum) - #lambda <- v / sv - - ## - #Lambda[,,n] <- matrix(1,r,1) %*% lambda[n,] * Q - - #for (j in (1:r)){ Lambda[j,,n]=lambda[n,j] %*% Q[j,] } - for (j in (1:r)){ Lambda[j,,n]=lambda[n,j]*Q[j,] } - - gf=gammabar*f - - - for(i in 1:(n-1)){ - M <- Q*( matrix(1,r,1) *(eta[i,]))%*%( matrix( 1,1,r) * gf[i+1,] ) - c <- sum(apply(M,2,sum)) - Lambda[,,i] <- M/c - } - - ## - nu <- apply(lambda,2,mean) - - munew=matrix(0,1,r) - sigmanew=matrix(0,1,r) - Qnew=matrix(0,r,r) - - - for (j in 1:r) - { - w=lambda[,j]/nu[j]/n - #yy <- matrix(1,1,r) %*% y - munew[j]<- sum( y * w ) - sigmanew[j] <- sqrt( sum( y^2 * w) - munew[j]^2 ) - Qnew[j,] <- apply((Lambda[j,,]),1,mean) / nu[j] - } - - - - # w <- lambda / matrix(1,n,1)%*%nu /n - #yy <- matrix(1,1,r) %*% y - #munew <- sum( yy * w ) - #sigmanew <- sqrt( apply( yy^2 * w,2,sum) - munew^2 ) - #Qnew <- apply(Lambda,3,mean) / matrix(1,1,r)%*% t(nu) - - return(list(nu=nu, munew=munew, sigmanew=sigmanew, Qnew=Qnew, eta=eta, gammabar=gammabar, lambda=lambda, Lambda=Lambda, LL=LL)) -} diff --git a/man/EstHMM1d.Rd b/man/EstHMM1d.Rd index 27c23db..be492a3 100644 --- a/man/EstHMM1d.Rd +++ b/man/EstHMM1d.Rd @@ -1,54 +1,54 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/EstHMM1d.R -\name{EstHMM1d} -\alias{EstHMM1d} -\title{Estimation of a univariate Gaussian Hidden Markov Model (HMM)} -\usage{ -EstHMM1d(y, reg, max_iter = 10000, prec = 1e-04) -} -\arguments{ -\item{y}{(nx1) vector of data} - -\item{reg}{number of regimes} - -\item{max_iter}{maximum number of iterations of the EM algorithm; suggestion 10 000} - -\item{prec}{precision (stopping criteria); suggestion 0.0001.} -} -\value{ -\item{mu}{estimated mean for each regime} - -\item{sigma}{stimated standard deviation for each regime} - -\item{Q}{(reg x reg) estimated transition matrix} - -\item{eta}{(n x reg) probabilities of being in regime k at time t given observations up to time t} - -\item{lambda}{(n x reg) probabilities of being in regime k at time t given all observations} - -\item{cvm}{Cramér-von Mises statistic for the goodness-of-fit test} - -\item{W}{Pseudo-observations that should be uniformly distributed under the null hypothesis of a Gaussian HMM} - -\item{LL}{Log-likelihood} -} -\description{ -This function estimates parameters (mu, sigma, Q) of a univariate Hidden Markov Model. -It computes also the probability of being in each regime, given the past observations (eta) -and the whole series (lambda). The conditional distribution given past observations is applied to -obtains pseudo-observations W that should be uniformly distributed under the null hypothesis. -A Cramér-von Mises test statistic is then computed. -} -\examples{ -Q <- matrix(c(0.8, 0.3, 0.2, 0.7),2,2); mu <- c(-0.3 ,0.7) ; sigma <- c(0.15,0.05) -data <- Sim.HMM.Gaussian.1d(mu,sigma,Q,eta0=1,100)$x -est <- EstHMM1d(data, 2, max_iter=10000, prec=0.0001) - -} -\references{ -Chapter 10.2 of B. Rémillard (2013). Statistical Methods for Financial Engineering, -Chapman and Hall/CRC Financial Mathematics Series, Taylor & Francis. -} -\author{ -Bouchra R Nasri and Bruno N Rémillard, January 31, 2019 -} +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/EstHMM1d.R +\name{EstHMM1d} +\alias{EstHMM1d} +\title{Estimation of a univariate Gaussian Hidden Markov Model (HMM)} +\usage{ +EstHMM1d(y, reg, max_iter = 10000, eps = 1e-04) +} +\arguments{ +\item{y}{(nx1) vector of data} + +\item{reg}{number of regimes} + +\item{max_iter}{maximum number of iterations of the EM algorithm; suggestion 10 000} + +\item{eps}{precision (stopping criteria); suggestion 0.0001.} +} +\value{ +\item{mu}{estimated mean for each regime} + +\item{sigma}{stimated standard deviation for each regime} + +\item{Q}{(reg x reg) estimated transition matrix} + +\item{eta}{(n x reg) probabilities of being in regime k at time t given observations up to time t} + +\item{lambda}{(n x reg) probabilities of being in regime k at time t given all observations} + +\item{cvm}{Cramér-von Mises statistic for the goodness-of-fit test} + +\item{U}{Pseudo-observations that should be uniformly distributed under the null hypothesis of a Gaussian HMM} + +\item{LL}{Log-likelihood} +} +\description{ +This function estimates parameters (mu, sigma, Q) of a univariate Hidden Markov Model. +It computes also the probability of being in each regime, given the past observations (eta) +and the whole series (lambda). The conditional distribution given past observations is applied to +obtains pseudo-observations W that should be uniformly distributed under the null hypothesis. +A Cramér-von Mises test statistic is then computed. +} +\examples{ +Q <- matrix(c(0.8, 0.3, 0.2, 0.7),2,2); mu <- c(-0.3 ,0.7) ; sigma <- c(0.15,0.05) +data <- Sim.HMM.Gaussian.1d(mu,sigma,Q,eta0=1,100)$x +est <- EstHMM1d(data, 2, max_iter=10000, eps=0.0001) + +} +\references{ +Chapter 10.2 of B. Rémillard (2013). Statistical Methods for Financial Engineering, +Chapman and Hall/CRC Financial Mathematics Series, Taylor & Francis. +} +\author{ +Bouchra R Nasri and Bruno N Rémillard, January 31, 2019 +} diff --git a/man/EstRegime.Rd b/man/EstRegime.Rd index a0855b8..1a6d2c7 100644 --- a/man/EstRegime.Rd +++ b/man/EstRegime.Rd @@ -37,7 +37,7 @@ and probabilities of being in regime k at time t given observations up to time t Q <- matrix(c(0.8, 0.3, 0.2, 0.7),2,2); mu <- c(-0.3 ,0.7) ; sigma <- c(0.15,0.05); data <- Sim.HMM.Gaussian.1d(mu,sigma,Q,eta0=1,100)$x t=c(1:100); -est <- EstHMM1d(data, 2, max_iter=10000, prec=0.0001) +est <- EstHMM1d(data, 2) EstRegime(t,data,est$lambda, est$eta) } diff --git a/man/GofHMM1d.Rd b/man/GofHMM1d.Rd index c4ab6e4..aa04d31 100644 --- a/man/GofHMM1d.Rd +++ b/man/GofHMM1d.Rd @@ -4,7 +4,7 @@ \alias{GofHMM1d} \title{Goodness-of-fit test of a univariate Gaussian Hidden Markov Model} \usage{ -GofHMM1d(y, reg, max_iter = 10000, prec = 1e-04, n_sample = 1000, n_cores) +GofHMM1d(y, reg, max_iter = 10000, eps = 1e-04, n_sample = 1000, n_cores) } \arguments{ \item{y}{(n x 1) data vector} @@ -13,7 +13,7 @@ GofHMM1d(y, reg, max_iter = 10000, prec = 1e-04, n_sample = 1000, n_cores) \item{max_iter}{maxmimum number of iterations of the EM algorithm; suggestion 10 000} -\item{prec}{precision (stopping criteria); suggestion 0.0001} +\item{eps}{eps (stopping criteria); suggestion 0.0001} \item{n_sample}{number of bootstrap samples; suggestion 1000} @@ -46,7 +46,7 @@ using parametric bootstrap. \donttest{ Q <- matrix(c(0.8, 0.3, 0.2, 0.7),2,2); mu <- c(-0.3 ,0.7) ; sigma <- c(0.15,0.05) data <- Sim.HMM.Gaussian.1d(mu,sigma,Q,eta0=1,100)$x -gof <- GofHMM1d(data, 2, max_iter=10000, prec=0.0001, n_sample=100,n_cores=2) +gof <- GofHMM1d(data, 2, max_iter=10000, eps=0.0001, n_sample=100,n_cores=2) } } \references{ diff --git a/man/em.step.Rd b/man/em.step.Rd deleted file mode 100644 index 460d84f..0000000 --- a/man/em.step.Rd +++ /dev/null @@ -1,38 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/em.step.R -\name{em.step} -\alias{em.step} -\title{Function to perform the E-M steps for the estimation of the paramters} -\usage{ -em.step(y, mu, sigma, Q) -} -\arguments{ -\item{y}{points at which the density function is comptuted (mx1);} - -\item{mu}{vector of means for each regime (r x 1);} - -\item{sigma}{vector of standard deviations for each regime (r x 1);} - -\item{Q}{transition probality matrix (r x r);} -} -\value{ -\item{f}{values of the density function at time n+k} - -\item{w}{weights of the mixture} -} -\description{ -This function perform the E-M steps for the estimation of the parameters -of a univariate Gaussian HMM. -} -\examples{ -mu <- c(-0.3 ,0.7) ; sigma <- c(0.15,0.05); Q <- matrix(c(0.8, 0.3, 0.2, 0.7),2,2) ; -data <- Sim.HMM.Gaussian.1d(mu,sigma,Q,eta0=1,100)$x -out <- em.step(data,mu,sigma,Q) -} -\references{ -Chapter 10.2 of B. Rémillard (2013). Statistical Methods for Financial Engineering, -Chapman and Hall/CRC Financial Mathematics Series, Taylor & Francis. -} -\author{ -Bouchra R Nasri and Bruno N Rémillard, January 31, 2019 -} diff --git a/src/HMM_est_1d_new.c b/src/HMM_est_1d_new.c new file mode 100644 index 0000000..9d4265b --- /dev/null +++ b/src/HMM_est_1d_new.c @@ -0,0 +1,609 @@ +#include +#include +#include + +#define cte 0.3989423 + + double maxi(double a, double b) + { + if (a > b) + return a; + else + return b; + } + +double sumx( double *x, int n) + { + int i; + double sum = 0.0; + + for(i=0;i // for NULL +#include + +/* FIXME: + Check these declarations against the C/Fortran source code. +*/ + + /* .C calls */ +extern void est_hmm_1d_new(void *, void *, void *, void *,void *, void *, void *, void *, void *, void *, void *, void *, void *, void *); +static const R_CMethodDef CEntries[] = { + {"est_hmm_1d_new", (DL_FUNC) &est_hmm_1d_new, 14}, + {NULL, NULL, 0} +}; + +void R_init_GaussianHMM1d(DllInfo *dll) +{ + R_registerRoutines(dll, CEntries, NULL, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +}