# 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:

- Poisson
- Negative Binomial
- Binomial
- Beta Binomial
- Zero-Inflated Poisson (ZIP)
- Zero-Inflated Binomial (ZIB)
- Zero-Inflated Beta-Binomial (ZIBB)

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:

- TMET — Time Series Minimax Exponential Tilting
- GHK — Geweke–Hajivassiliou–Keane simulator
- CE — Continuous Extension approximation

In addition, the package offers:

- Simulation of count time series with ARMA latent structure
- Maximum likelihood estimation using the above likelihoods
- Residual diagnostic tools
- One-step-ahead prediction

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,
                         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, family = "gaussian",
                            seed = 1)
    y <- sim_data$y

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

    # Likelihood approximation
    llk_tmet <- pmvn(lower, upper, tau, od = arma_order, method ="TMET",
                     pm = 30, QMC = TRUE)
    llk_ghk  <- pmvn(lower, upper, tau, od = arma_order, method = "GHK",
                     QMC = TRUE)
    llk_ce  <-  pmvn(lower, upper, tau, od = arma_order, method ="CE")

    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 <- 1000
    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, data= data.frame(y),
      marginal = poisson.marg(link = "log"),
      cormat   = arma.cormat(p = 1, q = 1),
      method   = "CE", family ="gaussian"
    )

    summary(fit)
    #>
    #> --- Summary of Copula Time Series Model ---
    #> Call:
    #>  gctsc(formula = y ~ 1, data = data.frame(y), marginal = poisson.marg(link = "log"),      cormat = arma.cormat(p = 1, q = 1), method = "CE", family = "gaussian")
    #>
    #> Marginal Model Coefficients:
    #>              Estimate Std. Error z value Pr(>|z|)
    #> mu_intercept  2.29551    0.01899   120.9   <2e-16 ***
    #>
    #> Copula (Dependence) Coefficients:
    #>     Estimate Std. Error z value Pr(>|z|)
    #> ar1  0.55001    0.03535  15.557   <2e-16 ***
    #> ma1  0.12136    0.04818   2.519   0.0118 *
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Model Fit Statistics:
    #>   Log-likelihood: -2342.0674
    #>   AIC:             4690.1347
    #>   BIC:             4704.8580

### 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,
      newdata   = data.frame(y = 10)
    )

    prediction
    #> $mean
    #> [1] 8.372444
    #>
    #> $median
    #> [1] 8
    #>
    #> $mode
    #> [1] 8
    #>
    #> $variance
    #> [1] 5.500033
    #>
    #> $p_y
    #>  [1] 6.808993e-06 1.832288e-04 1.718823e-03 8.675555e-03 2.814953e-02
    #>  [6] 6.475330e-02 1.123686e-01 1.535625e-01 1.705589e-01 1.577209e-01
    #> [11] 1.237679e-01 8.370223e-02 4.941072e-02 2.573487e-02 1.193475e-02
    #> [16] 4.967343e-03 1.868283e-03 6.388474e-04 1.996720e-04 5.731676e-05
    #> [21] 1.517602e-05 3.720767e-06 8.476867e-07 1.800371e-07 3.575114e-08
    #> [26] 6.655685e-09 1.164536e-09 1.919431e-10 2.986621e-11 4.395828e-12
    #> [31] 6.131360e-13 8.117524e-14 1.507941e-14 1.507941e-14 1.507941e-14
    #> [36] 1.507941e-14
    #>
    #> $y_grid
    #>  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
    #> [26] 25 26 27 28 29 30 31 32 33 34 35
    #>
    #> $lower
    #> [1] 4
    #>
    #> $upper
    #> [1] 13
    #>
    #> $CRPS
    #> [1] 1.017134
    #>
    #> $LOGS
    #> [1] 2.089347
    #>
    #> $y_true
    #> [1] 10

### Example 7: Fitting a t Copula Count Model

    n <- 1000
    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, data = data.frame(y),
      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, data = data.frame(y), 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  2.29377    0.01794   127.9   <2e-16 ***
    #>
    #> Copula (Dependence) Coefficients:
    #>     Estimate Std. Error z value Pr(>|z|)
    #> ar1  0.49381    0.03825  12.909   <2e-16 ***
    #> ma1  0.12946    0.04897   2.644   0.0082 **
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Model Fit Statistics:
    #>   Log-likelihood: -2535.3621
    #>   AIC:             5076.7242
    #>   BIC:             5091.4474

