STGARCH: STGARCH

#' Simulation of 1D STARCH processes
#'
#' @param param Numeric vector of length 2, consiting of \eqn{\alpha_0} and \eqn{\alpha_1}
#' @param dim Numeric vector of length 2, indicating dimension size on form spatial x temporal
#' @param burnin Integer, size of temporal burnin
#' @param cutoff Integer, number of spatial points to be "cut off" in each spatial direction in order
#' to make the simulated process non-circular. If \code{cutoff=0}, the process is circular.
#' @return Simulated STARCH process in form of a matrix of dimension \code{dim}.
#' @examples
#' data<-simulation_1d_arch(param=c(0.3, 0.1),
#'                          dim=c(10, 100),
#'                          burnin=100,
#'                          cutoff=0)
#' @export
simulation_1d_arch<-function(param=NULL, dim=NULL, burnin=100, cutoff=0){

  #Check for correct input:
  if(length(dim)!=2) stop("Incorrect dimension specification")
  if(length(param)!=2) stop("Incorrect parameter specification")

  #Creating constants from vectors to ease notation
  a0<-param[1]
  a1<-param[2]
  X<-dim[1]+2*cutoff
  n<-dim[2]+burnin

  #Pre-allocating:
  data<-array(NA_real_,dim=c(X,n))
  data[,1]=rnorm(X,0,1)
  S<-array(NA_real_, dim=c(X,n))
  S[,1]<-abs(rnorm(X,0,1))
  # Circular simulation routine
  for(t in 2:n){
    a.tmp<-data[,t-1]^2
    a.tmp<-a.tmp[c(X,1:(X-1))]+a.tmp[1:X]+a.tmp[c(2:X,1)]
    S[,t]<-sqrt(a0+a1*a.tmp)
    data[,t]<-S[,t]*rnorm(X,0,1)
  }
  #Returning the simulated dataset with cutoff and burn-in removed
  return(data[(1+cutoff):(X-cutoff),(burnin+1):n])
}

#' Simulation of 2D STARCH-process
#'
#' @param dim Numeric vector of length 3, indicating dimension size on form spatial1 x spatial2 x temporal
#' @inheritParams simulation_1d_arch
#' @param seed For reproducibility, you may set a seed
#' @return Simulated STARCH process in form of a matrix of dimension \code{dim}.
#' @examples
#' data<-simulation_2d_arch(param=c(0.3, 0.05),
#'                          dim=c(10, 10, 100),
#'                          burnin=100,
#'                          cutoff=0)
#' @export


simulation_2d_arch<-function(param=NULL, dim=NULL, burnin=100, cutoff=0, seed=NULL){
  #Check for correct input:
  if(length(dim)!=3) stop("Incorrect dimension specification")
  if(length(param)!=2) stop("Incorrect parameter specification")

  if(!is.null(seed)) set.seed(seed)
  #Creating constants from vectors to ease notation
  a0<-param[1]; a1<-param[2];
  X<-dim[1]+2*cutoff; Y<-dim[2]+2*cutoff;n<-dim[3]+burnin;

  #Pre-allocating:
  data<-array(NA_real_,dim=c(X,Y,n))
  data[,,1]=rnorm(X*Y,0,1)
  S<-array(NA_real_, dim=c(X,Y,n))
  S[,,1]<-rnorm(X*Y,0,1)^2

  # Circular simulation routine
  for(t in 2:n){
    a.tmp<-data[,,t-1]^2
    a.tmp<-a.tmp[c(X,1:(X-1)),]+a.tmp[1:X,]+a.tmp[c(2:X,1),]
    a.tmp<-a.tmp[,c(Y,1:(Y-1))]+a.tmp[,1:Y]+a.tmp[,c(2:Y,1)]
    S[,,t]<-sqrt(a0+a1*a.tmp)
    data[,,t]<-S[,,t]*matrix(rnorm(n=X*Y),nrow=X, ncol=Y)

  }
  #Returning the simulated dataset with cutoff and burn-in removed
  return(data[(1+cutoff):(X-cutoff),(1+cutoff):(Y-cutoff),(burnin+1):n])
}

