data_generation.R

# Generalized Poisson AR(1) - GPAR(1) Simulation
#
# X_t = S_t(X_{t-1}) + e_t, t=1, 2 , ...
#
# e_t    ~ GP(q*mu, xi)
# S_t(.) ~ QB(alpha, xi/mu, .)
# X_t    ~ GP(mu, xi)
#
# Constraints:
#   0 < alpha < 1
#   mu > 0
#   xi >= 0

library(gp)

graphics.off()  # close all plots
rm(list=ls())	# clear the (global) environment variables
cat("\014") 	# clear the console (Ctrl+L)
set.seed(40)

# Print auxiliary function
printf <- function(...) {
    x = paste(sprintf(...),"\n")
    return(cat(x))
}

# Process parameters and outliers
alpha <- 0.5                    # autocorrelation parameter
mu    <- 2
xi    <- 0.5                    # dispersion parater
q     <- 1-alpha                # auxiliary variable
tau   <- c(7, 26, 60, 90, 91)   # outliers location
eta   <- 6                      # outliers magnitude
filename <- "s9"                # csv file to save the data

# Theoretical mean, variance and sd
mean_x <- mu/(1-xi)
var_x <- mu/(1-xi)^3
sd_x <- sqrt(var_x)

# Simulation period
T <- 120

# Quasi-binomial thinning pmf
dqb <- function(x, alpha, xi, n){
    q <- 1-alpha
    num <- alpha*q*choose(n,x) * (alpha + x*xi)^(x-1) * (q + (n-x)*xi)^(n-x-1)
    den <- (1 + n*xi)^(n-1)
    y <- num/den
    return(y)
}

# Quasi-binomial thinning cdf
pqb <- function(x, alpha, xi, n){
    x_ <- 0:x
    y <- sum(sapply(x_, dqb, alpha=alpha, xi=xi, n=n))
    return(y)
}

# Draw from quasi-binomial
rqb <- function(size, alpha, xi, n){
    x <- 0:n
    w <- sapply(x, dqb, alpha=alpha, xi=xi, n=n)
    y <- sample(x, size=size, replace=TRUE, prob=w)
    return(y)
}

# GPAR(1) simulation
x <- NULL
x[1] <- mean_x
for (t in 2:T){
    # Innovations
    e <- rgp(n=2, q*mu, xi, method="Branching")[1]

    # Quasi-binomial thinning operator
    s <- rqb(1, alpha, xi/mu, x[t-1])
    x[t] <- s + e
}

# Insert the outliers
y <- x
y[tau] <- y[tau] + eta

#####
# Theoretical vs simulated mean and variance
printf("Theoretical mean: %.2f", mean_x)
printf("Sample mean: %.2f", mean(x))
printf("Theoretical var: %.2f", var_x)
printf("Theoretical sd: %.2f", sd_x)
printf("Sample var: %.2f", var(x))
printf("mu + 2*sd: %.2f", mean_x + 2*sd_x)
printf("mu + 5*sd: %.2f", mean_x + 5*sd_x)

#####
# Plot
par(mar = c(2, 2, 0.2, 0.2))
plot(y, type='s')
points(tau, y[tau], col="red")
abline(h=mean_x, col="gray", lty=2)
abline(h=mean_x + 3*sd_x, col="gray", lty=3)
legend(x="topright", lty=c(2,3), col="gray", legend=c("E[X]", "E[X] + 3sigma[X]"))


#####
# Save the data
write.table(y, file = paste("data/", filename, ".csv", sep=""), row.names=FALSE, col.names=FALSE)