#==================================================================================================#
# Date: 14/04/2020
# Description: Function for estimate number of regimes when unknown. Default method is
# Metropolized Carlin and Chib, when NAIC is TRUE number of regimes is choose only with this result.
# Function:
#==================================================================================================#
mtarnumreg = function(ini_obj, level = 0.95, burn_m = NULL, niter_m = 1000,
iterprev = 500, chain_m = FALSE, list_m = FALSE, NAIC = FALSE,
ordersprev = list(maxpj = 2,maxqj = 0,maxdj = 0), parallel = TRUE){
compiler::enableJIT(3)
if (!is.logical(chain_m)) {stop('chain_m must be a logical object')}
if (!is.logical(parallel)) {stop('parallel must be a logical object')}
if (!is.logical(NAIC)) {stop('NAIC must be a logical object')}
# Checking
if (!inherits(ini_obj, 'regime_inipars')) {
stop('ini_obj must be a regime_inipars object')
}
if (is.null(ini_obj$tsregime_obj$Zt)) {
stop('Threshold process must be enter in ini_obj for evaluate l0_max number of regimes')}
if (is.null(ini_obj$l0_min)) {
message('l0_min NULL default 2 \n')
l0_min = 2
}else{
l0_min = ini_obj$l0_min
}
if (is.null(ini_obj$l0_max)) {
message('l0_max NULL default 3 \n')
l0 = 3
}else{
l0 = ini_obj$l0_max
}
if (is.null(ini_obj$method)) {
message('method NULL default KUO \n')
method = 'KUO'
}else{
method = ini_obj$method
}
if(is.null(ini_obj$tsregime_obj$Xt) & ordersprev$maxqj > 0 ){
warning('Xt does not exist, ordersprev $ maxqj = 0')
ordersprev$maxqj = 0
}
# data
Yt = ini_obj$tsregime_obj$Yt
Ut = cbind(ini_obj$tsregime_obj$Zt,ini_obj$tsregime_obj$Xt)
k = ini_obj$tsregime_obj$k
N = ini_obj$tsregime_obj$N
nu = ini_obj$tsregime_obj$nu
if (is.null(nu)) {nu = 0}
maxpj = ifelse(is.null(ordersprev$maxpj),2,ordersprev$maxpj)
maxqj = ifelse(is.null(ordersprev$maxqj),0,ordersprev$maxqj)
maxdj = ifelse(is.null(ordersprev$maxdj),0,ordersprev$maxdj)
# Code
burn_m = ifelse(is.null(burn_m),round(0.3*niter_m),burn_m)
Yt = t(Yt)
Zt = Ut[,1]
if (nu == 0) {
Xt = matrix(0,ncol = N,nrow = 1)
}else{
Xt = t(ini_obj$tsregime_obj$Xt)
}
rini = ini_obj$init$r
a = ifelse(is.null(rini$za),min(Zt),stats::quantile(Zt,probs = rini$za))
b = ifelse(is.null(rini$zb),max(Zt),stats::quantile(Zt,probs = rini$zb))
### FUNCTIONS
dmunif = function(r, a, b){
l = length(r) + 1
names(a) = names(b) = NULL
volume = ((b - a)^{l - 1})/(factorial(l - 1))
for (i in 1:{l - 1}) {
if (r[i] >= a & r[i] <= b) {
if (l <= 2) {prob = 1}
if (l > 2) {
prob = 1
for (j in 1:{l - 2}) {
if (r[j] < r[j + 1]) {prob = prob*1}else{prob = prob*0}
}
}
}else{prob = 0}
}
rj = matrix(nrow = 2,ncol = l)
rj[,1] = c(-Inf,r[1])
rj[,l] = c(rev(r)[1],Inf)
if (l > 2) {
for (i2 in 2:{l - 1}) {rj[,i2] = c(r[i2 - 1],r[i2])}
}
Ind = c()
for (j in 1:l) {
for (w in 1:N) {
if (Zt[w] > rj[1,j] & Zt[w] <= rj[2,j]) {
Ind[w] = j
}
}
}
Nrg = c()
for (lj in 1:l) {
Nrg[lj] = length(Ind[Ind == lj])
}
if (sum(Nrg/sum(Nrg) > 0.2) == l) {prob = 1*prob}else{prob = 0*prob}
return(prob/volume)
}
dmunif = compiler::cmpfun(dmunif)
rdunif = function(m,l0,...){
sec = l0_min:l0
sec = sec[sec != m]
if (length(sec) == 1) {
muestra = sec
}else{
muestra = sample(x = sec, size = 1)
}
return(muestra)
}
rdunif = compiler::cmpfun(rdunif)
### Function to create lists
lists = function(l, r, pjmax, qjmax, djmax, etam, ...){
rj = matrix(nrow = 2,ncol = l)
rj[,1] = c(-Inf,r[1])
rj[,l] = c(rev(r)[1],Inf)
if (l > 2) {for (j2 in 2:{l - 1}) {rj[,j2] = c(r[j2 - 1],r[j2])}}
# indicator variable for the regime
Ind = vector(mode = 'numeric',length = N)
for (j in 1:l) {
Ind[Zt > rj[1,j] & Zt <= rj[2,j]] = j
}
Nrg = vector(mode = 'numeric')
listaXj = listaYj = list()
length(listaXj) = length(listaYj) = l
Inj_X = function(ti,Yt,Zt,Xt,p,q,d){
yti = vector(mode = "numeric")
for (w in 1:p) {yti = c(yti,Yt[,ti - w])}
xti = vector(mode = "numeric")
for (w in 1:q) {xti = c(xti,Xt[,ti - w])}
zti = vector(mode = "numeric")
for (w in 1:d) {zti = c(zti,Zt[ti - w])}
if (q == 0 & d != 0) {
wtj = c(1,yti,zti)
}else if (d == 0 & q != 0) {
wtj = c(1,yti,xti)
}else if (d == 0 & q == 0) {
wtj = c(1,yti)
}else{
wtj = c(1,yti,xti,zti)}
return(wtj)
}
Inj_X = Vectorize(Inj_X,vectorize.args = "ti")
for (lj in 1:l) {
p = pjmax[lj]
q = qjmax[lj]
d = djmax[lj]
maxj = max(p,q,d)
Inj = which(Ind == lj)
Inj = Inj[Inj > maxj]
Nrg[lj] = length(Inj)
Yj = matrix(Yt[,Inj],nrow = k,ncol = Nrg[lj])
if (identical(Inj,integer(0))) {
Xj = matrix(nrow = 0,ncol = etam[lj]*k)
}else{
Wj = sapply(Inj,Inj_X,Yt = Yt,Zt = Zt,Xt = Xt,p = p,q = q,d = d)
Xj = t(Wj) %x% diag(k)[1,]
if (k != 1) {for (s in 2:k) {Xj = cbind(Xj,t(Wj) %x% diag(k)[s,])}}
}
listaXj[[lj]] = Xj
listaYj[[lj]] = Yj
}
return(list(Nrg = Nrg,listaX = listaXj,listaY = listaYj))
}
lists = compiler::cmpfun(lists)
### Likelihood function for Yt
fycond = function(ir, listar, gamma, theta, sigma){
acum = 0
l = length(listar$listaY)
Nrg = listar$Nrg
for (lj in 1:l) {
yj = c(listar$listaY[[lj]])
Xj = listar$listaX[[lj]]
acum = acum + t(yj - Xj %*% diag(gamma[[lj]][,ir]) %*% theta[[lj]][,ir]) %*%
{diag(Nrg[lj]) %x% solve(sigma[[lj]][[ir]])} %*%
(yj - Xj %*% diag(gamma[[lj]][,ir]) %*% theta[[lj]][,ir])
}
sigmareg = lapply(sigma,function(x){x[[ir]]})
cte = prodB(Brobdingnag::as.brob(sapply(sigmareg,function(x){
return(c(determinant(x,logarithm = FALSE)$modulus))}))^{-Nrg/2})
val = cte*exp(-1/2*Brobdingnag::as.brob(acum))
return(val)
}
fycond = compiler::cmpfun(fycond)
### Previous iterations
fill = function(m, iter = 500, burn = 1000, ...){
cat('m =',m,'... \n')
i = 1
ordersm = list(pj = rep(maxpj,m),qj = rep(maxqj,m),dj = rep(maxdj,m))
etam = 1 + ordersm$pj*k + ordersm$qj*nu + ordersm$dj
ini_obj_m = mtarinipars(tsregime_obj = ini_obj$tsregime_obj,
list_model = list(pars = list('l' = m),orders = ordersm),method = method)
#Parameters priori
theta0Pm = lapply(ini_obj_m$init$Theta,function(x){x$theta0j})
sigma0Pm = lapply(ini_obj_m$init$Theta,function(x){x$cov0j})
S0Pm = lapply(ini_obj_m$init$Sigma,function(x){x$S0j})
nu0Pm = lapply(ini_obj_m$init$Sigma,function(x){x$nu0j})
pij0Pm = ini_obj_m$init$Gamma
# previous iterations
par = mtarstr(ini_obj = ini_obj_m, niter = iter, chain = TRUE,burn = burn)
#Parameters pseudo Theta
theta0Sm = lapply(par$estimates$Theta,function(x){x[,2]})
sigma0Sm = lapply(par$Chain$Theta,function(x){stats::cov(t(x))})
#Parameters pseudo Sigma
S0Sm = lapply(par$regime,function(x){x$sigma})
nu0Sm = lapply(1:m,function(x){1000})
#Parameters pseudo Gamma
pij0Sm = lapply(par$Chain$Gamma,function(x){apply(x,1,mean)})
#Parameters pseudo R
if (m != 1) {
rmean0Sm = par$estimates$r[,2]
rcov0Sm = stats::cov(t(par$Chain$r))
}else{
rmean0Sm = rcov0Sm = NULL
}
# cadenas y primeros valores
theta_iter = sigma_iter = gam_iter = Dj = Rj = vector('list',m)
if (method == 'SSVS') {
tauij = itauij = cij = vector('list',m)
}
if (m != 1) {
r_iter = matrix(ncol = niter_m + burn_m,nrow = m - 1)
r_iter[,1] = c(stats::quantile(Zt, probs = 1/m*(1:{m - 1})))
r0iter = r_iter[,1]
}else{
r_iter = NULL
r0iter = 0
}
for (lj in 1:m) {
theta_iter[[lj]] = gam_iter[[lj]] = matrix(ncol = niter_m + burn_m,nrow = k*etam[lj])
sigma_iter[[lj]] = vector('list',length = niter_m + burn_m)
gam_iter[[lj]][,1] = rep(1,k*etam[lj])
if (method == 'KUO') {
theta_iter[[lj]][,1] = mvtnorm::rmvnorm(1,mean = theta0Pm[[lj]],sigma = sigma0Pm[[lj]])
}
if (method == 'SSVS') {
cij[[lj]] = ini_obj_m$init$Theta[[paste0('R',lj)]]$Cij
tauij[[lj]] = ini_obj_m$init$Theta[[paste0('R',lj)]]$Tauij
itauij[[lj]] = cij[[lj]]*tauij[[lj]]
Dj[[lj]] = diag(itauij[[lj]])
Rj[[lj]] = diag(k*etam[lj])
theta_iter[[lj]][,1] = mvtnorm::rmvnorm(1,mean = rep(0,k*etam[lj]),sigma = Dj[[lj]] %*% Rj[[lj]] %*% Dj[[lj]])
}
sigma_iter[[lj]][[1]] = MCMCpack::riwish(v = nu0Pm[[lj]],S = S0Pm[[lj]])
}
# LISTS
if (method == 'SSVS') {
iniP = list(Theta = list(mean = theta0Pm, cov = sigma0Pm), Sigma = list(cov = S0Pm,gl = nu0Pm),
Gamma = list(prob = pij0Pm, itauij = itauij, tauij = tauij, cij = cij, Rj = Rj))
}else{
iniP = list(Theta = list(mean = theta0Pm, cov = sigma0Pm), Sigma = list(cov = S0Pm,gl = nu0Pm),
Gamma = list(prob = pij0Pm))
}
iniS = list(Theta = list(mean = theta0Sm,cov = sigma0Sm), Sigma = list(cov = S0Sm,gl = nu0Sm),
Gamma = list(prob = pij0Sm), r = list(mean = rmean0Sm, cov = rcov0Sm))
listchain = list(Theta = theta_iter, Sigma = sigma_iter,
Gamma = gam_iter, r = r_iter)
listr = lists(l = m,r = r0iter,pjmax = ordersm$pj,qjmax = ordersm$qj,djmax = ordersm$dj, etam)
return(list(i = i,orders = ordersm,Priori = iniP,Pseudo = iniS,Chain = listchain,listr = listr,par = par))
}
fill = compiler::cmpfun(fill)
### Function to compute posterioris
updatelist = function(l, ...){
rgamber = function(pos, reg, ig, ...){
gam_j = gam_iter
gam_j[[reg]][pos,ig] = 1
pycond1 = fycond(ig,listaj,gam_j,theta_iter,sigma_iter)
gam_j[[reg]][pos,ig] = 0
pycond0 = fycond(ig,listaj,gam_j,theta_iter,sigma_iter)
if (method == 'KUO') {
aij = pycond1*pij[[reg]][pos]
bij = pycond0*(1 - pij[[reg]][pos])
}else if (method == 'SSVS') {
Xj = listaj$listaX[[reg]]
yj = c(listaj$listaY[[reg]])
gam_j[[reg]][pos,ig] = 1
itauij[[reg]][gam_j[[reg]][,ig] == 0] = tauij[[reg]][gam_j[[reg]][,ig] == 0]
itauij[[reg]][gam_j[[reg]][,ig] == 1] =
cij[[reg]][gam_j[[reg]][,ig] == 1]*tauij[[reg]][gam_j[[reg]][,ig] == 1]
Dj[[reg]] = diag(itauij[[reg]])
pthetacond1 = dmnormB(x = theta_iter[[reg]][,ig],mean = rep(0,k*etam[reg]),sigma = Dj[[reg]] %*% Rj[[reg]] %*% Dj[[reg]])
aij = pycond1*pthetacond1*pij[[reg]][pos]
gam_j[[reg]][pos,ig] = 0
itauij[[reg]][gam_j[[reg]][,ig] == 0] = tauij[[reg]][gam_j[[reg]][,ig] == 0]
itauij[[reg]][gam_j[[reg]][,ig] == 1] =
cij[[reg]][gam_j[[reg]][,ig] == 1]*tauij[[reg]][gam_j[[reg]][,ig] == 1]
Dj[[reg]] = diag(itauij[[reg]])
pthetacond0 = dmnormB(x = theta_iter[[reg]][,ig],mean = rep(0,k*etam[reg]),sigma = Dj[[reg]] %*% Rj[[reg]] %*% Dj[[reg]])
bij = pycond0*pthetacond0*(1 - pij[[reg]][pos])
}
return(stats::rbinom(1,size = 1,prob = as.numeric((aij)/(aij + bij))))
}
rgamber = Vectorize(rgamber,"pos")
listPr = listm[[paste0('m',l)]]
Dj = vector('list')
Rj = vector('list')
i2 = listPr$i
#creando temporales
pjmax = listPr$orders$pj
qjmax = listPr$orders$qj
djmax = listPr$orders$dj
etam = 1 + pjmax*k + qjmax*nu + djmax
theta_iter = listPr$Chain$Theta
sigma_iter = listPr$Chain$Sigma
gam_iter = listPr$Chain$Gamma
r_iter = listPr$Chain$r
theta0j = listPr$Priori$Theta$mean
sigma0j = listPr$Priori$Theta$cov
S0j = listPr$Priori$Sigma$cov
nu0j = listPr$Priori$Sigma$gl
pij = listPr$Priori$Gamma$prob
if (method == 'SSVS') {
itauij = listPr$Priori$Gamma$itauij
tauij = listPr$Priori$Gamma$tauij
cij = listPr$Priori$Gamma$cij
Rj = listPr$Priori$Gamma$Rj
}
#iterations update
if (!is.null(r_iter)) {
listaj = lists(l,r_iter[,i2],pjmax,qjmax,djmax,etam)
}else{
listaj = lists(l,0,pjmax,qjmax,djmax,etam)
}
for (lj in 1:l) {
Xj = listaj$listaX[[lj]]
Yj = listaj$listaY[[lj]]
yj = c(Yj)
Nj = listaj$Nrg[lj]
if (method == 'SSVS') {
itauij[[lj]][gam_iter[[lj]][,i2] == 0] = tauij[[lj]][gam_iter[[lj]][,i2] == 0]
itauij[[lj]][gam_iter[[lj]][,i2] == 1] = cij[[lj]][gam_iter[[lj]][,i2] == 1]*tauij[[lj]][gam_iter[[lj]][,i2] == 1]
Dj[[lj]] = diag(itauij[[lj]])
theta0j = list()
theta0j[[lj]] = rep(0,k*etam[lj])
}else if (method == 'KUO') {
Dj[[lj]] = diag(k*etam[lj])
Rj[[lj]] = sigma0j[[lj]]}
Vj = solve(t(diag(gam_iter[[lj]][,i2])) %*% t(Xj) %*% {diag(Nj) %x% solve(sigma_iter[[lj]][[i2]])} %*% Xj %*% diag(gam_iter[[lj]][,i2]) + solve(Dj[[lj]] %*% Rj[[lj]] %*% Dj[[lj]]))
thetaj = Vj %*% {t(diag(gam_iter[[lj]][,i2])) %*% t(Xj) %*% {diag(Nj) %x% solve(sigma_iter[[lj]][[i2]])} %*% yj + solve(sigma0j[[lj]]) %*% theta0j[[lj]]}
theta_iter[[lj]][,i2 + 1] = mvtnorm::rmvnorm(1,mean = thetaj,sigma = Vj)
Hj = ks::invvec({Xj %*% diag(gam_iter[[lj]][,i2]) %*% theta_iter[[lj]][,i2 + 1]},nrow = k,ncol = Nj)
Sj = (Yj - Hj) %*% t(Yj - Hj)
sigma_iter[[lj]][[i2 + 1]] = MCMCpack::riwish(v = Nj + nu0j[[lj]],S = Sj + S0j[[lj]])
gam_iter[[lj]][,i2 + 1] = gam_iter[[lj]][,i2]
}
#gamma
for (jj in 1:l) {
gam_iter[[jj]][,i2 + 1] = rgamber(pos = 1:{k*etam[jj]},reg = jj,ig = i2 + 1)
}
#r
if (l != 1) {
if (i2 < 70) {
ek = mvtnorm::rmvnorm(1,mean = rep(0,l - 1),sigma = 0.5*diag(l - 1))
}else{
ek = stats::runif(l - 1,-abs(rini$val_rmh),abs(rini$val_rmh))
}
# ek = runif(l - 1,-abs(rini$val_rmh),abs(rini$val_rmh))
rk = r_iter[,i2] + ek
listark = lists(l,rk,pjmax,qjmax,djmax,etam)
pr = dmunif(rk,a,b)*fycond(i2 + 1,listark,gam_iter,theta_iter,sigma_iter)
px = dmunif(r_iter[,i2],a,b)*fycond(i2 + 1,listaj,gam_iter,theta_iter,sigma_iter)
alpha = min(1,as.numeric(pr/px))
if (alpha >= stats::runif(1)) {
r_iter[,i2 + 1] = rk
listr = listark
}else{
r_iter[,i2 + 1] = r_iter[,i2]
listr = listaj
}
}else{
listr = lists(l,0,pjmax,qjmax,djmax,etam)
}
listPr$Chain = list(Theta = theta_iter, Sigma = sigma_iter,Gamma = gam_iter, r = r_iter)
listPr$listr = listr
listPr$i = i2 + 1
return(listPr)
}
updatelist = compiler::cmpfun(updatelist)
rpseudo = function(l,...){
listPr = listm[[paste0('m',l)]]
i2 = listPr$i + 1
pjmax = listPr$orders$pj
qjmax = listPr$orders$qj
djmax = listPr$orders$dj
etam = 1 + pjmax*k + qjmax*nu + djmax
theta_iter = listPr$Chain$Theta
sigma_iter = listPr$Chain$Sigma
gam_iter = listPr$Chain$Gamma
r_iter = listPr$Chain$r
theta0jS = listPr$Pseudo$Theta$mean
sigma0jS = listPr$Pseudo$Theta$cov
S0jS = listPr$Pseudo$Sigma$cov
nu0jS = listPr$Pseudo$Sigma$gl
pijS = listPr$Pseudo$Gamma$prob
rmeanS = listPr$Pseudo$r$mean
rcovS = listPr$Pseudo$r$cov
for (lj in 1:l) {
theta_iter[[lj]][,i2] = mvtnorm::rmvnorm(1,mean = theta0jS[[lj]],sigma = sigma0jS[[lj]])
sigma_iter[[lj]][[i2]] = 1/nu0jS[[lj]]*MCMCpack::rwish(v = nu0jS[[lj]] ,S = S0jS[[lj]])
sigma_iter[[lj]][[i2]] = sigma_iter[[lj]][[i2]] %*% sigma_iter[[lj]][[i2]]
for (iga in 1:{k*etam[lj]}) {
gam_iter[[lj]][iga,i2] = stats::rbinom(n = 1,size = 1,prob = pijS[[lj]][iga])
}
}
if (l != 1) {
r_iter[,i2] = mvtnorm::rmvnorm(1,mean = rmeanS, sigma = as.matrix(rcovS))
listPr$listr = lists(l,r_iter[,i2],pjmax,qjmax,djmax,etam)
}else{
listPr$listr = lists(l,0,pjmax,qjmax,djmax,etam)
}
listPr$Chain = list(Theta = theta_iter, Sigma = sigma_iter,Gamma = gam_iter, r = r_iter)
listPr$i = i2
return(listPr)
}
rpseudo = compiler::cmpfun(rpseudo)
### Function to compute quotient
prodA = function(thetaym, thetaymp){
pgammaPn = pthetaPn = psigmaPn = Brobdingnag::as.brob(1)
pgammaPd = pthetaPd = psigmaPd = Brobdingnag::as.brob(1)
pgammaSn = pthetaSn = psigmaSn = Brobdingnag::as.brob(1)
pgammaSd = pthetaSd = psigmaSd = Brobdingnag::as.brob(1)
lm = length(thetaym$listr$Nrg)
lmp = length(thetaymp$listr$Nrg)
theta_iterm = thetaym$Chain$Theta
sigma_iterm = thetaym$Chain$Sigma
gam_iterm = thetaym$Chain$Gamma
r_iterm = thetaym$Chain$r
iAm = thetaym$i
iAmp = thetaymp$i
theta_itermp = thetaymp$Chain$Theta
sigma_itermp = thetaymp$Chain$Sigma
gam_itermp = thetaymp$Chain$Gamma
r_itermp = thetaymp$Chain$r
theta0jmp = thetaymp$Priori$Theta$mean
sigma0jmp = thetaymp$Priori$Theta$cov
S0jmp = thetaymp$Priori$Sigma$cov
nu0jmp = thetaymp$Priori$Sigma$gl
pijmp = thetaymp$Priori$Gamma$prob
theta0jm = thetaym$Priori$Theta$mean
sigma0jm = thetaym$Priori$Theta$cov
S0jm = thetaym$Priori$Sigma$cov
nu0jm = thetaym$Priori$Sigma$gl
pijm = thetaym$Priori$Gamma$prob
theta0jSm = thetaym$Pseudo$Theta$mean
sigma0jSm = thetaym$Pseudo$Theta$cov
S0jSm = thetaym$Pseudo$Sigma$cov
nu0jSm = thetaym$Pseudo$Sigma$gl
pijSm = thetaym$Pseudo$Gamma$prob
rmeanSm = thetaym$Pseudo$r$mean
rcovSm = thetaym$Pseudo$r$cov
theta0jSmp = thetaymp$Pseudo$Theta$mean
sigma0jSmp = thetaymp$Pseudo$Theta$cov
S0jSmp = thetaymp$Pseudo$Sigma$cov
nu0jSmp = thetaymp$Pseudo$Sigma$gl
pijSmp = thetaymp$Pseudo$Gamma$prob
rmeanSmp = thetaymp$Pseudo$r$mean
rcovSmp = thetaymp$Pseudo$r$cov
for (lj in 1:lmp) {
pgammaPn = pgammaPn*prodB(Brobdingnag::as.brob(stats::dbinom(gam_itermp[[lj]][,iAmp],size = 1,prob = pijmp[[lj]])))
pthetaPn = pthetaPn*dmnormB(theta_itermp[[lj]][,iAmp],mean = theta0jmp[[lj]], sigma = sigma0jmp[[lj]])
psigmaPn = psigmaPn*dwishartB(sigma_itermp[[lj]][[iAmp]], nu = nu0jmp[[lj]],S = solve(as.matrix(S0jmp[[lj]])))
pgammaSd = pgammaSd*prodB(Brobdingnag::as.brob(stats::dbinom(gam_itermp[[lj]][,iAmp],size = 1,prob = pijSmp[[lj]])))
pthetaSd = pthetaSd*dmnormB(theta_itermp[[lj]][,iAmp],mean = theta0jSmp[[lj]], sigma = sigma0jSmp[[lj]])
psigmaSd = psigmaSd*dwishartB(expm::sqrtm(sigma_itermp[[lj]][[iAmp]]), nu = nu0jSmp[[lj]],S = as.matrix(S0jSmp[[lj]]))
}
fn = fycond(iAmp,thetaymp$listr,thetaymp$Chain$Gamma,thetaymp$Chain$Theta,thetaymp$Chain$Sigma)
for (lj in 1:lm) {
pgammaPd = pgammaPd*prodB(Brobdingnag::as.brob(stats::dbinom(gam_iterm[[lj]][,iAm],size = 1,prob = pijm[[lj]])))
pthetaPd = pthetaPd*dmnormB(theta_iterm[[lj]][,iAm],mean = theta0jm[[lj]], sigma = sigma0jm[[lj]])
psigmaPd = psigmaPd*dwishartB(sigma_iterm[[lj]][[iAm]], nu = nu0jm[[lj]],S = solve(as.matrix(S0jm[[lj]])))
pgammaSn = pgammaSn*prodB(Brobdingnag::as.brob(stats::dbinom(gam_iterm[[lj]][,iAm],size = 1,prob = pijSm[[lj]])))
pthetaSn = pthetaSn*dmnormB(theta_iterm[[lj]][,iAm],mean = theta0jSm[[lj]], sigma = sigma0jSm[[lj]])
psigmaSn = psigmaSn*dwishartB(expm::sqrtm(sigma_iterm[[lj]][[iAm]]), nu = nu0jSm[[lj]],S = as.matrix(S0jSm[[lj]]))
}
if (lmp != 1) {
prPn = dmunif(r_itermp[,iAmp],a,b)
prSd = dmnormB(r_itermp[,iAmp],mean = rmeanSmp, sigma = as.matrix(rcovSmp))
}else{
prPn = prSd = 1
}
if (lm != 1) {
prPd = dmunif(r_iterm[,iAm],a,b)
prSn = dmnormB(r_iterm[,iAm],mean = rmeanSm, sigma = rcovSm)
}else{
prPd = prSn = 1
}
fd = fycond(iAm,thetaym$listr,thetaym$Chain$Gamma,thetaym$Chain$Theta,thetaym$Chain$Sigma)
# Calculo para el valor del cociente
vald = fd*(pgammaPd*pthetaPd*psigmaPd*prPd)*(pgammaSd*pthetaSd*psigmaSd*prSd)
valn = fn*(pgammaPn*pthetaPn*psigmaPn*prPn)*(pgammaSn*pthetaSn*psigmaSn*prSn)
val = valn/vald
return(val)
}
prodA = compiler::cmpfun(prodA)
# Previous runs
message('Running previous chains ... \n')
listm = vector('list')
if (parallel) {
n_dis = parallel::detectCores()
if (l0 > {n_dis - 1}) {
parallel = FALSE
message('Cannot use paralell (l0_max > n_dis) ... \n')
}else {
n_cores = l0
}
}
if (parallel) {
cat('Using parallel','\n')
cat('Number of CPU cores used:',n_cores,'\n')
micluster = parallel::makeCluster(n_cores)
doParallel::registerDoParallel(micluster)
funcParallel = function(i,iterprev){return(fill(m = i,iter = iterprev, burn = round(0.3*iterprev)))}
parallel::clusterEvalQ(micluster, library(BMTAR))
obj_S = list('ini_obj','burn_m','niter_m','chain_m','list_m',
'ordersprev','k','N','nu','method','fill','maxpj','maxqj','maxdj',
'Zt','Yt','Xt','Ut','lists','fycond','rdunif','dmunif','l0_min')
parallel::clusterExport(cl = micluster,varlist = obj_S,envir = environment())
s = parallel::parLapply(micluster,as.list(l0_min:l0), funcParallel, iterprev = iterprev)
cc = 1
for (o in l0_min:l0) {
listm[[paste0('m',o)]] = s[[cc]]
cc = cc + 1
}
parallel::stopCluster(micluster)
}else{
for (o in l0_min:l0) {
listm[[paste0('m',o)]] = fill(m = o,iter = iterprev, burn = round(0.3*iterprev))
}
}
if (NAIC) {
results = list(tsregim = ini_obj$tsregime_obj,list_m = listm)
for (i in l0_min:l0) {
results$NAIC[[paste0('m',i)]] = mtarNAIC(listm[[paste0('m',i)]]$par)
}
x_NAIC = unlist(sapply(results$NAIC,function(x){x$NAIC}))
results$NAIC_final_m = as.numeric(gsub('m','',names(x_NAIC[which.min(x_NAIC)])))
return(results)
}
# Code
m_iter = vector('numeric')
m_iter[1] = sample(l0_min:l0,1)
acepm = 0
pm_im = matrix(nrow = l0,ncol = niter_m + burn_m - 1)
message('Estimating number of regimes ... \n')
pb = pbapply::timerProgressBar(min = 2, max = niter_m + burn_m,style = 1)
for (im in 2:{niter_m + burn_m}) {
m_iter[im] = rdunif(m_iter[im - 1],l0)
# Generamos valor para Theta_m y Theta_mp
listm[[paste0('m',m_iter[im - 1])]] = updatelist(l = m_iter[im - 1])
listm[[paste0('m',m_iter[im])]] = rpseudo(l = m_iter[im])
# Calculo para la probabilidad alpha
val = prodA(thetaym = listm[[paste0('m',m_iter[im - 1])]],
thetaymp = listm[[paste0('m',m_iter[im])]])
# Evaluacion de el criterio
alpham = min(1,as.numeric(val))
if (alpham >= stats::runif(1)) {
m_iter[im] = m_iter[im]
acepm = acepm + 1
}else{
m_iter[im] = m_iter[im - 1]
}
pm_im[,im - 1] = table(factor(m_iter,levels = 1:l0))/im
pbapply::setTimerProgressBar(pb,im)
}
close(pb)
message('Saving results ... \n')
# exits
m_iter = m_iter[-c(1:burn_m)]
vecm = sort(table(m_iter)/length(m_iter)*100,decreasing = TRUE)
namest = c('first freq','second freq','third freq','fourth freq')
ls = as.numeric(names(vecm))
names(vecm) = paste(namest[1:length(vecm)],paste0(round(vecm,2),'%'))
vecm[1:length(vecm)] = ls[1:length(vecm)]
if (chain_m & list_m) {
results = list(tsregim = ini_obj$tsregime_obj,
list_m = listm ,m_chain = m_iter,
prop = pm_im,
estimates = vecm,final_m = vecm[[1]])
}else if (chain_m & !list_m) {
results = list(tsregim = ini_obj$tsregime_obj,
m_chain = m_iter,
prop = pm_im,
estimates = vecm,final_m = vecm[1])
}else if (!chain_m & list_m) {
results = list(tsregim = ini_obj$tsregime_obj,
list_m = listm,
prop = pm_im,
estimates = vecm,final_m = vecm[1])
}else if (!chain_m & !list_m) {
results = list(tsregim = ini_obj$tsregime_obj,
prop = pm_im,
estimates = vecm,
final_m = vecm[1])
}
for (o in l0_min:l0) {
results$NAIC[[paste0('m',o)]] = mtarNAIC(listm[[paste0('m',o)]]$par)
}
compiler::enableJIT(0)
class(results) = 'regime_number'
return(results)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.