Fls <-
function(x, ord)
{
        p <- ord
        n <- length(x)
        xf <- x[(p + 1):n]
        FF <- matrix(0, n - p, p)
        for(i in 1:p) {
                FF[, i] <- x[(p - i + 1):(n - i)]
        }
        pf <- solve(t(FF) %*% FF) %*% t(FF) %*% xf
        varf <- xf - FF %*% pf
        list(ar = pf[,1], var.pred = (t(varf) %*% varf)[1,1]/(n - 2 * p), order=ord)
}


Bls <-
function(x, ord)
{
        p <- ord
        n <- length(x)
        xf <- x[1:(n - p)]
        FF <- matrix(0, n - p, p)
        for(i in 1:p) {
                FF[, i] <- x[(i + 1):(n - p + i)]
        }
        pf <- solve(t(FF) %*% FF) %*% t(FF) %*% xf
        varf <- xf - FF %*% pf
        list(ar = pf[,1], var.pred = (t(varf) %*% varf)[1,1]/(n - 2 * p), order=ord)
}

FBls <-
function(x, ord)
{
        p <- ord
        n <- length(x)
        xf <- x[(p + 1):n]
        FF <- matrix(0, n - p, p)
        for(i in 1:p) {
                FF[, i] <- x[(p - i + 1):(n - i)]
        }
        xb <- x[1:(n - p)]
        BB <- matrix(0, n - p, p)
        for(i in 1:p) {
                BB[, i] <- x[(i + 1):(n - p + i)]
        }
        pf <- solve(t(FF) %*% FF + t(BB) %*% BB) %*% (t(FF) %*% xf + t(BB) %*% 
                xb)
        varf <- xf - FF %*% pf
        varb <- xb - BB %*% pf
        varres <- sum(varf^2) + sum(varb^2)
        list(ar = pf[,1], var.pred = varres/(2 * n - 3 * p), order=ord)
}


lev.dur <-
function(x, ord)
{
        n <- length(x)
        s <- rep(0, n)
        for(lag in 0:(n - 1))
                s[lag + 1] <- (1/n) * sum(x[1:(n - lag)] * x[(1 + lag):n])
        p <- matrix(0, n, n)
        P <- rep(0, n)
        p[1, 1] <- s[2]/s[1]
        P[1] <- s[1] * (1 - p[1, 1] * p[1, 1])
        for(k in 2:ord) {
                p[k, k] <- (s[k + 1] - sum(p[1:(k - 1), (k - 1)] * s[k:2]))/P[k -
                        1]
                for(j in 1:(k - 1))
                        p[j, k] <- p[j, k - 1] - p[k, k] * p[k - j, k - 1]
                P[k] <- P[k - 1] * (1 - p[k, k] * p[k, k])
        }
        list(ar = p[1:ord, ord], var.pred = P[ord], order=ord)
}

lev.dur.tap <-
function(x, ord)
{
        n <- length(x)
        s <- rep(0, n)
        h <- spec.taper(rep(1, n), p = 0.25)
        h <- as.vector(h)/sqrt(sum(h^2))
        x <- h * x
        for(lag in 0:(n - 1))
                s[lag + 1] <- sum(x[1:(n - lag)] * x[(1 + lag):n])
        p <- matrix(0, n, n)
        P <- rep(0, n)
        p[1, 1] <- s[2]/s[1]
        P[1] <- s[1] - p[1, 1] * s[2]
        for(k in 2:ord) {
                p[k, k] <- (s[k + 1] - sum(p[1:(k - 1), (k - 1)] * s[k:2]))/P[k -
                        1]
                for(j in 1:(k - 1))
                        p[j, k] <- p[j, k - 1] - p[k, k] * p[k - j, k - 1]
                P[k] <- P[k - 1] * (1 - p[k, k] * p[k, k])
        }
        list(ar = p[1:ord, ord], var.pred = P[ord], order=ord)
}

YW8.plot <-
function(x1 = ar4.64, x2 = ar4.256, x3 = ar4.1024, ar4 = ar4.real)
{
        par(mfcol = c(3, 1), mar = c(4, 4, 0, 2), mgp = c(1.8, 0.85, 0))
        n <- 64
        x <- x1
        ar.ld <- lev.dur(x, ord = 8)
        ar.specld <- spec.ar(ar.ld)
        ar.spec <- spec.ar(list(ar = ar4, var.pred = 1,order=4))
        plot(ar.specld$freq, ar.specld$spec, type = "l", axes = F, xlab = "", ylab = "dB", ylim = range(-40, 60), lwd = 2)
        lines(ar.spec$freq, ar.spec$spec)
        text(0.5, 55, "N=64, Yule-Walker AR(8)", adj = 1, cex = 0.8)
        axis(2, at = c(-40, -20, 0, 20, 40, 60))
        box()
        n <- 256
        x <- x2
        ar.ld <- lev.dur(x, ord = 8)
        ar.specld <- spec.ar(ar.ld)
        ar.spec <- spec.ar(list(ar = ar4, var.pred = 1,order=4))
        f <- seq(0, 0.5, l = length(ar.specld))
        plot(ar.specld$freq, ar.specld$spec, type = "l", axes = F, xlab = "", ylab = "dB", ylim = range(-40, 60), lwd = 2)
        lines(ar.spec$freq, ar.spec$spec)
        text(0.5, 55, "N=256, Yule-Walker AR(8)", adj = 1, cex = 0.8)
        axis(2, at = c(-40, -20, 0, 20, 40, 60))
        box()
        n <- 1024
        x <- x3
        ar.ld <- lev.dur(x, ord = 8)
        ar.specld <- spec.ar(ar.ld)
        ar.spec <- spec.ar(list(ar = ar4, var.pred = 1,order=4))
        f <- seq(0, 0.5, l = length(ar.specld))
        plot(ar.specld$freq, ar.specld$spec, type = "l", axes = F, xlab = "", ylab = "dB", ylim = range(-40, 60), lwd = 2)
        lines(ar.spec$freq, ar.spec$spec)
        text(0.5, 55, "N=1024, Yule-Walker AR(8)", adj = 1, cex = 0.8)
        axis(2, at = c(-40, -20, 0, 20, 40, 60))
        axis(1, at = seq(0, 0.5, by = 0.1))
        box()
}