R/SimARTAp.R

In itsoukal/anySim: Stochastic Simulation of Processes with any Marginal Distribution and Correlation Structure

#' @title Simulation of stationary process via ARTA(p) model.
#'
#' @description Simulation of stationary process with target marginal distribution and autocorrelation structure via ARTA(p) model.
#'
#' @param ARTApar A list containing the parameters of the model. The list is generated by function "EstARTAp".
#' @param burn A scalar specifying the length of burn-out sample.
#' @param steps A scalar specifying the length of the time series to be generated.
#' @param stand A boolean (T or F) indicating whether to standardize (or not) the auxiliary Gaussian time series prior to their mapping to the actual domain. The default value is FALSE.
#'
#' @return A list of the 3 generated time series (in vector format):
#' X: The final time series at the target domain with the target marginal distribution and correlation structure;
#' Z: The auxiliary Gaussian time series at the Gaussian domain and;
#' U: The auxiliary uniform time series at the Copula domain (i.e., in [0,1]).
#'
#' @export
#'
#' @examples
#' ## Simulation of univariate stationary process with Gamma marginal distribution 
#' ## and autocorrelation structure given by the product of a CAS and a periodic ACS.
#'\dontrun{
#' set.seed(12)
#' 
#' # Define the target autocorrelation structure.
#' acsS=csCAS(param=c(3,0.6),lag=1000) # Stationary CAS with b=3 and k=0.6.
#' acsP=csPeriodic(param=c(12,1.5),lag=1000) # Periodic ACS with p=12 and l=1.5.
#' ACS=csP*acsS # The target ACS as product of the two previous ones.
#' 
#' # Define the target distribution function (ICDF).
#' FX='qgamma' # the Gamma distribution
#' 
#' # Define the parameters of the target distribution.
#' pFX=list(shape=5,scale=1)
#' 
#' # Estimate the parameters of the auxiliary Gaussian AR(p) model.
#' ARTApar=EstARTAp(ACF=ACS,maxlag=0,dist=FX,params=pFX,NatafIntMethod='GH')
#' 
#' # Generate a synthetic series of 10000 length. 
#' SynthARTAcont=SimARTAp(ARTApar=ARTApar,steps=10^5)
#'}
#'
#' ## Simulation of univariate stationary process with discrete marginal distribution
#' ## (Beta-Binomial) and autocorrelation structure given by CAS.
#'\dontrun{
#' set.seed(16)
#' 
#' # Define the target autocorrelation structure.
#' ACS=acsCAS(param=c(1.5,0.3),lag=1000) # CAS with b=1.5 and k=0.3.
#' 
#' # Define the target distribution function (ICDF).
#' require(TailRank)
#' FX='qbb' # the Beta-Binomial distribution.
#' 
#' # Define the parameters of the target distribution.
#' pFX=list(N=10,u=3,v=10)
#' 
#' # Estimate the parameters of the auxiliary Gaussian AR(p) model.
#' ARTApar=EstARTAp(ACF=ACS,maxlag=0,dist=FX,params=pFX,NatafIntMethod="MC")
#' 
#' # Generate a synthetic series of 10000 length. 
#' SynthARTAdiscr=SimARTAp(ARTApar=ARTApar,steps=10^5)
#'}
#'
#' ## Simulation of univariate stationary process with zero-inflated marginal distribution 
#' ## (Gen. Gamma for the continuous part) and autocorrelation structure given by CAS.
#'\dontrun{
#' set.seed(18)
#' 
#' # Define the target autocorrelation structure.
#' ACS=acsCAS(param=c(0.91,1.09),lag=1000) # CAS with b=0.91 and k=1.09.
#' 
#' # Define the target distribution function (ICDF).
#' FX='qmixed' # Define that distribution is of zero-inflated type.
#' 
#' # Define the distribution for the continuous part of the process.
#' # Here, a re-parameterized version of Gen. Gamma distribution is used.
#' qgengamma=function(p,scale,shape1,shape2){
#'   require(VGAM)
#'   X=qgengamma.stacy(p=p,scale=scale,k=(shape1/shape2),d=shape2)
#'   return(X)
#' } 
#' 
#' # Define the parameters of the zero-inflated distribution function.
#' pFX=list(Distr=qgengamma,p0=0.8,scale=0.25,shape1=1.16,shape2=0.54) 
#' 
#' # Estimate the parameters of the auxiliary Gaussian AR(p) model.
#' ARTApar=EstARTAp(ACF=ACS,dist=FX,params=pFX,NatafIntMethod="GH",NoEval=9,polydeg=0)
#' 
#' # Generate a synthetic series of 10000 length. 
#' SynthARTAzi=SimARTAp(ARTApar=ARTApar,steps=10^5)
#'}
SimARTAp=function(ARTApar, burn=1000, steps=10000, stand=F){
  phi=ARTApar$phi
  sigma=ARTApar$simga
  dist=ARTApar$dist
  params=ARTApar$params

  lag=length(phi)
  Start=rnorm(lag,mean = 0,sd = 1)
  Z=mat.or.vec((burn+steps),1)
  Z[1:lag]=Start

  totalSteps=burn+steps
  for (i in (lag+1):totalSteps){
    Z[i]=Z[i-c(1:lag)]%*%phi+rnorm(1,mean = 0,sd=sigma)
  }
  Z=Z[(burn+1):totalSteps]

  if (stand==T){
    U=pnorm(Z,mean = mean(Z), sd=sd(Z))
    Z=qnorm(U)
  } else {
    U=pnorm(Z)
  }

  funcall <- as.call(c(as.name(dist), list(U), params))
  X=eval(funcall)

  return(list('X'=X,'Z'=Z,'U'=U))
}