GetL: Obtain likelihood estimates of gappy Gaussian time series

In bcpa: Behavioral Change Point Analysis of Animal Movement

Obtain likelihood estimates of gappy Gaussian time series

Description

Obtain likelihood of gappy standardized Gaussian time series "x" sampled at times "t" given parameter "rho" (autocorrelation). Alternatively computes the characteristic time scale "tau".

Arguments

x	Time series
t	Sampling times
rho	Auto-correlation
tau	logical: Whether or not to compute characteristic time scale instead of rho.

Value

Returns the log-likelihood of the data.

Author(s)

Eliezer Gurarie

See Also

Core function of BCPA, used directly in GetRho

Examples

# simulate autocorrelated time series
  rho.true <- 0.8
  x.full <- arima.sim(1000, model=list(ar = rho.true))
  t.full <- 1:1000
  
# subsample time series
  keep <- sort(sample(1:1000, 200))
  x <- x.full[keep]
  t <- t.full[keep]
  plot(t,x, type="l")
  
# Obtain MLE of rho
  rhos <- seq(0,.99,.01)
  L <- sapply(rhos, function(r) GetL(x, t, r))
  rho.hat <- rhos[which.max(L)]
  plot(rhos, L, type = "l")
  abline(v = c(rho.true, rho.hat), lty=3:2, lwd=2)
  legend("bottomleft", legend=c("true value","MLE"), lty=3:2, lwd=2, 
         title = expression(rho))
         
# Why tau is better
  tau.true <- -1/log(rho.true)
  taus <- seq(1,10,.1)
  L <- sapply(taus, function(r) GetL(x, t, r, tau = TRUE))
  tau.hat <- taus[which.max(L)]

  plot(taus, L, type = "l")
  abline(v = c(tau.true, tau.hat), lty=3:2, lwd=2)
  legend("bottomleft", legend=c("true value","MLE"), lty=3:2, lwd=2, 
         title = expression(tau))

bcpa documentation built on May 30, 2022, 5:07 p.m.