R/sim.resp.R

In KironmoyDas/KD-STAT0035-GMupdate: Generalised Joint Regression Modelling

sim.resp <- function(margin, rsim, eta, sigma2, nu, setseed = TRUE){


if(margin == "ZTP") rZTP <- function(n, mu) qpois(runif(n, dpois(0, mu), 1), mu)

if(margin %in% c("DGP","DGPII")){

       rDGP <- function(n, mu, sigma, margin){ # mu is eta here

                minp <- distrHsATDiscr(0, mu, sigma, NULL, margin, NULL, robust = TRUE, min.dn = 1e-323, min.pr = 1e-16, max.pr = 0.999999)$pdf2
                
                p    <- runif(rsim)

                indv <- p <= minp
                

                if(margin == "DGP")   q <- ifelse(indv == TRUE, 0, ceiling( sigma/mu*(     (1 - p)^(-mu)   - 1   )) - 1 )
                if(margin == "DGPII") q <- ifelse(indv == TRUE, 0, ceiling( sigma/(mu^2)*( (1 - p)^(-mu^2) - 1   )) - 1 )
                
                q
                
                                      }
}


if(setseed == TRUE) set.seed(1)

if(margin == "N")     y <- rNO(   rsim,    mu =     eta,     sigma = sigma2) 
if(margin == "GU")    y <- rGU(   rsim,    mu =     eta,     sigma = sigma2) 
if(margin == "rGU")   y <- rRG(   rsim,    mu =     eta,     sigma = sigma2) 
if(margin == "LO")    y <- rLO(   rsim,    mu =     eta,     sigma = sigma2) 
if(margin == "LN")    y <- rLOGNO(rsim,    mu =     eta,     sigma = sigma2) 
if(margin == "WEI")   y <- rWEI(  rsim,    mu = exp(eta),    sigma = sigma2) 
if(margin == "iG")    y <- rIG(   rsim,    mu = exp(eta),    sigma = sigma2) 
if(margin == "GA")    y <- rGA(   rsim,    mu = exp(eta),    sigma = sigma2) 
if(margin == "GAi")   y <- rGA(   rsim,    mu =     eta,     sigma = sigma2)






if(margin == "TW"){    


     if(rsim == 1) y <- rTweedie(mu = exp(eta), phi = sigma2, p = nu)
     
     if(rsim > 1){
           y <- NA
           if(length(sigma2) == 1) sigma2 <- rep(sigma2, rsim) 
           if(length(nu)     == 1) nu     <- rep(nu, rsim) 
           for(i in 1:rsim) y[i] <- rTweedie(mu = exp(eta[i]), phi = sigma2[i], p = nu[i])                
                 }     
     
}


if(margin %in% c("GP","GPII","GPo")){

             if(rsim == 1){ 
             
               if(margin == "GP")    y <- rgpd(  rsim, loc = 0, shape = eta,            scale = sigma2)
               if(margin == "GPII")  y <- rgpd(  rsim, loc = 0, shape = exp(eta) - 0.5, scale = sigma2)
               if(margin == "GPo")   y <- rgpd(  rsim, loc = 0, shape = exp(eta) - 0.5, scale = sigma2/(1+(exp(eta) - 0.5)))
               
             
             
             }
             
             if(rsim > 1){
                y <- NA
                if(length(sigma2) == 1) sigma2 <- rep(sigma2, rsim) 
                for(i in 1:rsim){
                  if(margin == "GP")     y[i] <- rgpd(  1, loc = 0, shape = eta[i],            scale = sigma2[i])
                  if(margin == "GPII")   y[i] <- rgpd(  1, loc = 0, shape = exp(eta[i]) - 0.5, scale = sigma2[i])
                  if(margin == "GPo")    y[i] <- rgpd(  1, loc = 0, shape = exp(eta[i]) - 0.5, scale = sigma2[i]/( 1 + (exp(eta[i]) - 0.5) )  )
                }
             
             
             }
}


if(margin == "DAGUM") y <- rGB2(  rsim,    mu = exp(eta),    sigma = sigma2, nu = nu, tau = 1) 
if(margin == "SM")    y <- rGB2(  rsim,    mu = exp(eta),    sigma = sigma2, nu = 1 , tau = nu) 
if(margin == "BE")    y <- rBE(   rsim,    mu = plogis(eta), sigma = sigma2) 
if(margin == "FISK")  y <- rGB2(  rsim,    mu = exp(eta),    sigma = sigma2, nu = 1 , tau = 1)
if(margin == "NBI")   y <- rNBI(  rsim,    mu = exp(eta),    sigma = sigma2) 
if(margin == "NBII")  y <- rNBII( rsim,    mu = exp(eta),    sigma = sigma2) 
if(margin == "PIG")   y <- rPIG(  rsim,    mu = exp(eta),    sigma = sigma2) 
if(margin == "PO")    y <- rPO(   rsim,    mu = exp(eta)) 
if(margin == "ZTP")   y <- rZTP(  rsim,    mu = exp(eta)) 

if(margin %in% c("DGP", "DGPII"))   y <- rDGP(  rsim,    mu = eta,         sigma = sigma2, margin = margin) # in dgpII eta^2 is done internally

if(margin %in% c("probit", "logit", "cloglog"))  y <- rbinom(rsim, 1, prob = probm(eta, margin, min.dn = 1e-40, min.pr = 1e-16, max.pr = 0.999999)$pr )
                                                

if(setseed == TRUE) rm(list = ".Random.seed", envir = globalenv()) 

y

}