R/mixARcalc.R

In GeoBosh/mixAR: Mixture Autoregressive Models

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

cond_loglik <- function(model, x, index){                 # todo: the same for the derivative?
    if(missing(index))                                    #       or for the derivative of the
        index <- (max(model@order) + 1) : length(x)       #       one in the M-step of EM.

    wrk <- mix_pdf(model, x, index)
    wrk <- log(wrk)
    wrk <- sum(wrk)
    wrk
}

## For use with seasonal coefficients.
## todo: inefficient, could adapt mix_pdf etc.. to make it quicker
cond_loglikS <- function(model, x, index){
    n <- length(x)
    g <- length(model@prob)
    if(missing(index))
        index <- (max(model@arcoef@ps * model@arcoef@s)+1):length(x)
    yhat <- matrix(0, nrow = length(index) , ncol = g) # rows contain same time,
                                                   # columns contain same component
    for(t in seq_along(index)){  ## Calculate predictions
        for(k in 1:g){
            yhat[t,k] <- model@shift[k] +
                sum(model@arcoef@a[[k]] * x[(index[t] - 1):(index[t] - model@order[k])]) +
                sum( model@arcoef@as[[k]] * ui(x,index[t], model@arcoef@s *
                                                             (1:model@arcoef@ps[k])) )
        }
    }
    pdf <- numeric(length(index))  ## Calculate density
    for(t in seq_along(index)){
        for(k in 1:g){
            pdf[t] <- pdf[t] + model@prob[k] * dnorm(x[index[t]], yhat[t, k], model@scale[k])
        }
    }
    ll <- log(pdf)  ## loglikelihood
    sum(ll)
}

lik_params_bounds <- function(model){
    k <- .nmix(model)            # number of mixture components
    n <- length(lik_params(model))     # effective # of parameters

    lower <- rep(-Inf, n)              # lower bounds
    lower[1:(k-1)] <- rep(0,k-1)       # probabilities are >= 0
    lower[k:(k+k-1)] <- rep(0,k)       # st.dev's are >=0

    upper <- rep(Inf, n)
    upper[1:(k-1)] <- rep(1,k-1)       # probabilities are <= 1

    list(lower = lower, upper = upper)
}

## rdist - list of functions, 
## z     - vector of positive integers ('regimes')
## return value - numeric vector of same length as z
mixARnoise_sim <- function(rdist, z){ 
    g <- length(rdist)          # g must be positive
    n <- length(z)
    if(g == 1){
        do.call(rdist[[1]], list(n)) # the dist is one and the same, z not needed in this case
    }else{
        wrk <- rep(as.numeric(NA), n)
        for(i in 1:g){
            ni <- sum(z == i)
            wrk[z == i] <- do.call(rdist[[i]], list(ni))
        }
        wrk
    }
}

## todo: tova e nay-dobre da stane generic.
## todo: change the default for nskip to 0?
##       set default for init using the unconditional mean?
mixAR_sim <- function(model, n, init, nskip = 100, flag = FALSE){
    sigma <- model@scale
    shift <- model@shift
    zdist <- model@prob
    K <- length(zdist) # number of components

    ## 'seas_flag' is used below at places were the seasonal case needs its own treatment
    seas_flag <- inherits(model@arcoef, "raggedCoefS")
    
    if(seas_flag){     
        s <- model@arcoef@s
        order <- matrix(0, nrow = K, ncol = max(model@order + model@arcoef@ps))
        for(i in seq_len(K)){
            order[i, 1:(model@order[i] + model@arcoef@ps[i])] <-
                c(seq_len(model@order[i]), s * seq_len(model@arcoef@ps[i]))
        }
    }else
        order <- model@order


    if(length(init) < max(order))
        stop("The length of argument 'init' must be at least equal to the maximal order.")
    resinit <- tail(init, max(order))   # lastn(init, max(order))
                                        # use only the last max(order) values  in init.
    ntot     <- nskip + n
    nresinit <- length(resinit)

    z <- sample.int(K, size = ntot, replace = TRUE, prob = zdist)

    noisedist <- noise_rand(model)
    res <- mixARnoise_sim(noisedist, z)

    res <- c(resinit, res)                # now prepend the initial values
    z   <- c(numeric(nresinit), z)

    ar <- model@arcoef
    ## needs efficient implementation, this is extremely inefficient (but straightforward)
    for(i in nresinit + (1:ntot)){
        k <- z[i]
        wrk <- sigma[k] * res[i] + shift[k]
                                        # 2018-12-09: added simulation from seasonal AR
        if(seas_flag){
            filt <- c(ar@a[[k]], ar@as[[k]])
            for(j in seq_along(order[k, which(order[k, ] !=0)])){
                wrk <- wrk + filt[j] * res[i - order[k, j]]
            }
        } else{
            filt <- ar[[k]]
            for(j in seq_len(order[k])){       # 2011-07-12: was 1:order[k]
                wrk <- wrk + filt[j] * res[i - j]
            }
        }
        res[i] <- wrk
    }
    if(flag) structure(lastn(res, n), regimes = lastn(z, n)) # return the last n values (could
    else     lastn(res, n)                                   # use tail()); if 'flag', attach
}                                                            # 'z' as attribute to the result.

