arima-simulation.R

library(polynom)


#' Generation of AR coefficients
#'
#' @param order : AR order.
#'
#' @return `numeric` array of AR coefficients.
#' @export
generate.ar.coef <- function(order) {
    if (order <= 0) {stop('Order coefficients must be greater than 0')}
    
    valid <- FALSE 
    
    # Initialize a loop to find random coefficients which fulfill stationary 
    # and causality conditions
    while (!valid) {
        # First generate values for the coefficients (zero mean)
        if (order > 1) {
            sign <- ifelse(runif(order-1) < 0.2, -1, 1)
            first_values <- sign*rnorm(order-1, 0.3, 0.3)
        } else {
            first_values <- c()
        }
        
        # Genreate the value of the last coefficient (non zero mean since 
        # it must be significative)
        sign <- ifelse(runif(1) < 0.1, -1, 1)
        last_value <- sign*rnorm(1, 0.8, 0.05) 
        
        # Concatenate arrays
        ar <- append(first_values, last_value)
        
        # Check stationary and causality conditions
        z <- solve(polynomial(coef=append(1, -ar)), 0)
        
        if (!any(abs(z) <= 1)) {
            valid <- TRUE
        }
    }
    return(ar)  
}

#' Generation of MA coefficients
#'
#' @param order : MA order.
#'
#' @return `numeric` array of MA coefficients.
#' @export
generate.ma.coef <- function(order) {
    if (order <= 0) {stop('MA order must be grater than 0')}
    valid <- FALSE
    
    while (!valid) {
        if (order > 1) {
            sign <- ifelse(runif(order-1) < 0.8, -1, 1)
            first_values <- sign*rnorm(order-1, 0.3, 0.07)
        } else {
            first_values <- c()
        }
        
        sign <- ifelse(runif(1) < 0.1, -1, 1)
        last_value <- sign*rnorm(1, 0.6, 0.05)
        
        ma <- append(first_values, last_value)
        
        z <- solve(polynomial(coef=append(1, ma)), 0)
        
        if (!any(abs(z) <= 1)) {
            valid <- TRUE
        }
    }
    return(ma)
}


comb <- function(n, x) {
    return(factorial(n) / factorial(n-x) / factorial(x))
}

B.binom <- function(X, d, t) {
    binom_coeff <- comb(d, 1:d)
    Xwindow <- rev(window(X, start=t-d, end=t-1))
    return(sum(binom_coeff*((-Xwindow)**(1:d))))
}


#' Simulate time series generated by an ARIMA process
#'
#' @param model : Order (p,d,q) of the ARIMA process.
#' @param n: Number of observations.
#' @param with.constant: Either to use constant.
#' @param x_ini_mean: Mean to generate a random constant.
#'
#' @return `ts` Time series.
#' @export
sim.arima <- function(model=list(p=0, d=0, q=0), n=1000, with.constant=FALSE, 
                      x_ini_mean=0) {
    if ('ar' %in% names(model)) {
        model$p <- length(model[['ar']])
    }
    if ('ma' %in% names(model)) {
        model$q <- length(model[['ma']])
    }
    if (!'d' %in% names(model)) {
        model$d <- 0
    }
    
    if (with.constant && (model$d > 0)) {stop('Si d > 0 no puede haber constante')}
    
    # Generation of white noise
    a <- ts(rnorm(n+model$q, 0, 0.05), start=-model$q, end=n-1)
    
    # Generation of a constant
    c <- ifelse(with.constant, rnorm(1, x_ini_mean, 0.1), 0)
    
    # In case no order is specified, return gaussian noise
    if ((model$p==0) && (model$q==0) && (model$d==0)) {
        x <- c + a
        return(list(X=x, ar=NULL, ma=NULL, c=c, d=0))
    } 
    
    # Generation of the first values of the time series
    if (model$p > 0) {
        X <- ts(rnorm(max(c(model$p, model$d)), 0, 0.5), end=-1)
        # X <- ts(rep(0, max(c(model$p, model$d))), end=-1)
    } else {
        X <- c()
    }
    
    if (!'ar' %in% names(model) && (model$p>0)) {
        model$ar <- generate.ar.coef(model$p)
    }
    
    if (model$p > 0) {
        apply.ar <- function(X, t) {
            Xwindow <- rev(window(X, start=t-model$p, end=t-1))
            return(sum(model$ar*Xwindow))
        }
    } else { apply.ar <- function(X, t) {return(0)} }
    
    if (!'ma' %in% names(model) && (model$q > 0)) {
        model$ma <- generate.ma.coef(model$q)
    }
    
    
    if (model$q > 0) {
        apply.ma <- function(a, t) {
            awindow <- rev(window(a, start=t-model$q, end=t-1))
            return(sum(model$ma*awindow))
        }
    } else { apply.ma <- function(x, t) {return(0)} }
    
    
    # Update values
    for (t in 0:(n-1)) {
        ar_part <- apply.ar(X, t)
        ma_part <- apply.ma(a, t)
        
        Xt <- c +  ar_part + ma_part + window(a, start=t, end=t)
        X <- ts(append(X, Xt), end=t)
    }
    
    X <- window(X, start=0, end=n-1)
    
    # Since X ~ proceso ARMA(p, q), we need to obtain Z such that X_t = Z_t - Z_{t-1}
    if (model$d > 0) {
        for (i in 1:model$d) {
            Z <- ts(rnorm(1, 0, 0.01), start=-1)
            for (t in 0:(n-1)) {
                Zt <-  window(X, start=t, end=t) - B.binom(Z, 1, t)
                Z <- ts(append(Z, Zt), start=-1, end=t)
            }
            X <- window(Z, start=0, end=n-1) 
        }
    }
    
    result <- list(ar=model$ar, ma=model$ma, c=c, X=X, d=model$d)
    return(result)
    
    
    # Testing
    ajuste <- auto.fit.arima(result$X, show_info=F)
     
    if (is_valid(ajuste)) {
       valid <- (model$p == ajuste$arma[1]) & (model$q == ajuste$arma[2]) & (model$d == ajuste$arma[6]) & 
           (('intercept' %in% names(ajuste$coef)) == with.constant)
    } else {
       valid <- F
    }
    
    
    if (valid) {
       return(result)
    } else {
       return(sim.arima(model, n, with.constant, x_ini_mean))
    }
}