Gaussian and Student t Copula Time Series Models for Count Data

The gctsc package provides fast and scalable tools for fitting Gaussian and Student t copula models for count time series, supporting a wide range of discrete marginals:

These models are constructed from a latent ARMA process combined with a discrete marginal distribution. Likelihood evaluation requires computing multivariate probabilities under Gaussian or Student–t dependence structures.

The package provides several likelihood approximation methods:

In addition, the package offers:

A key feature of the package is the TMET method, which exploits the conditional structure of ARMA processes to compute likelihood approximations efficiently even for long time series. The package also includes the classical GHK simulator and the Continuous Extension approximation for comparison.

This vignette illustrates the main functionality of the gctsc package, including simulation, likelihood approximation, model fitting, diagnostics, and prediction for copula-based count time series models.

Installation

The package can be installed from CRAN using

install.packages("gctsc")

The development version is available at

remotes::install_github("QNNHU/gctsc")

Package overview

The package provides tools for simulation, likelihood evaluation, estimation, diagnostics, and prediction for copula-based count time series models with Gaussian or Student–t copula dependence.

The package implements several likelihood approximation methods, including Time Series Minimax Exponential Tilting (TMET), the classical Geweke–Hajivassiliou–Keane (GHK) simulator, and the Continuous Extension (CE) approximation.

Main functions

Function Description
sim_*() Simulate copula count time series with latent ARMA dependence
pmvn_*() / pmvt_*() Compute multivariate rectangle probabilities for Gaussian and Student–t copulas
gctsc() Fit Gaussian or Student–t copula count time series models
predict() Compute one-step-ahead predictive distributions and moments
plot() Produce residual diagnostic plots
summary() Display parameter estimates and model fit statistics

The following sections illustrate the main functionality of the package, including simulation, likelihood approximation, model estimation, diagnostics, and prediction.

Example 1: Simulating Poisson AR(1) counts under a Gaussian copula

n  <- 300
mu <- 10
phi <- 0.2
arma_order <- c(1, 0)
tau <- c(phi)

sim_data <- sim_poisson(
  mu = mu,
  tau = tau,
  arma_order = arma_order,
  nsim = n,
  family = "gaussian",
  seed = 7
)

y <- sim_data$y
plot(y, type = "l", main = "Simulated Poisson AR(1) Counts: Gaussian Copula")

Example 2: Simulating Poisson AR(1) counts under a Student–t copula

n  <- 300
mu <- 10
phi <- 0.2
arma_order <- c(1, 0)
tau <- c(phi)

sim_data <- sim_poisson(
  mu = mu,
  tau = tau,
  arma_order = arma_order,
  nsim = n,
  family  = "t",
  df = 5,
  seed = 7
)

y <- sim_data$y
plot(y, type = "l", main = "Simulated Poisson AR(1) Counts: t Copula")

Example 3: Simulating Zero–Inflated Beta–Binomial (ZIBB) AR(1) Count Time Series

n <- 300
phi <- 0.5
tau <- c(phi)

beta_0 <- 1.2
prob   <- plogis(beta_0)
rho    <- 0.1
pi0    <- 0.8
size   <- 24

set.seed(7)
sim_data <- sim_zibb(prob = prob * rep(1,n),
                     rho = rho, pi0 = pi0,
                     size = size, tau = tau,
                     arma_order = c(1,0),
                     family = "gaussian",
                     nsim = n)
y <- sim_data$y

plot(y, type = "l", main = "Simulated ZIBB AR(1) Counts: gaussian Copula")

Example 4: Likelihood approximation (TMET vs GHK vs CE)

n <- 300
mu <- 10
phi <- 0.5
theta <- 0.2
arma_order <- c(1,1)
tau <- c(phi, theta)

# Simulate Poisson count data
sim_data <- sim_poisson(mu, tau, arma_order, nsim = n, seed = 1)
y <- sim_data$y

# Compute bounds for truncated MVN using marginal object
marg <- poisson.marg()
lower <- qnorm(ppois(y - 1, lambda = mu))
upper <- qnorm(ppois(y, lambda = mu))

# Likelihood approximation
llk_tmet <- pmvn_tmet(lower, upper, tau, od = arma_order, pm = 30, QMC = TRUE)
llk_ghk  <- pmvn_ghk(lower, upper, tau, od = arma_order, QMC = TRUE)
llk_ce  <- pmvn_ce(lower, upper, tau, od = arma_order)

c(TMET = llk_tmet, GHK = llk_ghk, CE = llk_ce)
#>      TMET       GHK        CE 
#> -704.8721 -704.8980 -704.8964

