# R/DLLoglikelihood.R In ltsa: Linear Time Series Analysis

````DLLoglikelihood` <-
function (r, z, useC = TRUE)
{
EPS <- .Machine\$double.eps # 1+EPS==1, machine epsilon
if (useC) {
if (length(r) != length(z))
stop("Error: arguments must be same length")
stanQ <- 1
out <- .C("DLResid", as.double(z), res = as.double(numeric(length(z))),
as.integer(length(z)), as.double(r), LogL = as.double(1), EPS, as.integer(stanQ),
fault = as.integer(1), PACKAGE = "ltsa")
fault <- out\$fault
if (fault == 1)
stop("error: sequence not p.d.")
LogL <- out\$LogL
}
else {
n <- length(z)
if (n != length(r))
stop("arguments have unequal length")
error <- numeric(n)
sigmasq <- numeric(n)
error[1] <- z[1]
sigmasq[1] <- r[1]
phi <- r[2]/r[1]
error[2] <- z[2] - phi * z[1]
sigmasqkm1 <- r[1] * (1 - phi^2)
sigmasq[2] <- sigmasqkm1
logg <- log(r[1]) + log(sigmasqkm1)
for (k in 2:(n - 1)) {
if (sigmasqkm1 < 0 || abs(sigmasqkm1) < EPS)
stop("r is not a p.d. sequence")
phikk <- (r[k + 1] - phi %*% rev(r[2:k]))/sigmasqkm1
sigmasqk <- sigmasqkm1 * (1 - phikk^2)
phinew <- phi - phikk * rev(phi)
phi <- c(phinew, phikk)
sigmasqkm1 <- sigmasqk
logg <- logg + log(sigmasqk)
error[k + 1] <- z[k + 1] - crossprod(phi, rev(z[1:k]))
sigmasq[k + 1] <- sigmasqk
}
S <- sum((error * error)/sigmasq)
LogL <- c(-n/2 * log(S/n) - logg/2)
}
LogL
}
```

## Try the ltsa package in your browser

Any scripts or data that you put into this service are public.

ltsa documentation built on May 2, 2019, 4:01 a.m.