R/armaFit.R

In fArma: Rmetrics - Modelling ARMA Time Series Processes

# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General
# Public License along with this library; if not, write to the
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
# MA  02111-1307  USA


################################################################################
# FUNCTION:               PARAMETER ESTIMATION:
#  armaFit                 Fits parameters for ARMA Time Series process
#  .arFit                  Internal function called by armaFit
#  .arimaFit               Internal function called by armaFit
#  .arfimaFit              Internal function called by armaFit
################################################################################


armaFit <- 
    function(formula, data, method = c("mle", "ols"), include.mean = TRUE, 
    fixed = NULL, title = NULL, description = NULL, ...)
{   
    # A function implemented by Diethelm Wuertz

    # Description:
    #   Fits Model Parameters for an ARMA Time Series Process

    # Arguments:
    #   method - "mle", "ols"

    # Notes:
    #   Valid formulas are:
    #       "ar" "ma" "arma" "arima" "arfima" (not documented "fracdiff")
    #       Note, "arma(p,q)" uses arima(p,0,q)

    # Details:
    #   R-base:
    #       arima(
    #           x,
    #           order = c(0, 0, 0),
    #           seasonal = list(order = c(0, 0, 0), period = NA),
    #           xreg = NULL,
    #           include.mean = TRUE,
    #           transform.pars = TRUE,
    #           fixed = NULL,
    #           init = NULL,
    #           method = c("CSS-ML", "ML", "CSS"),
    #           n.cond,
    #           optim.control = list(),
    #           kappa = 1e+06)
    #   Compare with SPlus:
    #       arima.mle(
    #           x,
    #           model,
    #           n.cond,
    #           xreg=NULL,
    #           ...)

    # FUNCTION:

    # Check for Method:
    # ar.method       = c("yw", "burg1", "burg2", "ols", "mle")
    # arma.method     = c("CSS")
    # arima.method    = c("CSS-ML", "ML", "CSS")
    # fracdiff.method = NA
    method = method[1]  # Don't use match.arg(methods)

    # Call:
    fit = NULL
    call = match.call()

    # Add to Fracdiff: h and M default Settings
    mf = match.call(expand.dots = TRUE)
    m = match("h", names(mf), 0)
    if (m == 0) h = -1 else h = eval(mf[[m]])
    fracdiff.h = h
    m = match("M", names(mf), 0)
    if (m == 0) M = 100 else M = eval(mf[[m]])
    fracdiff.M = M

    # Get Series:
    # DW 2005-09-03
    mf = match.call(expand.dots = FALSE)
    m = match(c("formula", "data"), names(mf), 0)
    mf = mf[c(1, m)]
    mf[[1]] = as.name(".modelSeries")
    mf$fake = FALSE
    mf$lhs = TRUE
    if (missing(data)) data = eval(parse(text = search()[2]), parent.frame())
    mf$data = as.matrix(data)
    x = eval(mf, parent.frame())
    x = ts = as.vector(x[, 1])

###     # Allow for univariate 'timeSeries' Objects:
###     # Added 2004-09-04 DW
###     if (class(ts) == "timeSeries") ts = as.vector(ts)

    # Which Model?
    # DW 2006-02-21
    K = length(formula)
    end = regexpr("\\(", as.character(formula[K]))-1
    tsmodel =  substr(as.character(formula[K]), 1, end)

    # Valid Model?
    valid = FALSE
    if (tsmodel == "ar" ) valid = TRUE
    if (tsmodel == "ma" ) valid = TRUE
    if (tsmodel == "arma") valid = TRUE
    if (tsmodel == "arima") valid = TRUE
    if (tsmodel == "arfima") valid = TRUE
    if (tsmodel == "fracdiff") valid = TRUE
    if (!valid) stop("Invalid Formula Specification")

    # Which Order?
    start = regexpr("\\(", as.character(formula[K]))+1
    end   = regexpr("\\)", as.character(formula[K]))-1
    order = substr(as.character(formula[K]), start, end)

    if (tsmodel == "arfima" || tsmodel == "fracdiff") {
        # method will be ignored ...
        pos = regexpr(",", order)
        p = as.integer(substr(order, 1, pos-1))
        q = as.integer(substr(order, pos+1, nchar(order)))
        order = c(p, q)
        tsmodel = "arfima"
    }

    if (tsmodel == "arima") {
        if (method == "mle") method = "CSS-ML"
        pos = regexpr(",", order)
        p = as.integer(substr(order, 1, pos-1))
        order = substr(order, pos+2, nchar(order))
        d = as.integer(substr(order, 1, pos-1))
        q = as.integer(substr(order, pos+1, nchar(order)))
        order = c(p = p, d = d, q = q)
    }

    if (tsmodel == "arma") {
        if (method == "mle") method = "CSS-ML"
        # "arma" uses "arima"
        pos = regexpr(",", order)
        p = as.integer(substr(order, 1, pos-1))
        q = as.integer(substr(order, pos+1, nchar(order)))
        order = c(p = p, d = 0, q = q)
        tsmodel = "arima"
    }

    if (tsmodel == "ar") {
        # if method is CSS-ML, CSS, or ML, then "ar" uses "arima":
        order = as.integer(order)
        if (method == "mle") method = "CSS-ML"
        if (method == "CSS-ML" | method == "CSS" | method == "ML") {
            p = order
            order = c(p = p , d = 0, q = 0)
            tsmodel = "arima"
        }
    }

    if (tsmodel == "ma") {
        # if method is CSS-ML, CSS, or ML, then "ma" uses "arima":
        if (method == "mle") method = "CSS-ML"
        order = as.integer(order)
        order = c(p = 0 , d = 0, q = order)
        tsmodel = "arima"
    }

    # Which Function?
    fun = match.fun(paste(".", tsmodel, "Fit", sep = ""))

    # Fit:
    fit = fun(x = ts, order = order, include.mean = include.mean,
        method = method[1], fixed = fixed, ...)

    # "ols" specific:
    if (method == "ols") {
        se.coef = unlist(fit$se.coef)
        if (include.mean){
            ols.mean = se.coef[1]
            fit$se.coef = c(se.coef[-1], ols.mean) }
    }
    fit$call = call
    fit$tsmodel = tsmodel
    fit$class = "fARMA"
    class(fit) = "list"

    # Add title and desription:
    if (is.null(title)) title = "ARIMA Modelling"
    if (is.null(description)) description = description()


    # Parameters:
    parameter = list(include.mean = include.mean, fixed = fixed)
    if (tsmodel == "arfima" || tsmodel == "fracdiff") {
        parameter$M = fracdiff.M
        parameter$h = fracdiff.h
    }

    # Return Value:
    new("fARMA",
        call = as.call(match.call()),
        formula = as.formula(formula),
        method = as.character(method),
        parameter = parameter,
        data = list(x = x),
        fit = fit,
        residuals = list(residuals = fit$residuals),
        fitted = list(fitted = fit$fitted.values),
        title = as.character(title),
        description = as.character(description) )
}


