GetL {bcpa} | R Documentation |
Obtain likelihood of gappy standardized Gaussian time series "x" sampled at times "t" given parameter "rho" (autocorrelation). Alternatively computes the characteristic time scale "tau".
x |
Time series |
t |
Sampling times |
rho |
Auto-correlation |
tau |
logical: Whether or not to compute characteristic time scale instead of rho. |
Returns the log-likelihood of the data.
Eliezer Gurarie
Core function of BCPA, used directly in GetRho
# 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))