# chex44.s                         J F Monahan, June 2001
 # likelihood computation for ARMA model
 phi <- .9
 th  <- .5
 gam <- vector("numeric",7)
 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:7) gam[j] <- phi*gam[j-1]
 gam
 d <- 1 + abs(matrix(1:4,4,4) - t(matrix(1:4,4,4)))
 A <- matrix(gam[d],4,4)
 "covariance matrix"
 A

 "Cholesky factor"
 Lt <- chol(A)
 Lt
 det <-prod(diag(Lt))
 "determinant"
 det*det
 z <- c(-2, 1, 0, 1)
 one <- rep(1,4)
 LiBz <- solve(t(Lt),z)
 LiBone <- solve(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
 rm(list=ls())