# ------------------------------------------------------------------------------


.arFit <- 
    function(x, order, include.mean, fixed = NULL,
    method = c("yw", "burg1", "burg2", "ols", "mle"), M = NULL, h = NULL, ...)
{   
    # A function implemented by Diethelm Wuertz

    # Description:
    #   Fits an AR time series model

    # Note:
    #   Calls ar() from R-stats.

    # FUNCTION:

    # Fit:
    call = match.call()
    var.method = as.integer(substring(method[1], 5, 5))
    method = substr(method[1], 1, 4)
    fit = ar(x = x, aic = FALSE, order.max = order, method = method,
        var.method = var.method, demean = include.mean,
        intercept = include.mean, ...)

    # Add and Modify:
    fit$call = call
    fit$tstitle = paste("AR(",
        as.character(order), ") with method: ", method, sep = "")
    fit$order = order

    # Residuals:
    fit$residuals = fit$resid
    fit$fitted.values = x - fit$resid
    fit$sigma2 = fit$var.pred

    # Coefficients:
    if (method == "ols") {
        fit$coef = fit$ar[,,1]
    } else {
        fit$coef = fit$ar
    }
    names(fit$coef) = c(paste("ar", 1:order, sep=""))

    # Mean:
    if (include.mean) {
        coeff = c(fit$coef, fit$x.mean)
        names(coeff) = c(names(fit$coef), "intercept")
        fit$coef = coeff
    }

    # Standard Errors:
    if (method == "ols") {
        fit$se.coef = fit$asy.se.coef
        n = sqrt(length(as.vector(fit$se.coef)))
        fit$var.coef = matrix(rep(NA, times = n*n), ncol = n)
    } else {
        fit$var.coef = fit$asy.var.coef
        fit$se.coef = sqrt(diag(fit$asy.var.coef))
        if (include.mean) {
            m = dim(fit$asy.var.coef)[1] + 1
            var.coef = matrix(rep(NA, times = m*m), m, m)
            for ( i in 1:(m-1) ) {
                for ( j in 1:(m-1) ) {
                    var.coef[i,j] = fit$var.coef[i,j]
                }
            }
            fit$var.coef = var.coef
            fit$se.coef = c(fit$se.coef, NA)
        }
    }

    # Add Data:
    fit$x = x

    # Return Value:
    fit
}