### Example 7: Real Data Analysis — Campylobacter Time Series

    ## Load weekly Campylobacter incidence data
    data("campyl", package = "gctsc")
    y <- campyl[,"breisgau_hochschwarzwald_county"]
    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(
      sin52 = sin(2 * pi * time / 52),
      cos52 = cos(2 * pi * time / 52),
      sin52_2 = sin(2 * pi * 2 * time  / 52),
      cos52_2 = cos(2 * pi * 2 * time  / 52)
    )

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

    ## Use the first 1000 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
    ##   - CE likelihood approximation
    ##
    ## The marginal mean model is:
    ##   E[Y_t] = exp(beta0 + beta1 sin52_t + beta2 cos52_t
    ##                    + beta3 sin52_2_t + beta4 cos52_2_t)
    ##

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

    summary(fit)
    #>
    #> --- Summary of Copula Time Series Model ---
    #> Call:
    #>  gctsc(formula = Y ~ sin52 + cos52 + sin52_2 + cos52_2, data = data_train,      marginal = negbin.marg(link = "log"), cormat = arma.cormat(p = 1,          q = 1), method = "CE", family = "gaussian", options = gctsc.opts(seed = 1,          hessian = TRUE))
    #>
    #> Marginal Model Coefficients:
    #>              Estimate Std. Error z value Pr(>|z|)
    #> mu_intercept  1.57892    0.08239  19.164  < 2e-16 ***
    #> mu_sin52     -0.32650    0.02940 -11.106  < 2e-16 ***
    #> mu_cos52     -0.28217    0.02767 -10.197  < 2e-16 ***
    #> mu_sin52_2    0.04366    0.02407   1.814   0.0697 .
    #> mu_cos52_2    0.11412    0.02422   4.711 2.46e-06 ***
    #> dispersion    0.12735    0.03757   3.389   0.0007 ***
    #>
    #> Copula (Dependence) Coefficients:
    #>      Estimate Std. Error z value Pr(>|z|)
    #> ar1  0.989182   0.008509  116.25   <2e-16 ***
    #> ma1 -0.920228   0.020377  -45.16   <2e-16 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Model Fit Statistics:
    #>   Log-likelihood: -2305.9243
    #>   AIC:             4627.8485
    #>   BIC:             4667.1106

    ## ---------------------------------------------------------------
    ## 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_{1001} using fitted model
    new_obs_index <- train_end + 1
    pred <- predict(
      fit,
      newdata = data_df[new_obs_index, ]
    )

    pred
    #> $mean
    #> [1] 3.542668
    #>
    #> $median
    #> [1] 3
    #>
    #> $mode
    #> [1] 3
    #>
    #> $variance
    #> [1] 4.061811
    #>
    #> $p_y
    #>  [1] 3.132897e-02 1.146452e-01 1.885134e-01 2.068569e-01 1.754718e-01
    #>  [6] 1.241418e-01 7.665220e-02 4.255314e-02 2.168267e-02 1.029419e-02
    #> [11] 4.605515e-03 1.958686e-03 7.973530e-04 3.124275e-04 1.183681e-04
    #> [16] 4.352542e-05 1.558281e-05 5.446350e-06 1.862584e-06 6.245014e-07
    #> [21] 2.056372e-07 6.659951e-08 2.124285e-08 6.680858e-09 2.806416e-09
    #> [26] 2.806416e-09 2.806416e-09 2.806416e-09 2.806416e-09 2.806416e-09
    #> [31] 2.806416e-09
    #>
    #> $y_grid
    #>  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
    #> [26] 25 26 27 28 29 30
    #>
    #> $lower
    #> [1] 0
    #>
    #> $upper
    #> [1] 8
    #>
    #> $CRPS
    #> [1] 0.7898279
    #>
    #> $LOGS
    #> [1] 1.668586
    #>
    #> $y_true
    #> [1] 2

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.