# 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()