# 2012-12-03: new
# quick fix for Shahadat, consolidiray later.
mixAny_sim <- function(model, n, init, nskip=100, flag = FALSE,
                       theta, # ma coef, a list
                       galpha0,   # alpha0[k], k=1,...,g.
                       galpha,   # garch alpha
                       gbeta    # garch beta
                       ){
    order <- model@order
    sigma <- model@scale
    shift <- model@shift
    zdist <- model@prob

    g <- length(model@prob)
    thetaorder <- sapply(theta, length)
    alphaorder <- sapply(galpha, length)
    betaorder <- sapply(gbeta, length)
    baseind <- max(betaorder) + 1   # baseind - j is tauinv[t-j]
                                                           # column k is for the k-th component
    tauinvmat <- matrix(rep(galpha0), byrow = TRUE, ncol = g, nrow = max(betaorder))


    K <- length(zdist) # number of components

    if(length(init) < max(order))
        stop("The length of argument 'init' must be at least equal to the maximal order.")
    resinit <- tail(init, max(order,alphaorder,thetaorder))   # lastn(init, max(order))
                                        # use only the last max(order) values  in init.

    ntot     <- nskip + n
    nresinit <- length(resinit)

    z   <- sample.int(K, size=ntot, replace=TRUE, prob=zdist)

    noisedist <- noise_rand(model)
    res <- mixARnoise_sim(noisedist, z)

    eps <- c(numeric(nresinit), res)  # needed for the moving averages

    res <- c(resinit, res)                # now prepend the initial values
    z   <- c(numeric(nresinit), z)

    ar <- model@arcoef
    ## needs efficient implementation, this is extremely inefficient (but straightforward)
    for(i in nresinit+(1:ntot)){
        k <- z[i]

        tauinv <- galpha0  # compute new tauinv              # eq. (5.2), Shahadat
        for(icomp in 1:g){
            for(j in seq_len(alphaorder[icomp])){
                tauinv[icomp] <- tauinv[icomp] + galpha[[icomp]][j]*eps[i-j]^2
            }

            for(j in seq_len(betaorder[icomp])){
                tauinv[icomp] <- tauinv[icomp] +
                                                gbeta[[icomp]][j]*tauinvmat[baseind - j, icomp]
            }
        }
        sigma[k] <- sqrt(tauinv[k])  # multiply epsilon by this below

        tauinvmat <- rbind(tauinvmat[-1, ], tauinv) # need only the last few tauinv's
                                                    # so drop the oldest

        filt <- ar[[k]]
        wrk <- sigma[k]*res[i] + shift[k]
        for(j in seq_len(order[k])){       # 2011-07-12: was 1:order[k]
            wrk <- wrk + filt[j]*res[i-j]
        }

        maco <- theta[[k]]           # ma part
        for(j in seq_len(length(maco))){
            wrk <- wrk + maco[j]*eps[i-j]
        }

        res[i] <- wrk
    }
    if(flag) structure(lastn(res,n), regimes = lastn(z,n))   # return the last n values (could
    else     lastn(res,n)                                    # use tail()); if 'flag', attach
}                                                            # 'z' as attribute to the result.