R/WACSsimul.R

In WACS: Multivariate Weather-State Approach Conditionally Skew-Normal Generator

###################################################################
#
# This function is part of WACSgen V1.0
# Copyright © 2013,2014,2015, D. Allard, BioSP,
# and Ronan Trépos MIA-T, INRA
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details. http://www.gnu.org
#
###################################################################

#' Performs simulations based on estimated parameters of the WACS model
#' 
#' @import mnormt
#' @import tmvtnorm
#' 
#' @export
#' 
#' @param wacspar Parameters of the WACS model estimated with WACSestim
#' @param from    Starting date of the simulation (format: "yyyy-mm-dd")
#' @param to      Ending date of the simulation (format: "yyyy-mm-dd")
#' @param first.day Conditioning values for first day (optional)              
#' @param REJECT  Boolean indicating whether a rejection technique is used to guarantee variables within bounds.   
#'                Default is FALSE. In this case, values outside bounds are forced to the bounds.
#'
#' @return A list containing the simulation results
#'
#' @examples
#' \dontrun{
#'   ## Simple example
#'   data(ClimateSeries)
#'   ThisData = WACSdata(ClimateSeries)
#'   ThisPar  = WACSestim(ThisData)
#'   ThisSim  = WACSsimul(ThisPar, from="1995-01-01", to="2012-12-31")
#' }
#' @note
#' 
#' Variables are simulated sequentially: day d is simulated conditionally on the values at day (d-1).
#' If \code{REJECT=TRUE}, a rejection technique is used to force simulated variables within the bounds. 
#' If \code{REJECT=FALSE}, variables that could have been simulated outside the bounds are forced to the limits.
#' The rejection technique tends to produce biases. \code{REJECT=FALSE} is thus recommended



WACSsimul=function(wacspar, from, to, first.day= NULL, REJECT=FALSE){    
  
  # Checking wacspar
  if (class(wacspar) != "WACSpar") {
    stop ("[WACSsimul] Parameters should be of class 'WACSpar', as generated by calling WACSestim")
  }
  
  ###############################################
  #    
  # If too many years, we run several shorter simulations, and concatenate the result at the end 
  #
  ###############################################
  
  Nv     = length(wacspar$varnames)
  Nyears = ( wacs.year(to) - wacs.year(from) ) + 1
  if (Nyears < 40){
    Xsimul = wacs.simul_innercall(wacspar, from, to, first.day, REJECT)
  }else{
    start.day   = wacs.day(from)
    start.month = wacs.month(from)
    start.year  = wacs.year(from)
    Nstep       = Nyears %/% 20 
    
    # First step
    innerto.year  = start.year + 20
    innerto       = as.Date(paste(innerto.year,"-",start.month,"-",start.day,sep=""))
    Xout          = wacs.simul_innercall(wacspar, from, innerto, first.day, REJECT)
    last          = dim(Xout)[1]
    last.day      = as.numeric(Xout[last,][6:(5+Nv)])
    Xsimul        = Xout[-last,]
    
    
    # Steps within loop
    for (istep in 2: Nstep){
      innerfrom.year = start.year + (istep-1)*20
      innerto.year   = start.year + istep*20
      innerfrom = as.Date(paste(innerfrom.year,"-",start.month,"-",start.day,sep=""))
      innerto   = as.Date(paste(innerto.year,"-",start.month,"-",start.day,sep=""))
      first.day  = last.day
      Xout       = wacs.simul_innercall(wacspar, innerfrom, innerto, first.day, REJECT)
      last       = dim(Xout)[1]
      last.day   = as.numeric(Xout[last,][6:(5+Nv)])
      Xsimul     = rbind(Xsimul,Xout[-last,])
    }
    
    # Last step
    
    if (as.numeric(difftime(to,innerto))>0){
      first.day = last.day
      Xout      = wacs.simul_innercall(wacspar, innerto, to, first.day, REJECT)
      Xsimul    = rbind(Xsimul,Xout)
    }
  }
  res = list(sim=Xsimul)
  class(res) = "WACSsimul"
  return(res)
}