# R/tsglm.interv_decompose.R In tscount: Analysis of Count Time Series

#### Defines functions tsglm.interv_decompose

```tsglm.interv_decompose <- function(ts, model, xreg, link, paramvec, isolate){
#Decomposes a time series in the effect of one or several covariates and an INGARCH model with the remaining covariates
#isolate: Integer vector. Giving the indices of the covariates whose effect should be separated from the time series.
##############################

##############
#Checks and preparations:
n <- length(ts)
p <- length(model\$past_obs)
P <- seq(along=numeric(p)) #sequence 1:p if p>0 and NULL otherwise
p_max <- max(model\$past_obs, 0)
q <- length(model\$past_mean)
Q <- seq(along=numeric(q)) #sequence 1:q if q>0 and NULL otherwise
q_max <- max(model\$past_mean, 0)
Q_max <- seq(along=numeric(q_max))
r <- max(ncol(xreg), 0)
R <- seq(along=numeric(r)) #sequence 1:r if r>0 and NULL otherwise
param <- list( #transform parameter vector to a list
intercept=paramvec[1],
past_obs=paramvec[1+P],
past_mean=paramvec[1+p+Q],
xreg=paramvec[1+p+q+R]
)
stopifnot(
r>0,
1<=isolate & isolate<=r
)
elim <- 1:r %in% isolate
##############

##############
#Initialisation:
#The code for this initialisation is taken from the function tsglm.condmean. Do not modify it here without modifying it there.
#Initialisation by stationary solution (and its partial derivatives):
denom <- (1-sum(param\$past_obs)-sum(param\$past_mean))[[1]]
nu_stationary <- (param\$intercept/denom)[[1]]
nu <- c(rep(nu_stationary, q_max), numeric(n))
mu <- numeric(q_max+n)
contamination <- integer(p_max+n)
cleaned <- y
X <- matrix(0, nrow=q_max+n, ncol=r) #regressors are initalised by zero
X[q_max+(1:n), ] <- xreg
##############

##############
#Recursion:
for(t in 1:n){
nu[t+q_max] <- param\$intercept + sum(param\$past_obs*y[(t-model\$past_obs)+p_max]) + sum(param\$past_mean*nu[(t-model\$past_mean)+q_max]) + sum(param\$xreg*X[t+q_max, ]) - sum(param\$past_mean*colSums(model\$external*param\$xreg*t(X[(t-model\$past_mean)+q_max, , drop=FALSE])))
mu[t+q_max] <- sum(param\$past_obs*contamination[(t-model\$past_obs)+p_max]) + sum(param\$past_mean*mu[(t-model\$past_mean)+q_max]) + sum(param\$xreg[elim]*X[t+q_max, elim]) - sum(param\$past_mean*colSums(model\$external[elim]*param\$xreg[elim]*t(X[(t-model\$past_mean)+q_max, elim, drop=FALSE])))
contamination[t+p_max] <- min(round(mu[t+q_max]/nu[t+q_max]*y[t+p_max]),y[t+p_max]) #contamination may not be larger than the original observation
}
cleaned <- y-contamination
}
for(t in 1:n){
mu[t+q_max] <- sum(param\$past_obs*log(1+contamination[(t-model\$past_obs)+p_max]/(cleaned[(t-model\$past_obs)+p_max]+1))) + sum(param\$past_mean*mu[(t-model\$past_mean)+q_max]) + sum(param\$xreg[elim]*X[t+q_max, elim]) - sum(param\$past_mean*colSums(model\$external[elim]*param\$xreg[elim]*t(X[(t-model\$past_mean)+q_max, elim, drop=FALSE])))
cleaned[t+p_max] <- round(y[t+p_max]/exp(mu[t+q_max]))
contamination[t+p_max] <- y[t+p_max]-cleaned[t+p_max]
}
}

result <- list(ts=ts, cleaned=cleaned[p_max+(1:n)], contamination=contamination[p_max+(1:n)])
return(result)
}
```

## Try the tscount package in your browser

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

tscount documentation built on Nov. 25, 2017, 3:01 a.m.