.packageName <- "mAr"
"mAr.eig" <-
function (A, C = NULL, ...) 
{
    m = dim(A)[1]
    p = (dim(A)[2])/m
    At = matrix(nrow = m * p, ncol = m * p)
    if (p == 1) 
        At = A
    else {
        At[seq(1, m), seq(1, m * p)] = A
        At[seq(m + 1, m * p), seq(1, m * p - m)] = diag(1, (p - 
            1) * m)
        At[seq(m + 1, m * p), seq(m * p - m + 1, m * p)] = 0
    }
    l = eigen(At)$values
    V = eigen(At)$vectors
    if ((any(Mod(l) > 1))) 
        warning("unstable AR model")
    a = matrix(nrow = 1, ncol = dim(V)[2])
    b = matrix(nrow = 1, ncol = dim(V)[2])
    St = matrix(nrow = dim(V)[2], ncol = dim(V)[2])
    for (j in seq(1, dim(V)[2])) {
        a = Re(V[, j])
        b = Im(V[, j])
        ph = 0.5 * atan(2 * sum(a * b)/(b %*% b - a %*% a))
        na = sqrt(sum((cos(ph) * a - sin(ph) * b)^2))
        nb = sqrt(sum((sin(ph) * a + sin(ph) * b)^2))
        if (nb > na && ph < 0) {
            ph = ph - pi/2
        }
        if (nb > na && ph > 0) {
            ph = ph + pi/2
        }
        St[, j] = V[, j] %*% exp((0 + (0+1i)) * ph)
    }
    S = St[seq(1 + (p - 1) * m, p * m), ]
    StInv = solve(St)[, seq(1, m)]
    tau = matrix(nrow = 1, ncol = m * p)
    per = matrix(nrow = 1, ncol = m * p)
    exctn = matrix(nrow = 1, ncol = m * p)
    for (i in seq(1, m * p)) {
        tau[i] = -2/log((abs(l[i]))^2)
        a = Re(l[i])
        b = Im(l[i])
        if (identical(b, 0) && a >= 0) {
            per[i] = Inf
        }
        if (identical(b, 0) && a < 0) {
            per[i] = 2
        }
        else {
            per[i] = 2 * pi/abs(atan2(b, a))
        }
    }
    M = cbind(periods = as.vector(per), dampTime = as.vector(tau))
    result = M[c(order(M[, 2], decreasing = TRUE)), ]
    return(list(modes = result, eigv = S))
}
"mAr.est" <-
function (x, p, ...) 
{
    x = as.matrix(x)
    n = dim(x)[1]
    m = dim(x)[2]
    ne = n - p
    np = m * p + 1
    if (ne <= np) {
        stop("time series too short")
    }
    K = matrix(nrow = ne, ncol = np + m)
    K[, 1] = rep(1, ne)
    for (j in 1:p) {
        K[, seq(2 + m * (j - 1), 1 + m * j)] = data.matrix(x[seq(p - 
            j + 1, n - j), 1:m])
    }
    K[, seq(np + 1, np + m)] = data.matrix(x[seq(p + 1, n), 1:m])
    q = ncol(K)
    delta = (q^2 + q + 1) * (.Machine$double.eps)
    scale = sqrt(delta) * sqrt(apply(K^2, 2, sum))
    R = qr.R(qr((rbind(K, diag(scale)))), complete = TRUE)
    R22 = R[seq(np + 1, np + m), seq(np + 1, np + m)]
    logdp = 2 * log(abs(prod(diag(R22))))
    sbc = logdp/m - log(ne) * (ne - np)/ne
    R11 = R[seq(1, np), seq(1, np)]
    R12 = R[seq(1, np), seq(np + 1, np + m)]
    R22 = R[seq(np + 1, np + m), seq(np + 1, np + m)]
    R11[, 1] = R11[, 1] * max(scale[2:(np + m)])/scale[1]
    B = t(solve(R11) %*% R12)
    w = B[, 1] * max(scale[2:(np + m)])/scale[1]
    A = B[, 2:np]
    C = (t(R22) %*% R22)/(ne - np)
    t = seq(1, (n - p))
    res = matrix(nrow = (n - p), ncol = m)
    res[t, seq(1, m)] = x[t + p, ] - (rep(1, n - p) %*% t(w))
    for (j in seq(1, p)) {
        res[t, seq(1, m)] = res[t, seq(1, m)] - (x[(t - j + p), 
            ] %*% t(A[, seq(m * j - m + 1, j * m)]))
    }
    return(list(SBC = sbc, wHat = w, AHat = A, CHat = C, resid = res))
}
"mAr.pca" <-
function (x, p, k = dim(x)[2], ...) 
{
    require(mAr)
    n = dim(x)[1]
    m = dim(x)[2]
    if (k <= 1) {
        stop("EOF subspace dimension <=1")
    }
    xcentred = matrix(nrow = n, ncol = m)
    xscaled = matrix(nrow = n, ncol = m)
    for (i in 1:m) {
        xcentred[, i] = x[, i] - apply(x, 2, mean)[i]
    }
    for (i in 1:m) {
        xscaled[, i] = xcentred[, i]/apply(xcentred, 2, sd)[i]
    }
    d = svd(xscaled)$d
    fraction.variance = cumsum(d^2)/sum(d^2)
    D = diag(svd(xscaled)$d[1:k])
    U = svd(xscaled)$u[, 1:k]
    V = svd(xscaled)$v[, 1:k]
    scores = U %*% D
    loadings = V
    mArModel = mAr.est(scores, p)
    A = mArModel$AHat
    C = mArModel$CHat
    SBC = mArModel$SBC
    modes = mAr.eig(A)$modes
    P = V %*% mAr.eig(A, C)$eigv
    return(list(order = p, SBC = SBC, fraction.variance = fraction.variance[1:k], 
        resid = mArModel$resid, eigv = P, modes = modes))
}
"mAr.sim" <-
function (w, A, C, N, ...) 
{
    m = dim(A)[1]
    p = (dim(A)[2])/m
    At = matrix(nrow = m * p, ncol = m * p)
    if (p == 1) {
        At = A
    }
    else {
        At[seq(1, m), seq(1, m * p)] = A
        At[seq(m + 1, m * p), seq(1, m * p - m)] = diag(1, (p - 
            1) * m)
        At[seq(m + 1, m * p), seq(m * p - m + 1, m * p)] = 0
    }
    l = (eigen(At, only.values = TRUE))$values
    if ((any(Mod(l) > 1))) 
        warning("unstable AR model")
    nd = 1000
    U = chol(C)
    require(MASS)
    noisevec = mvrnorm(nd + N, rep(0, m), C)
    matw = rep(1, nd + N) %*% t(w)
    vec = noisevec + matw
    if (any(w != 0)) {
        B = diag(1, m)
        for (j in seq(1, p)) {
            B = B - A[, seq(m * j - m + 1, j * m)]
        }
        mproc = as.vector(solve(B) %*% w)
        xi = (matrix(1, nrow = p, ncol = 1)) %*% mproc
    }
    else {
        xi = matrix(nrow = p, ncol = m)
        xi[, ] = 0
    }
    u = matrix(nrow = p + nd + N, ncol = m)
    u[seq(1, p), seq(1, m)] = xi
    u[seq(p + 1, p + nd + N), seq(1, m)] = 0
    Atr = t(A)
    x = matrix(ncol = m, nrow = p)
    for (k in seq(p + 1, nd + N + p)) {
        for (j in seq(1, p)) {
            x[j, ] = u[k - j, ] %*% Atr[seq(m * j - m + 1, m * 
                j), ]
        }
        u[k, ] = as.matrix(apply(x, 2, sum) + vec[k - p, ])
    }
    v = u[seq(nd + p + 1, nd + p + N), ]
    simulated = data.frame(v[, seq(1, m)])
    return(simulated)
}
