# chex45.r J F Monahan, May 2010 # # demonstration of Levinson-Durbin algorithm # likelihood computation for ARMA(1,1) model # # set parameter values phi <- .9 th <- .5 n <- 4 # table covariance function -- gamma(0) stored in gam[1] gam <- rep(0,n+1) gam[1] <- (1+th*th-2*phi*th)/(1-phi*phi) gam[2] <- (1-phi*th)*(phi-th)/(1-phi*phi) for(j in 3:(n+1)) gam[j] <- phi*gam[j-1] "covariance function" gam # Toeplitz covariance matrix d <- 1 + abs(matrix(1:n,n,n) - t(matrix(1:n,n,n))) A <- matrix(gam[d],n,n) "covariance matrix" A # make correlations r <- gam[1+(1:n)]/gam[1] "correlations" r "right hand side" y <- c(-2,1,0,1) y # start Levinson-Durbin d <- 1 b <- -r[1] x <- y[1] det <- 1 # set up for Cholesky factor LAinv <- diag(rep(1,n)) for(k in(2:n)) { # first the scalars d <- 1 + sum(r[1:(k-1)]*b[1:(k-1)]) LAinv[((n-k+1):n),n-k+1] <- c(1,b)/sqrt(d) e <- r[k] + sum(r[rev(1:(k-1))]*b[1:(k-1)]) f <- y[k] - sum(r[rev(1:(k-1))]*x[1:(k-1)]) print("d,e,f") print(c(d,e,f)) # then the vectors x <- c(x+(f/d)*b[rev(1:(k-1))],f/d) print(x) b <- c(b-(e/d)*b[rev(1:(k-1))],-e/d) print(b) det <- det*d } "solution" x/gam[1] "determinant" det*(gam[1]^n) "Cholesky factor of inverse" LAinv "Is it identity? (except for gamma(0))" t(LAinv) %*% A %*% LAinv "Cholesky factor" Lt <- chol(A) Lt det <-prod(diag(Lt)) "determinant" det*det z <- c(-2, 1, 0, 1) one <- rep(1,4) LiBz <- forwardsolve(t(Lt),z) "solution" backsolve(Lt,LiBz) LiBone <- forwardsolve(t(Lt),one) zAiz <- sum( LiBz * LiBz ) zAione <- sum( LiBz * LiBone ) oneAione <- sum( LiBone * LiBone ) muhat <- zAione / oneAione "muhat" muhat "quadratic form" qf <- zAiz - zAione*muhat qf # clean up rm(list=ls()) q()