Example 4: Fitting a Gaussian Copula Count Model


n <- 300
mu <- 10
phi <- 0.5
theta <- 0.2
arma_order <- c(1,1)
tau <- c(phi, theta)

# Simulate Poisson count data
sim_data <- sim_poisson(mu, tau, arma_order, nsim = n, seed = 1, family ="gaussian")
y <- sim_data$y

fit <- gctsc(
  formula  = y ~ 1,
  marginal = poisson.marg(),
  cormat   = arma.cormat(p = 1, q = 1),
  method   = "CE", family ="gaussian"
)

summary(fit)
#> 
#> --- Summary of Copula Time Series Model ---
#> Call:
#>  gctsc(formula = y ~ 1, marginal = poisson.marg(), cormat = arma.cormat(p = 1,      q = 1), method = "CE", family = "gaussian") 
#> 
#> Marginal Model Coefficients:
#>                Estimate Std. Error z value Pr(>|z|)    
#> mu_(Intercept)  10.1915     0.3254   31.32   <2e-16 ***
#> 
#> Copula (Dependence) Coefficients:
#>     Estimate Std. Error z value Pr(>|z|)    
#> ar1  0.45000    0.07367   6.108 1.01e-09 ***
#> ma1  0.29137    0.09428   3.090    0.002 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Model Fit Statistics:
#>   Log-likelihood: -704.2735 
#>   AIC:             1414.5469 
#>   BIC:             1425.6583

Example 5: Residual diagnostics

oldpar <- par(no.readonly = TRUE)

par(mfrow = c(2,3))
plot(fit)


par(oldpar)

Example 6: One-step-ahead prediction

prediction <- predict(
  fit,
  y_obs   = 10
)

prediction
#> $mean
#> [1] 9.409088
#> 
#> $median
#> [1] 9
#> 
#> $mode
#> [1] 9
#> 
#> $variance
#> [1] 5.827867
#> 
#> $p_y
#>  [1] 6.686625e-07 2.714838e-05 3.567118e-04 2.416803e-03 1.022103e-02
#>  [6] 2.996922e-02 6.511687e-02 1.097833e-01 1.485445e-01 1.655231e-01
#> [11] 1.550210e-01 1.240559e-01 8.599216e-02 5.222636e-02 2.806394e-02
#> [16] 1.345521e-02 5.798285e-03 2.260343e-03 8.016669e-04 2.599981e-04
#> [21] 7.746185e-05 2.128819e-05 5.416863e-06 1.280548e-06 2.821226e-07
#> [26] 5.809223e-08 1.120942e-08 2.031859e-09 3.467630e-10 5.583605e-11
#> [31] 8.499461e-12 1.225362e-12 1.534571e-13 8.980257e-14 8.980257e-14
#> 
#> $lower
#> [1] 5
#> 
#> $upper
#> [1] 14
#> 
#> $CRPS
#> [1] 0.625926
#> 
#> $LOGS
#> [1] 1.864194

Example 7: Fitting a t Copula Count Model

n <- 100
mu <- 10
phi <- 0.5
theta <- 0.2
arma_order <- c(1,1)
tau <- c(phi, theta)

# Simulate Poisson count data
sim_data <- sim_poisson(mu, tau, arma_order, nsim = n, family = "t", df= 15, seed = 1)
y <- sim_data$y

fit <- gctsc(
  formula  = y ~ 1,
  marginal = poisson.marg(),
  cormat   = arma.cormat(p = 1, q = 1),
  method   = "CE", family ="t", df= 15
)

summary(fit)
#> 
#> --- Summary of Copula Time Series Model ---
#> Call:
#>  gctsc(formula = y ~ 1, marginal = poisson.marg(), cormat = arma.cormat(p = 1,      q = 1), method = "CE", family = "t", df = 15) 
#> 
#> Marginal Model Coefficients:
#>                Estimate Std. Error z value Pr(>|z|)    
#> mu_(Intercept)   9.8799     0.5222   18.92   <2e-16 ***
#> 
#> Copula (Dependence) Coefficients:
#>     Estimate Std. Error z value Pr(>|z|)    
#> ar1   0.4728     0.1409   3.355 0.000793 ***
#> ma1   0.3067     0.1628   1.884 0.059518 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Model Fit Statistics:
#>   Log-likelihood: -217.1995 
#>   AIC:             440.3991 
#>   BIC:             448.2146

Example 7: Real Data Analysis — Campylobacter Time Series

## Load weekly Campylobacter incidence data
data("campyl", package = "gctsc")
y <- campyl[,"X1"]
n <- length(y)