#' Simulation of 1D STGARCH-process
#'
#'
#' @inheritParams simulation_1d_arch
#' @param param Numeric vector of length 3, consiting of \eqn{\alpha_0}, \eqn{\alpha_1} and \eqn{\beta_1}
#' #' @return Simulated STARCH process in form of a matrix of dimension \code{dim}.
#' @examples
#' data<-simulation_1d_garch(param=c(0.3, 0.1, 0.05),
#'                           dim=c(10, 100),
#'                           burnin=100,
#'                           cutoff=0)
#' @export

simulation_1d_garch<-function(param=NULL, dim=NULL, burnin=100, cutoff=0){
  #Check for fitting values:
  if(length(dim)!=2) stop("Incorrect dimension specification")
  if(length(param)!=3) stop("Incorrect parameter specification")

  #creating constants from vectors to ease notation
  a0<-param[1];a1<-param[2]; b<-param[3];
  X<-dim[1]+cutoff;n<-dim[2]+burnin;

  #Pre-allocating:
  data<-array(NA_real_,dim=c(X,n))
  data[,1]=rnorm(X,0,1)
  S<-array(NA_real_, dim=c(X,n))
  S[,1]<-abs(rnorm(X,0,1))

  # Circular simulation routine
  for(t in 2:n){
    b.tmp<-S[,t-1]^2
    b.tmp<-b.tmp[c(X,1:(X-1))]+b.tmp[1:X]+b.tmp[c(2:X,1)]
    a.tmp<-data[,t-1]^2
    a.tmp<-a.tmp[c(X,1:(X-1))]+a.tmp[1:X]+a.tmp[c(2:X,1)]
    S[,t]<-sqrt(a0+a1*a.tmp+b*b.tmp)
    data[,t]<-S[,t]*rnorm(X,0,1)
  }

  return(data[(1+cutoff):(X-cutoff),(burnin+1):n])
}

#' Simulation of 2D STGARCH-process
#'
#'
#' @inheritParams simulation_1d_arch
#' @param param Numeric vector of length 3, consiting of \eqn{\alpha_0}, \eqn{\alpha_1} and \eqn{\beta_1}
#' @param dim Numeric vector of length 3, indicating dimension size on form spatial1 x spatial2 x temporal
#' @return Simulated STARCH process in form of a matrix of dimension \code{dim}.
#' @examples
#' data<-simulation_2d_garch(param=c(0.3, 0.05, 0.03),
#'                           dim=c(10, 10, 100),
#'                           burnin=100,
#'                           cutoff=0)
#' @export

simulation_2d_garch<-function(param=NULL, dim=NULL, burnin=100, cutoff=0){
  #Check for fitting values:
  if(length(dim)!=3) stop("Incorrect dimension specification")
  if(length(param)!=3) stop("Incorrect parameter specification")

  #creating constants from vectors to ease notation
  a0<-param[1];a1<-param[2]; b<-param[3];
  X<-dim[1]+2*cutoff;Y<-dim[2]+2*cutoff;n<-dim[3]+burnin;

  #Pre-allocating:
  data<-array(NA_real_,dim=c(X,Y,n))
  data[,,1]=rnorm(X*Y,0,1)
  S<-array(NA_real_, dim=c(X,Y,n))
  S[,,1]<-abs(rnorm(X*Y,0,1))

  # Circular simulation routine
  for(t in 2:n){
    b.tmp<-S[,,t-1]^2
    b.tmp<-b.tmp[c(X,1:(X-1)),]+b.tmp[1:X,]+b.tmp[c(2:X,1),]
    b.tmp<-b.tmp[,c(Y,1:(Y-1))]+b.tmp[,1:Y]+b.tmp[,c(2:Y,1)]
    a.tmp<-data[,,t-1]^2
    a.tmp<-a.tmp[c(X,1:(X-1)),]+a.tmp[1:X,]+a.tmp[c(2:X,1),]
    a.tmp<-a.tmp[,c(Y,1:(Y-1))]+a.tmp[,1:Y]+a.tmp[,c(2:Y,1)]
    S[,,t]<-sqrt(a0+a1*a.tmp+b*b.tmp)
    data[,,t]<-S[,,t]*matrix(rnorm(X*Y,0,1),nrow=X, ncol=Y)

  }

  return(data[(1+cutoff):(X-cutoff),(1+cutoff):(Y-cutoff),(burnin+1):n])
}