#' Function poinar
#'
#' Fit in inar model to a univariate time series by yule-walker method
#'
#' @param x a numeric vector or time series.
#' @param order.max a one dimensional integer vestor giving the order of the model to fit. This value corresponds the INAR order.
#' @param series name for the series.
#'
#'@return Resultados
#'
#'@references
#'
#'Du, J.G. and Li,Y.(1991):
#'The integer-valued autorregressive (INAR(p)) model.
#'\emph{Journal of time series analysis} \bold{12}, 129--142.
#'
#'Freeland RK. Statistical analysis of discrete time series with applications to the analysis of workers compensation
#'claims data [unpublished doctoral dissertation]. Vancouver (Canada): University of British Columbia; 1998.
#'
#'@examples
#'data(claims)
#'mean(claims[,5])
#'var(claims[,5])
#'var(claims[,5])/mean(claims[,5]) # dispersion index
#'acf(claims[,5])
#'pacf(claims[,5])
#'poinar(claims[,5], 1)
#'
#' @export
poinar <-
function(x, order.max,series=NULL)
{
if (is.null(series))
series <- deparse(substitute(x))
xfreq <- frequency(x)
n <- length(x)
r <- acf(x, plot = FALSE)$acf[2:(order.max+1)]
R <- diag(order.max)
for(i in 1:order.max){
for(j in 1:order.max){
if(i!=j){
R[i,j] <- r[abs(i-j)]
}
}
}
residuals <- function(x,coef,p,lambda)
{
x<- x
e <- NULL
for(t in (p+1):length(x) )
{
e[t] <- x[t] - ( sum(coef*x[(t-1):(t-p)]) + lambda)
}
return(e)
}
coef <- round(solve(R, r), 4)
xbar <- mean(x) #mean of the serie
mu.e <- var.error <- xbar*(1-sum(coef)) #mean and variance of the error
mu.x <- mu.e/(1-sum(coef))
sum.var <- sum( coef*(1-coef))
Vp <- var.error + mu.x*sum.var
fitted <- poinar.sim(length(x), order.max = length(coef), alpha = coef,lambda = mu.e, n.start=200)
resid <- residuals(x,coef,order.max,mu.e)
rms <- sqrt(mean(resid*resid))
AICc. <- n*log(Vp) + n*((1+order.max/n)/(1-(order.max+2)/n))
AIC. <- -n*log(Vp) + 2*order.max
BIC. <- n*log(Vp) + (order.max/n)*log(n)
res <- list(order = order.max,
coef = coef,
lambda= mu.e,
x.mean = xbar,
rms = rms,
fitted.values <- fitted,
bic= BIC.,
aicc=AICc.,
aic = AIC.,
n.used = n,
order.max = order.max,
resid = na.remove(resid),
method = "yule-walker",
series = series,
frequency = xfreq,
call = match.call())
class(res) <- "inar"
res
}
coef.inar <-
function(object, ...)
{
if(!inherits(object, "inar"))
stop("method is only for arma objects")
return(object$coef)
}
residuals.inar <-
function(object, ...)
{
if(!inherits(object, "inar"))
stop("method is only for arma objects")
return(object$resid)
}
fitted.inar <-
function(object, ...)
{
if(!inherits(object, "inar"))
stop("method is only for arma objects")
return(object$fitted.values)
}
print.inar <-
function(x, digits = max(3, getOption("digits") - 3), ...)
{
if(!inherits(x, "inar"))
stop("method is only for inar objects")
cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
rq <- structure(quantile(x$resid),names=c("Min","1Q","Median","3Q","Max"))
print(rq,digits = digits,...)
cat("Coefficient(s):\n")
print.default(format(coef(x), digits = digits),
print.gap = 2, quote = FALSE)
cat("\nFit:\n")
cat("RMS=", format(x$rms, digits = digits),
", lambda estimated as = ", format(round(x$lambda, 2)),
", AIC = ", format(round(x$aic, 2)), "\n", sep = "")
cat("\n")
invisible(x)
}
plot.inar <-
function(x, ask = interactive(), ...)
{
if(!inherits(x, "inar"))
stop("method is only for inar objects")
op <- par()
par(ask = ask, mfrow = c(2, 1))
data <- eval.parent(parse(text = x$series))
if(any(is.na(data))) stop(paste("NAs in", x$series))
plot(data, main = x$series, ylab = "Series")
plot(x$resid, main = "Residuals", ylab = "Series")
acf(data, main = paste("ACF of", x$series))
acf(x$resid, main = "ACF of Residuals",
na.action = na.remove)
pacf(data, main = paste("PACF of", x$series))
pacf(x$resid, main = "PACF of Residuals",
na.action = na.remove)
par(ask = op$ask, mfrow = op$mfrow)
invisible(x)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.