R/mixSARfit.R

In GeoBosh/mixAR: Mixture Autoregressive Models

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

ui <- function(x, t, i){
    ## 2018-11-13 this is extremely inefficient, replaced by Georgi
    ##
    ## out <- numeric(length(i))
    ## for(k in 1:length(i)){
    ##     out[k] <- if (i[k] == 0)
    ##                   1
    ##               else
    ##                   x[t - i[k]]
    ## }
    ## out

    if(any(t <= i))
        stop("all elements of 't - i' must be positive")

    ## this is a one-liner but  may fail if t > length(x):
    ##     res <- ifelse(i == 0, 1, x[t - i])

    if(any(i == 0)){
        res <- rep.int(1, length(i))
        res[i != 0] <- x[t - i[i != 0]]
        res
    }else
        x[t - i]
}

## Version with intercept
mixSARfit <- function(y, model, est_shift = FALSE, tol = 10^-14){

    verbose <- interactive() && options("verbose")$verbose
    
    one_or_zero <- if(est_shift) 0 else 1

    n  <- length(y)

    spk <- model@arcoef@ps
    sp  <- max(spk)
    s   <- model@arcoef@s
    g   <- length(model@prob)
    pk  <- model@order
    p   <- max(pk)
    
    index <- ((sp * s) + 1):n
    oldlik <- cond_loglikS(model, y, index)
    if(verbose)
        cat("Initial vallogf: ", oldlik, "\n")

    nc <- sp + p + one_or_zero
    nck_all <- pk + spk + one_or_zero

    sigma0 <- model@scale
    pi0 <- model@prob

    diff <- 1
    count <- 0
    while(diff > tol){
        count <- count+1
        ## rows contain same time, columns contain same component
        e <- matrix(0, nrow = n - sp * s, ncol = g)
        for(k in 1:g){
            for(t in 1:(n - sp * s)){
                e[t, k] <- y[t + sp * s] -
                    model@shift[k] -
                    sum(model@arcoef@a[[k]] * y[(t + sp * s - 1):(t + sp * s - pk[k])]) -
                    sum(model@arcoef@as[[k]] * ui(y, t + sp * s, s * 1:spk[k]))
            }
        }
        dens.e <- dnorm(t( t(e) / sigma0 ))

        tau  <-  t(pi0 / sigma0 * t(dens.e)) / rowSums(t(pi0 / sigma0 * t(dens.e)))

        if(g == 2)
            pi <- c(sum(tau[ , 1]) / (n - sp * s), 1 - sum(tau[ , 1]) / (n - sp * s))
        else if(g > 2){
            pi <- colSums(tau[ , 1:(g - 1)]) / (n - sp * s)
                                        # @Davide: there was a space between '<' and '-'
                                        #     pi < -c(pi, 1 - sum(pi))
            pi <- c(pi, 1 - sum(pi))
        }

        sigma <- sqrt(colSums(tau * e ^ 2) / colSums(tau))

        phimat <- matrix(0, nrow = g, ncol = nc) # matrix to contain coefficients during computation

        for(k in 1:g){
            nck <- nck_all[k] 
            yy <- numeric(nck) 
            ## @Davide: is it guaranteed that pk[k] > 0 here?
            ## @Georgi: I assumed "MixAR" objects by default have minimum pk[k] = 1
            ##          I can rewrite this if needs to handle pk[k] = 0  
            ##
            ## @Davide: if your code assume this, raise an error - otherwise the program may
            ##          raise an error for unknown reason or, worse, silently produce wrong
            ##          results.
            ##
            ## This is yet another reason not to do everything everytime from scratch, since
            ## the provided functions take care of such things (even if they don't, they can
            ## be improved to do so.
            onepkk <- seq_len(pk[k]) # ==  1:pk[k] when pk[k] > 0
            sonespk <- s * seq_len(spk[k])

            ## Left hand side
            for(t in 1:(n - sp * s)){
                tsps <- t + sp * s
                tsps1mk <- tsps - onepkk # == (tsps - 1):(tsps - pk[k]) == tsps - 1:pk[k])
                instr <- y[tsps1mk]
                if(one_or_zero == 1)
                    instr <- c(1, instr)
                if(spk[k] > 0)  ## i.e. if seasonal coefficients
                    instr <- c(instr, ui(y, tsps, sonespk))

                stopifnot(all(tsps1mk == (tsps - 1):(tsps - pk[k]))) # for testing

                yy <- yy + tau[t, k] * y[tsps] * instr
            }

            ## Right hand side
            gamma <- matrix(0, nrow = nck, ncol = nck)
            phi <- numeric(nc)
            for(x in 1:nrow(gamma)){
                for(t in 1:(n - sp * s)){
                    tsps <- t + sp * s
                    tsps1mk <- tsps - onepkk # == (tsps - 1):(tsps - pk[k]) == tsps - 1:pk[k])
                    instr <- y[tsps1mk]
                    if(one_or_zero == 1)
                        instr <- c(1, instr)
                    if(spk[k] > 0)
                        instr <- c(instr, ui(y, tsps, sonespk))

                    if(x == one_or_zero) # i.e. x == 1 and there is intercept
                        gamma[x, ] <- gamma[x, ] + tau[t, k] * instr
                    else if(x <= pk[k] + one_or_zero)
                        gamma[x, ] <- gamma[x, ] +
                            tau[t, k] * y[tsps - x + one_or_zero] * instr
                    else if(spk[k] > 0)
                        gamma[x, ] <- gamma[x, ] +
                            tau[t, k] * (
                                if(one_or_zero == 1)
                                    y[tsps - (x - pk[k] - one_or_zero) * s]
                                else
                                    ui(y, tsps, (x - pk[k] - one_or_zero) * s)
                            ) * instr
                }
            }
            phi[1:nck] <- pseudoInverse(gamma) %*% yy
            phimat[k, ] <- phi
        }

        ## Formatting results into "MixAR" object
        ar1 <- ar2 <- vector(mode = "list", g) # list()
        for(k in 1:g){
            ar1[[k]] <- phimat[k, (1 + one_or_zero):(pk[k] + one_or_zero)]
            ar2[[k]] <- phimat[k, (pk[k] + 1 + one_or_zero):(pk[k] + spk[k] + one_or_zero)]
        }

        ar <- new("raggedCoefS", a = ar1, as = ar2, s = s)

        newmod <- if(one_or_zero == 1)
                      new("MixARGaussian", prob = pi, scale = sigma, arcoef = ar,
                          shift = phimat[, 1])
                  else
                      new("MixARGaussian", prob = pi, scale = sigma, arcoef = ar)

        newlik <- cond_loglikS(newmod, y, index)

        if(count %% 25 == 0  && verbose)
            cat("niter: ", count, "\tvallogf: ", newlik, "\n")
        if(is.nan(newlik)  && verbose)
            cat("!!!! Log-likelihood is NaN, maybe due to singularity.\n")

        ## Update stopping criterion and model
        diff <- abs((oldlik - newlik) / oldlik)
        oldlik <- newlik
        model <- newmod
        sigma0 <- model@scale
        pi0 <- model@prob
        if(count == 200)
            break  ## max number of iterations
    }
    if(verbose)
        cat("Final niter: ", count, "\tvallogf: ", newlik, "\n")
    
    list("model" = newmod, "vallogf" = newlik)
}