## -------------------------------
## Example: Zero-Inflated Beta-Binomial AR(1) model
## -------------------------------

library(gctsc)

## --- Parameter setup ---
n    <- 500
size <- 24
phi  <- 0.5
tau  <- c(phi)
arma_order <- c(1, 0)

## True parameters
beta_mu  <- 1.2                 # logit-scale parameter for prob
prob     <- plogis(beta_mu)     # simulation-scale success probability
rho      <- 0.1                 # BB overdispersion (ICC)
pi0      <- 0.2                 # constant zero-inflation prob

## --- Simulate ZIBB time series ---
set.seed(7)
sim_data <- sim_zibb(
  prob       = rep(prob, n),
  rho        = rho,
  pi0        = rep(pi0, n),
  size       = size,
  tau        = tau,
  arma_order = arma_order,
  family = "gaussian",
  nsim       = n
)
y <- sim_data$y

X <- list(mu = matrix(1, nrow = n), pi0 = matrix(1, nrow = n))

## --- Compute truncation bounds ---
lambda <- c(qlogis(prob), qlogis(rho), qlogis(pi0))

marg <- zibb.marg(link = "logit", size= size)
ab <- marg$bounds(y, X, lambda, family ="gaussian")

## --- Likelihood approximation ---
llk_tmet <- pmvn(lower = ab[,1], upper = ab[,2],
                 tau = tau, od = arma_order, 
                 pm = 30, QMC = TRUE,  method = "TMET")

llk_ghk  <- pmvn( lower = ab[,1], upper = ab[,2],
                  tau = tau, od = arma_order,
                  QMC = TRUE,  method = "GHK")

c(TMET = llk_tmet, GHK = llk_ghk)

## --- Fit ZIBB copula model ---
fit_zibb <- gctsc(
  formula  = list(mu = y ~ 1, pi0 = ~ 1), data = data.frame(y),
  marginal = zibb.marg(link = "logit", size = size),
  cormat   = arma.cormat(p = 1, q = 0),
  method   = "GHK", family = "gaussian",
  options  = gctsc.opts(seed = 7, M = 1000)
)

summary(fit_zibb)
plot(fit_zibb)
predict(fit_zibb)

## -------------------------------
## Example: ZIBB AR(1) with seasonal zero-inflation π₀(t)
## -------------------------------

library(gctsc)

## --- Parameter setup ---
n    <- 500
size <- 24
phi  <- 0.85
tau  <- c(phi)
arma_order <- c(1, 0)

prob <- 0.2        
rho  <- 0.16               

## Seasonal covariates for π₀(t)
X_pi <- cbind(
  1,
  sin(2*pi*(1:n)/12),
  cos(2*pi*(1:n)/12)
)
colnames(X_pi) <- c("int", "sin", "cos")

beta_pi <- c(1, 2, 5)
pi0 <- plogis(X_pi %*% beta_pi)

## --- Simulate ZIBB data ---
set.seed(1)
sim_data <- sim_zibb(
  prob       = rep(prob, n),
  rho        = rho,
  pi0        = pi0,
  size       = size,
  tau        = tau,
  arma_order = arma_order,
  family = "gaussian",
  nsim       = n
)
y <- sim_data$y


## --- Compute truncation bounds ---
lambda <- c(qlogis(prob), qlogis(rho), beta_pi)

marg <- zibb.marg(link = "logit", size= size)
ab <- marg$bounds(y, X, lambda, family ="gaussian")

## --- Likelihood approximation ---
llk_tmet <- pmvn(lower = ab[,1], upper = ab[,2],
                 tau = tau, od = arma_order, 
                 pm = 30, QMC = TRUE,  method  = "TMET")

llk_ghk  <- pmvn( lower = ab[,1], upper = ab[,2],
                  tau = tau, od = arma_order,
                  QMC = TRUE,  method = "GHK")

c(TMET = llk_tmet, GHK = llk_ghk)

## --- Fit using formula interface ---
df <- data.frame(y = y, X_pi)

fit_zibb_seasonal <- gctsc(
  formula = list(mu = y ~ 1,
                 pi0 = ~ sin + cos),
  data     = df,
  marginal = zibb.marg(link = "logit", size = size),
  cormat   = arma.cormat(p = 1, q = 0),
  method   = "TMET",family = "gaussian",
  options  = gctsc.opts(seed = 1, M = 1000)
)

summary(fit_zibb_seasonal)
plot(fit_zibb_seasonal)
predict(fit_zibb_seasonal,newdata = df[500, ])


## -------------------------------
## Example: ZIBB AR(1) with covariates in μ(t) and π₀(t)
## -------------------------------


library(gctsc)

## --- Parameter setup ---
n    <- 1000
size <- 10
phi  <- 0.5
tau  <- c(phi)
arma_order <- c(1, 0)
rho  <- 0.3

## --- Covariates for μ(t) ---
X_mu <- cbind(
  1,
  sin(2*pi*(1:n)/365),
  cos(2*pi*(1:n)/365)
)
colnames(X_mu) <- c("int", "sin", "cos")

beta_mu <- c(0.1,0.4,0.5)
prob <- plogis(X_mu %*% beta_mu)

## --- Covariates for π₀(t) ---
day_of_year <- rep(1:365, length.out = n)
season <- factor(ifelse(day_of_year < 100, "Winter",
                        ifelse(day_of_year < 180, "Spring",
                               ifelse(day_of_year < 270, "Summer", "Fall"))))
X_pi <- model.matrix(~ season)

alpha_pi <- c(0.2, 0.3, 0.5, 0.7)
pi0 <- plogis(X_pi %*% alpha_pi)

## --- Simulate ZIBB data ---
set.seed(1)
sim_data <- sim_zibb(
  prob       = prob,
  rho        = rho,
  pi0        = pi0,
  size       = size,
  tau        = tau,
  arma_order = arma_order,
  nsim       = n,
  family = "gaussian"
)
y <- sim_data$y

X <- list(mu = as.matrix(X_mu), pi0 = as.matrix(X_pi))

## --- Compute truncation bounds ---
lambda <- c(beta_mu, qlogis(rho), alpha_pi)

marg <- zibb.marg(link = "logit", size= size)
ab <- marg$bounds(y, X, lambda, family ="gaussian")

## --- Likelihood approximation ---
llk_tmet <- pmvn(lower = ab[,1], upper = ab[,2],
                      tau = tau, od = arma_order, 
                      pm = 30, QMC = TRUE, method ="TMET")

llk_ghk  <- pmvn( lower = ab[,1], upper = ab[,2],
                      tau = tau, od = arma_order,
                      QMC = TRUE, method  ="GHK")

c(TMET = llk_tmet, GHK = llk_ghk)


## --- Fit ZIBB copula model using formula interface ---
data_df <- data.frame(y = y, X_mu, X_pi)

fit_zibb_cov <- gctsc(
  formula = list(mu  = y ~ sin + cos,
                 pi0 = ~ seasonSpring + seasonSummer + seasonWinter),
  data     = data_df,
  marginal = zibb.marg(link = "logit", size = size),
  cormat   = arma.cormat(p = 1, q = 0),
  method   = "TMET",family = "gaussian",
  options  = gctsc.opts(seed = 1, M =1000)
)

summary(fit_zibb_cov)
plot(fit_zibb_cov)

predict(fit_zibb_cov, newdata = data_df[500, ])

