R/warSim.R

In WRI: Wasserstein Regression and Inference

#' Generate simulation data
#'
#' @description Generate WAR(p) simulation data sets: samples simulated from a WAR(3) model similar to the specification in Section 4 of the referenced paper.
#'
#' @return A list of:
#'\item{sample.ts}{one simulation run chosen from \code{sample.ts.full}}
#'\item{sample.ts.full}{1000 simulation runs, each of which consists of a sample time series (of length 100) of quantile functions generated by a WAR(3) model as specified by the reference paper}
#'\item{quantile.grid}{the grid over which the quantile functions in sample.ts.full are evaluated}
#'
#' @references
#' \cite{Wasserstein Autoregressive Models for Density Time Series, Chao Zhang, Piotr Kokoszka, Alexander Petersen, 2022}
#'
#' @export

warSim <- function() {

  # AR parameters
  phi.1 <- 0.825
  phi.2 <- -0.1875
  phi.3 <- 0.0125

  # the Wasserstein mean is unifrom (0,1)
  grid <- seq(0,1, length.out = 101)

  # burn in noise
  set.seed(131)
  burn.in.eta <- matrix(rnorm(100*100, 0, 1), nrow = 100, ncol = 100, byrow = TRUE)
  burn.in.delta <- matrix(runif(100*100, min = -0.2, max = 0.2), nrow = 100, ncol = 100, byrow = TRUE)
  burn.in.noise <- array(NA, dim = c(101, 100, 100),
                         dimnames = list(paste("grid", 1:101, sep = " "), paste("time", 1:100, sep = " "), paste("replicate", 1:100, sep = " ")))
  for(k in 1:100){
    temp.burn.in.eta <- burn.in.eta[k,]
    temp.burn.in.delta <- burn.in.delta[k,]
    for(j in 1:100){
      for(i in 1:101){
        burn.in.noise[i,j,k] <- temp.burn.in.eta[j] + sin(temp.burn.in.delta[j]*grid[i])
      }
    }
  }

  ################################
  ### Simulation case: n = 100 ###
  ################################

  # innovations
  set.seed(131)
  eta <- matrix(rnorm(100*100, 0, 1), nrow = 100, ncol = 100, byrow = TRUE)
  delta <- matrix(runif(100*100, min = -0.2, max = 0.2), nrow = 1000, ncol = 100, byrow = TRUE)
  inno <- array(NA, dim = c(101, 100, 100),
                dimnames = list(paste("grid", 1:101, sep = " "), paste("time", 1:100, sep = " "), paste("replicate", 1:100, sep = " ")))
  for(k in 1:100){
    temp.eta <- eta[k,]
    temp.delta <- delta[k,]
    for(j in 1:100){
      for(i in 1:101){
        inno[i,j,k] <- temp.eta[j] + sin(temp.delta[j]*grid[i])
      }
    }
  }

  # create an empty array to store simulated time series
  sample.array <- array(NA, dim = c(101,100,100))

  for(k in 1:100){
    for(i in 1:101){
      sample.array[i, ,k] <- arima.sim(model=list(order=c(3,0,0), ar=c(phi.1,phi.2,phi.3)), n=100, innov=inno[i, ,k],
                                       n.start=100, start.innov=burn.in.noise[i, ,k])
    }
  }

  warSimData <- list("sample.ts"=sample.array[, , 2],
                     "sample.ts.full"=sample.array,
                     "quantile.grid"=grid)
  return(warSimData)
}