# HW 3 Q8 ssab.R JLCY, 19 Dec 2013

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

ssab = function(yy, M) {

    #### SSA

    N = length(yy)

    m1 = M - 1
    Np = N - M + 1
    Np1 = Np + 1

    ## create the lagged matrix
    topl = matrix(0, Np, M)
    for (i in 1:M) {
        i1 = i - 1 + 1
        i2 = Np + i - 1
        topl[, i] = yy[i1:i2]
    }

    # zz=scale(topl)

    zs = var(topl)
    zsvd = svd(zs)

    # Eigen Values.. - fraction variance
    lambdas = (zsvd$d/sum(zsvd$d))

    ### PCs
    At = topl %*% zsvd$u


    #### Reconstructed Components RCs
    Rpc = matrix(0, N, M)

    for (ipc in 1:M) {

        ## t = 1,M-1 mt = 1/t; lt = 1, ut =t
        lt = 1
        m1 = M - 1

        for (i in 1:m1) {

            i1 = i - lt + 1
            i2 = i - i + 1

            Rpc[i, ipc] = sum(At[i1:i2, ipc] * zsvd$u[lt:i, ipc])/i

        }

        ### t = M, N', mt = 1/M, lt = 1, ut = M
        lt = 1
        ut = M

        for (i in M:Np) {

            i1 = i - lt + 1
            i2 = i - M + 1

            Rpc[i, ipc] = sum(At[i1:i2, ipc] * zsvd$u[lt:ut, ipc])/M

        }


        ### t = N'+1, N, mt = 1/(N-t+1), lt = t-N+M, ut = M
        lt = 1
        ut = M

        for (i in Np1:N) {
            lt = i - N + M
            i1 = i - lt + 1
            i2 = i - M + 1
            mt = (N - i + 1)

            Rpc[i, ipc] = sum(At[i1:i2, ipc] * zsvd$u[lt:ut, ipc])/mt

        }

    }

    out = c()
    out$lambdas = lambdas
    out$eof = zsvd$u
    out$At = At
    out$Rpc = Rpc

    as.list(out)

}