R/bayes_mixAR.R

In mixAR: Mixture Autoregressive Models

## Do not edit this file manually.
## It has been automatically generated from mixAR.org.

## Two functions to calculate component-specific error terms

err <- function(AR, mu, y, z, pk){
    n <- length(y)
    p <- max(pk)
    g <- length(pk)
    e <- matrix(0, n-p, g)
    for (t in 1 : (n - p)){
        for(k in 1:g){
            e[t, k] <-
                (y[t+p] - mu[k] -
                 (y[(t + p - 1) : (t + p - pk[k])] - mu[k]) %*% AR[[k]]) * z[t,k]
        }
    }
    e
}

err_k <- function(AR, mu, y, z, p, pk){
    n <- length(y)
    e <- numeric(n - p)
    for (t in 1 : (n - p)){
        e[t] <- (y[t + p] - mu - (y[(t + p - 1) : (t + p - pk)] - mu) %*% AR) * z[t]
    }
    e
}

## Main function
bayes_mixAR <- function(y, model, fix_shift = FALSE, a = 0.2, c = 2, tau, nsim, burnin){
    
    d      <- rep(1,length(model@prob))
    n      <- length(y)
    prob   <- model@prob
    sigma  <- model@scale
    prec   <- 1 / sigma ^ 2
    g      <- length(model@prob)
    pk     <- model@arcoef@p
    p      <- max(pk)
    shift  <-
    mu     <- model@shift
    AR     <- model@arcoef@a 

    bk <- sapply(AR, function(x) 1 - sum(x))
    
    mu[which(shift != 0)] <- shift[which(shift != 0)] / (bk[which(shift != 0)])

    if(length(tau) == 1) tau <- rep(tau, g)
    if(length(tau) != 1 & length(tau) != g){
        stop("tau should have length equal to 1 or g")
    }

    ## Hyperparameters
    Ry   <- max(y) - min(y)
    zeta <- min(y) + Ry / 2
    ka   <- 1 / Ry
    b    <- 100 * a / (c * Ry^2)

    ## Matrices to store output
    pid    <- 
    sigmad <- 
    precd  <- 
    mud    <- matrix(nrow = nsim, ncol = g)

    lnames <- NULL
    ard    <- list()
    for(k in 1:g){
        ard[[k]] <- matrix(0, nrow = nsim, ncol = pk[k])
        lnames   <- c(lnames, paste("Component_", k, sep = ""))
    }
    names(ard) <- lnames

    intd   <- matrix(nrow = nsim, ncol = g)
    count  <- rep(0, g)

    for(j in 1:nsim){
        ## Draw z_t's
        samp1 <- sampZpi(y, pk, prob, mu, AR, sigma, 1, d)
        z  <- samp1$latentZ
        nk <- samp1$nk

        ## Update proportions
        prob <- samp1$mix_weights
        pid[j, ] <- prob

        ## Update mu
        if(fix_shift == FALSE){
            samp2 <- sampMuShift(y, pk, prec, nk, shift, z, AR, 1)
            mu    <- samp2$mu
            shift <- samp2$shift
        }
        intd[j, ] <- shift
        mud[j, ]  <- mu

        ## Update precision
        samp3 <- sampSigmaTau(y, pk, prec, nk, AR, mu, z, a, c, 1)
        
        lambda <- samp3$lambda
        prec   <- samp3$precision
        sigma  <- samp3$scale

        sigmad[j, ] <- sigma
        precd[j, ]  <- prec

        ## Update phi
        ratio <- rep(0, g)
        phi_k <- list()

        for(k in 1:g){
            phi_k[[k]] <- rnorm(pk[k], AR[[k]], tau[k])
            ratio[k] <-
                min(1, exp(- prec[k]/2 *
                           sum(err_k(phi_k[[k]], mu[k], y, z[ ,k], p, pk[k]) ^ 2 -
                               err_k(AR[[k]], mu[k], y, z[ ,k], p, pk[k]) ^ 2
                               )))
        }

        cc <- rep(0, g)  ## counter vector
        u <- runif(g)   ## to accept or reject a set of new AR values

        for(k in 1:g){
            if(u[k] > ratio[k]){  ## if rejected
                phi_k[[k]] <- AR[[k]] ## phi_k is set to current values
            }else {
                if(j > burnin) cc[k] <- cc[k] + 1
            }
        }## does not count acceptance rate over burnin period

        A <- new("MixARGaussian", prob = as.numeric(prob), scale = as.numeric(sigma),
                 arcoef = phi_k, shift = as.numeric(shift))
        if(isStable(A)){  ## check for stability
            for(k in 1:g){
                ard[[k]][j, ] <- phi_k[[k]] ## add row to output
            }
            AR <- phi_k ## update current values
            if(j > burnin) {count <- count + cc}
        } else{
            for(k in 1:g){ ## if model is not stable, reject
                ard[[k]][j, ] <- as.matrix(AR[[k]])
                prob <- pid[j, ] <- pid[j-1, ]
                intd[j, ] <- shift <- intd[j-1, ]
                mud[j, ] <- mu <- mud[j-1, ]
                sigmad[j, ] <- sigma <- sigmad[j-1, ]
                precd[j, ]  <- prec <- precd[j-1, ]
            }
        }
        
    }
  
    ard <- lapply(ard, function(x) as.matrix(x[-c(1:burnin), ]))

    list(mix_weights = pid[-c(1:burnin), ],
         scale       = sigmad[-c(1:burnin), ],
         precision   = precd[-c(1:burnin), ],
         shift       = intd[-c(1:burnin), ],
         mu          = mud[-c(1:burnin), ],
         ARcoeff     = ard,  ## now returns a matrix for each K, also when pk=1 
         acc_rate    = count / (nsim - burnin),
         n_samp      = nsim - burnin,
         LatentZ     = z,
         ncomp       = g,
         fix_shift   = fix_shift)
}