## Plot the time series
ts_y <- ts(y, start = c(2001, 1), frequency = 52)
plot(ts_y, main = "Weekly Campylobacter Incidence",
     ylab = "Cases", xlab = "Year")


## ---------------------------------------------------------------
## Construct seasonal covariates
## ---------------------------------------------------------------
## Seasonal structure appears to have yearly periodicity,
## so we include sine/cosine terms with period = 52 weeks.

time <- 1:n
X <- data.frame(
  intercept = 1,
  sin52 = sin(2 * pi * time / 52),
  cos52 = cos(2 * pi * time / 52)
)

## Combine response and covariates
data_df <- data.frame(Y = y, X)

## Use the first 800 observations for model fitting
train_end <- 1000
data_train <- data_df[1:train_end, ]

## ---------------------------------------------------------------
## Fit a Negative Binomial Gaussian Copula model
## ---------------------------------------------------------------
## We use:
##   - Negative Binomial marginal (log link)
##   - ARMA(1,1) latent correlation structure
##   - GHK likelihood approximation
##
## The model is:
##     E[Y_t] = exp(β0 + β1 sin + β2 cos)
##
## Covariates enter only through the marginal mean.

fit <- gctsc(
  formula  = Y ~ sin52 + cos52,
  data     = data_train,
  marginal = negbin.marg(link = "log"),
  cormat   = arma.cormat(p = 1, q = 1),
  method   = "CE", family ="gaussian",
  options  = gctsc.opts(seed = 1)   # fixed seed for reproducibility
)

summary(fit)
#> 
#> --- Summary of Copula Time Series Model ---
#> Call:
#>  gctsc(formula = Y ~ sin52 + cos52, data = data_train, marginal = negbin.marg(link = "log"),      cormat = arma.cormat(p = 1, q = 1), method = "CE", family = "gaussian",      options = gctsc.opts(seed = 1)) 
#> 
#> Marginal Model Coefficients:
#>                Estimate Std. Error z value Pr(>|z|)    
#> mu_(Intercept)  1.57385    0.07968  19.753  < 2e-16 ***
#> mu_sin52       -0.31857    0.02911 -10.945  < 2e-16 ***
#> mu_cos52       -0.30412    0.02860 -10.634  < 2e-16 ***
#> dispersion      0.13623    0.03234   4.212 2.53e-05 ***
#> 
#> Copula (Dependence) Coefficients:
#>      Estimate Std. Error z value Pr(>|z|)    
#> ar1  0.988033   0.006538  151.13   <2e-16 ***
#> ma1 -0.914750   0.020137  -45.42   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Model Fit Statistics:
#>   Log-likelihood: -2318.2763 
#>   AIC:             4648.5526 
#>   BIC:             4677.9992

## ---------------------------------------------------------------
## Residual diagnostics
## ---------------------------------------------------------------
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))   # restore user settings

par(mfrow=c(2,3))
plot(fit)



## ---------------------------------------------------------------
## One-step-ahead prediction
## ---------------------------------------------------------------
## Predict Y_{801} using fitted model
new_obs_index <- train_end + 1
pred <- predict(
  fit,
  X_test = data_df[new_obs_index, ],
  y_obs  = data_df[new_obs_index, "Y"]
)

pred
#> $mean
#> [1] 4.022336
#> 
#> $median
#> [1] 4
#> 
#> $mode
#> [1] 3
#> 
#> $variance
#> [1] 4.855465
#> 
#> $p_y
#>  [1] 2.133817e-02 8.592509e-02 1.558846e-01 1.893890e-01 1.785033e-01
#>  [6] 1.407713e-01 9.716988e-02 6.045921e-02 3.460629e-02 1.849376e-02
#> [11] 9.330130e-03 4.481838e-03 2.063762e-03 9.159080e-04 3.935089e-04
#> [16] 1.642699e-04 6.683347e-05 2.656987e-05 1.034443e-05 3.951611e-06
#> [21] 1.483580e-06 5.482073e-07 1.996306e-07 7.172078e-08 2.544668e-08
#> [26] 8.924245e-09 3.009186e-09 3.009186e-09 3.009186e-09 3.009186e-09
#> 
#> $lower
#> [1] 1
#> 
#> $upper
#> [1] 9
#> 
#> $CRPS
#> [1] 1.066654
#> 
#> $LOGS
#> [1] 1.858639

Further examples (Negative Binomial, Binomial,Beta-Binomial, ZIP, ZIB, ZIBB, including cases with covariates) and another real data analysis are available in the inst/examples/ directory of the package source.