# ------------------------------------------------------------------------------


.arimaFit <- 
    function (x, order, include.mean, fixed,
    method = c("CSS-ML", "ML", "CSS"), M = NULL, h = NULL, ...)
{   
    # A function implemented by Diethelm Wuertz

    # Description:
    #   Fits an ARIMA time series model

    # Note:
    #   Calls arima() from R-stats.

    # FUNCTION:

    # Fit:
    call = match.call()
    fit = arima(x = x, order = order, method =  method[1],
        include.mean = include.mean, fixed = fixed, ...)

    # Add Title:
    fit$tstitle = paste("ARIMA(",
        as.character(order[1]), ",", as.character(order[2]), ",",
        as.character(order[3]), ") with method: ", method[1], sep = "")

    # Add Data:
    fit$x = x

    # Add Fitted Values:
    fit$fitted.values = fit$x - fit$residuals

    # Add Standard Errors:
    fit$se.coef = sqrt(diag(fit$var.coef))

    # Add Call:
    fit$call = call

    # Return Value:
    fit
}


# ------------------------------------------------------------------------------


.arfimaFit =
    function (x, order, include.mean, fixed, method = "arfima",
    M = 100, h = -1, ...)
{   
    # A function implemented by Diethelm Wuertz

    # Description:
    #   Fits an ARFIMA (FRACDIFF) time series model

    # Arguments:
    #   x - time series for the ARIMA model
    #   nar - number of autoregressive parameters
    #   nma - number of moving average parameters
    #   ar - initial autoregressive parameters
    #   ma - initial moving average parameters
    #   dtol - desired accurcay for d, by default (and if
    #       negative), (4th root of machine precision)
    #       is used.  dtol will be changed internally if
    #       necessary
    #   drange - interval over which the likelihood function is
    #       to be maximized as a function of d
    #   h - finite difference interval
    #   M - number of terms in the likelihood approximation
    #       (see Haslett and Raftery 1989)

    # Note:
    #   A Builtin Copy from R's fracdiff Package
    #   Calls fracdiff() from R-fracdiff

    # FUNCTION:

    # Settings:
    call = match.call()
    nar = order[1]
    nma = order[2]
    ar = rep(NA, max(order[1], 1))
    ma = rep(NA, max(order[2], 1))
    dtol = .Machine$double.eps^0.25 # ~ 1.22e-4
    drange = c(0, 0.5)

    # fracdiff:
    if (any(is.na(x)))
        stop("missing values not allowed in time series")
    if (is.matrix(x) && ncol(x) > 2)
        stop("multivariate time series not allowed")
    n = length(x)
    npq = nar + nma
    npq1 = npq + 1
    lwork = max(npq+2*(n+M), 3*n+(n+6)*npq+npq%/%2+1, (3+2*npq1)*npq1+1)
    ar[is.na(ar)] = 0
    ma[is.na(ma)] = 0

    # if dtol < 0: the fortran code will choose defaults
    result = .Fortran("fracdf", as.double(x), as.integer(n),
        as.integer(M), as.integer(nar), as.integer(nma),
        dtol = as.double(dtol), drange = as.double(drange),
        hood = double(1), d = double(1), ar = as.double(ar),
        ma = as.double(ma), w = double(lwork), as.integer(lwork),
        info = integer(1), .Machine$double.xmin,
        .Machine$double.xmax, .Machine$double.neg.eps,
        .Machine$double.eps, PACKAGE = "fArma")
    if (result$info) switch(result$info,
        stop("insufficient workspace"),
        stop("error in gamma function"),
        stop("invalid MINPACK input"),
        warning(" Warning in gamma function"),
        warning(" Optimization failure"),
        warning(" Optimization limit reached"))
    hess = .Fortran("fdhpq",
         as.double(x), hess = double(npq1 * npq1), as.integer(npq1),
         result$w, PACKAGE = "fArma")$hess
    temp = .Fortran("fdcov", as.double(x), as.double(result$d),
         h = as.double(if (missing(h)) -1 else h), hd = double(npq1),
         cov = hess, as.integer(npq1), cor = hess, as.integer(npq1),
         se = double(npq1), result$w, info = integer(1),
         PACKAGE = "fArma")
    if (temp$info) switch(temp$info,
         warning(" Warning in gamma function"),
         warning(" Singular Hessian matrix"),
         warning(" Unable to compute correlation matrix"),
         stop("error in gamma function"))
    if (npq == 0) {
        result$ar = NULL
        result$ma = NULL }
    nam = "d"
    if (nar) nam = c(nam, paste("ar", 1:nar, sep = ""))
    if (nma) nam = c(nam, paste("ma", 1:nma, sep = ""))
    hess = matrix(hess, nrow = npq1, ncol = npq1, dimnames = list(nam, nam))
    hess[1, ] = temp$hd
    hess[row(hess) > col(hess)] = hess[row(hess) < col(hess)]
    se.ok = temp$info != 0 || temp$info < 3

    # Fitting Result:
    fit = list(
        log.likelihood = result$hood,
        d = result$d,
        ar = result$ar, ma = result$ma,
        covariance.dpq = array(temp$cov, c(npq1, npq1), list(nam, nam)),
        stderror.dpq = if (se.ok) temp$se, # else NULL
        correlation.dpq = if (se.ok) array(temp$cor, c(npq1, npq1)), # else NULL
        h = temp$h, d.tol = result$dtol, M = M, hessian.dpq = hess)

    # Add ts Title:
    fit$tstitle = paste("FRACDIFF(", as.character(order[1]), ",",
        as.character(order[2]), ") with method: ", method[1], sep = "")

    # Add Series:
    fit$x = x

    # Add Coefficients:
    fit$coef = c(fit$d, fit$ar, fit$ma)
    namesCoef = "d"
    if (order[1] > 0) {
        names.ar = c(paste("ar", 1:order[1], sep=""))
        namesCoef = c(namesCoef, names.ar) }
    if (order[2] > 0) {
        names.ma = c(paste("ma", 1:order[2], sep=""))
        namesCoef = c(namesCoef, names.ma) }
    names(fit$coef) = namesCoef
    fit$var.coef = fit$correlation.dpq

    # Add Fitted Values:
    n = 0:fit$M
    w = lgamma(-fit$d+n) - (lgamma(-fit$d)+lgamma(n+1))
    w = exp(w)
    fit$fitted.values = filter(fit$x, w, sides = 1)

    # Add Residuals:
    fit$residuals = x - fit$fitted.values

    # Add Standard Errors:
    fit$se.coef = fit$stderror.dpq

    # Add fracdiff Parameters:
    fit$fracdiff = c(M, h)

    # Add Call:
    fit$call = call

    # Return Value:
    fit
}


################################################################################