version 1.1.1

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) <doi:10.1201/b14285>.
 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 <bouchra.nasri@umontreal.ca>
 Repository: CRAN
-Date/Publication: 2023-06-24 18:50:08 UTC
+Date/Publication: 2023-07-08 10:50:10 UTC

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

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)

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,111 +19,53 @@
 #'@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,
 #'Chapman and Hall/CRC Financial Mathematics Series, Taylor & Francis.
 #'
 #'@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)
 }
-

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

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){