swh:1:snp:3a44eb759780145deea094ac2a25c5049546a085
Raw File
Tip revision: c2f8f5dad0723f8b0e85fd9fe7602962777136ef authored by Han Lin Shang on 04 June 2010, 13:24:11 UTC
version 1.8
Tip revision: c2f8f5d
lik.r
lik = function (par, y, nstate, errortype, trendtype, seasontype, damped,
                par.noopt, lowerb, upperb, opt.crit, nmse, bounds, m, 
                pnames, pnames2)
{
    names(par) <- pnames
    names(par.noopt) <- pnames2
    alpha <- c(par["alpha"], par.noopt["alpha"])["alpha"]
    if (is.na(alpha))
        stop("alpha problem!")
    if (trendtype != "N") {
        beta <- c(par["beta"], par.noopt["beta"])["beta"]
        if (is.na(beta))
            stop("beta Problem!")
    }
    else beta <- NULL
    if (seasontype != "N") {
        gamma <- c(par["gamma"], par.noopt["gamma"])["gamma"]
        if (is.na(gamma))
            stop("gamma Problem!")
    }
    else gamma <- NULL
    if (damped) {
        phi <- c(par["phi"], par.noopt["phi"])["phi"]
        if (is.na(phi))
            stop("phi Problem!")
    }
    else phi <- NULL
    if (!check.param(alpha, beta, gamma, phi, lowerb, upperb, bounds, m))
        return(1e+12)
    np <- length(par)
    init.state <- par[(np - nstate + 1):np]
    if (seasontype != "N")
        init.state <- c(init.state, m * (seasontype == "M") -
            sum(init.state[(2 + (trendtype != "N")):nstate]))
    if (seasontype == "M") {
        seas.states <- init.state[-(1:(1 + (trendtype != "N")))]
        if (min(seas.states) < 0)
            return(1e+08)
    }
    e <- pegelsresid.C(y, frequency(y), init.state, errortype,
        trendtype, seasontype, damped, alpha, beta, gamma, phi)
    if (is.na(e$lik))
        return(1e+08)
    if (e$lik < -1e+10)
        return(-1e+10)
    if (opt.crit == "lik")
        return(e$lik)
    else if (opt.crit == "mse")
        return(e$amse[1])
    else if (opt.crit == "amse")
        return(mean(e$amse[1:nmse]))
    else if (opt.crit == "sigma")
        return(mean(e$e^2))
}
back to top