inst/doc/vignette.R

In CoSMoS: Complete Stochastic Modelling Solution

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(eval = TRUE,
                      echo = TRUE,
                      fig.width = 7, 
                      warning = FALSE,
                      message = FALSE)
library(CoSMoS)
library(plot3D)

## ---- warning=FALSE, message=FALSE--------------------------------------------
## (i) specifying the sample size
no <- 1000
## (ii) defining the type of marginal distribution and its parameters
marginaldist <- "ggamma"
param <- list(scale = 1, shape1 = .8, shape2 = .8)
## (iii) defining the desired autocorrelation
acf.my <- c(1, 0.8)
## (iv) simulating
ggamma_sim <- generateTS(n = no, margdist = marginaldist, margarg = param, acsvalue = acf.my)
## and (v) visually checking the generated time series
quickTSPlot(ggamma_sim[[1]]) 


## ---- warning=FALSE, message=FALSE--------------------------------------------
acf <- c(1, 0.6, 0.5, 0.4, 0.3) #up to lag-4
ggamma_sim <- generateTS(n = no, margdist = marginaldist, margarg = param, acsvalue = acf)
quickTSPlot(ggamma_sim[[1]])

## ---- warning=FALSE, message=FALSE--------------------------------------------
## specify lag
lags <- 0:10

## get the ACS
f <- acs(id = "fgn", t = lags, H = .75)
b <- acs(id = "burrXII", t = lags, scale = 1, shape1 = .6, shape2 = .4)
w <- acs(id = "weibull", t = lags, scale = 2, shape = 0.8)
p <- acs(id = "paretoII", t = lags, scale = 3, shape = 0.3)

## visualize the ACS
dta <- data.table(lags, f, b, w, p)
m.dta <- melt(data = dta, id.vars = "lags")

ggplot(data = m.dta, mapping = aes(x = lags, y = value, group = variable, colour = variable)) + 
  geom_point(size = 2.5) + geom_line(lwd = 1) +
  scale_color_manual(values = c("steelblue4", "red4", "green4", "darkorange"), 
  labels = c("FGN", "Burr XII", "Weibull", "Pareto II"), name = "") +
  labs(x = bquote(lag ~ tau), y = "ACS") + scale_x_continuous(breaks = lags) + theme_light()

## ---- warning=FALSE, message=FALSE--------------------------------------------
acf <- acs(id = "paretoII", t = 0:30, scale = 1, shape = .75)
ggamma_sim <- generateTS(n = no, margdist = marginaldist, margarg = param, acsvalue = acf)
dta <- data.frame(time = 1:no, value = ggamma_sim[[1]])

quickTSPlot(dta$value)

## ---- warning=FALSE, message=FALSE--------------------------------------------
my_acf <- exp(seq(0, -2, -0.1))
ggamma_sim <- generateTS(n = no, margdist = marginaldist, margarg = param, acsvalue = my_acf)
quickTSPlot(ggamma_sim[[1]])


## ---- fig.height = 7, warning=FALSE, message=FALSE----------------------------
prob_zero <- .9
## the argument `TSn = 5` enables the simulation of 5 timeseries.
ggamma_sim <- generateTS(n = no, margdist = marginaldist, margarg = param, acsvalue = acf, 
                         p0 = prob_zero, TSn = 5)
plot(x = ggamma_sim, main = "") + theme_light()

## ---- warning=FALSE, message=FALSE--------------------------------------------
checkTS(ggamma_sim)

## ---- fig.height = 7, warning=FALSE, message=FALSE----------------------------

## specify grid of spatial and temporal lags
d <- 51
st <- expand.grid(0:(d - 1), 0:(d - 1))

## get the STCS
wc <- stcfclayton(t = st[, 1], s = st[, 2], scfid = "weibull", tcfid = "weibull", copulaarg = 2,
                  scfarg = list(scale = 20, shape = 0.7), tcfarg = list(scale = 5.1,shape = 0.8))

## visualize the STCS
wc.m <- matrix(data = wc, nrow = d)
j <- tail(which(wc.m[1, ] > 0.15), 1)
i <- tail(which(wc.m[, 1] > 0.15), 1)
wc.m <- wc.m[1:i, 1:j]

persp3D(z = wc.m, x = 1: nrow(wc.m), y = 1:ncol(wc.m),
expand = 1, main = "", scale = TRUE, facets = TRUE,
xlab="Time lag", ylab = "Distance", zlab = "STCF", colkey = list(side = 4, length = 0.5),
phi = 20, theta = 120, resfac = 5,  col= gg2.col(100))

## ---- fig.height = 5, warning=FALSE, message=FALSE----------------------------
## set a sequence of hypothetical coordinates 
d <- 5
coord <- cbind(runif(d)*30, runif(d)*30)

## compute VAR model parameters  
fit <- fitVAR(spacepoints = coord, 
              p = 4, 
              margdist = "burrXII",
              margarg = list(scale = 3, shape1 = .9, shape2 = .2), 
              p0 = 0.8, 
              stcsid = "clayton",
              stcsarg = list(scfid = "weibull", tcfid = "weibull", copulaarg = 2,
                   scfarg = list(scale = 25, shape = 0.7), 
                   tcfarg = list(scale = 3.1, shape = 0.8) ) )

## generate correlated timeseries  
sim <- generateMTS(n = 500, STmodel = fit)

## visualize simulated timeseries
dta <- melt(data = data.table(time = 1:nrow(sim), sim[,1:d]), id.vars = "time")

ggplot(data = dta, mapping = aes(x = time, y = value)) + geom_line() +
       facet_grid(facets = variable ~ ., scales = "free_y") + theme_light()

## ---- fig.height = 5, warning=FALSE, message=FALSE----------------------------
## set a sequence of hypothetical coordinates
d <- 5
coord <- cbind(runif(d)*30, runif(d)*30)

## fit and generate correlated timeseries
sim <- generateMTSFast(n = 500, 
                       spacepoints = coord, 
                       p0 = 0.7, 
                       margdist ="burrXII",
                       margarg = list(scale = 3, shape1 = .9, shape2 = .2),  
                       stcsarg = list(scfid = "weibull", tcfid = "weibull",
                       scfarg = list(scale = 25, shape = 0.7),
                       tcfarg = list(scale = 3.1, shape = 0.8)) )

## visualize simulated timeseries
dta <- melt(data = data.table(time = 1:nrow(sim), sim[,1:d]), id.vars = "time")

ggplot(data = dta, mapping = aes(x = time, y = value)) + geom_line() +
       facet_grid(facets = variable ~ ., scales = "free_y") + theme_light()

## ---- warning=FALSE, message=FALSE--------------------------------------------
## compute VAR model parameters
## CPU time: ~15s 
fit <- fitVAR(spacepoints = 20, p = 3, margdist ="burrXII",
              margarg = list(scale = 3, shape1 = .9, shape2 = .2), p0 = 0.8, stcsid = "clayton",
              stcsarg = list(scfid = "weibull", tcfid = "weibull", copulaarg = 2,
              scfarg = list(scale = 20, shape = 0.7), tcfarg = list(scale = 1.1, shape = 0.8) ) )

## generate isotropic random fields with nonseparable correlation structure 
sim1 <- generateRF(n = 1000, STmodel = fit)

## fast simulation of isotropic random fields with separable correlation structure
sim2 <- generateRFFast(n = 1000, spacepoints = 20, p0 = 0.7, margdist ="paretoII",
                       margarg = list(scale = 1, shape = .3),
                       stcsarg = list(scfid = "weibull", tcfid = "weibull",
                       scfarg = list(scale = 20, shape = 0.7), 
                       tcfarg = list(scale = 1.1, shape = 0.8)))


## ---- fig.height = 6.5, warning=FALSE, message=FALSE--------------------------
## check random fields
## CPU time: ~20s 
checkRF(RF = sim1, nfields = 9*9, method = "stat")
checkRF(RF = sim1, nfields = 9*9, method = "statplot")
checkRF(RF = sim1, nfields = 9*9, method = "field")

## ---- fig.width = 4, fig.height = 4, warning=FALSE, message=FALSE-------------
## specify a grid of coordinates
m = 30
aux <- seq(0, m - 1, length = m)
coord <- expand.grid(aux, aux)

## transform the coordinate system
at <- anisotropyTswirl(spacepoints = coord, x0 = floor(m / 2), y0 = floor(m / 2), 
                       b = 10, alpha = 1.8 * pi)

## visualize transformed coordinate system
  aux = data.frame(lon = at[ ,1], lat = at[ ,2], id1 = rep(1:m, each = m), id2 = rep(1:m, m))
  ggplot(aux, aes(x = lon, y = lat)) + 
    geom_path(aes(group = id1)) + geom_path(aes(group = id2)) + geom_point(col = 2) + theme_light()


## ---- fig.height = 7, warning=FALSE, message=FALSE----------------------------
## compute VAR model parameters
fit <- fitVAR(spacepoints = m, p = 1, margdist = 'burrXII', 
              margarg = list(scale = 3, shape1 = .9, shape2 = .2), p0 = 0.1, stcsid = "clayton", 
              stcsarg = list(scfid = "weibull", tcfid = "weibull", copulaarg = 2,
                        scfarg = list(scale = 2, shape = 0.7), tcfarg = list(scale = 20.1, shape = 0.8)),
              anisotropyid = "swirl", 
              anisotropyarg = list(x0 = floor(m / 2), y0 = floor(m / 2), b = 10, alpha = 1.8 * pi) )

## generate
set.seed(9)
sim3 <- generateRF(n=25, STmodel = fit)

## check
checkRF(RF = sim3, nfields = 5*5, method = "field")


## ---- fig.height = 7.0, warning=FALSE, message=FALSE--------------------------
## compute VAR model parameters
fit <- fitVAR(spacepoints = 30, p = 1, margdist = 'burrXII', 
              margarg = list(scale = 3, shape1 = .9, shape2 = .2), p0 = 0.8, stcsid = "clayton", 
              stcsarg = list(scfid = "weibull", tcfid = "weibull", copulaarg = 2, 
                      scfarg = list(scale = 20, shape = 0.7), tcfarg = list(scale = 30.1, shape = 0.8)),
              anisotropyid = "affine", 
              anisotropyarg = list(phi1 = 2., phi2 = 0.5, phi12 = 0, theta = -pi/3),
              advectionid = "uniform", advectionarg = list(u = 2.5, v = 1.5) )

## generate
set.seed(10) # 1 # 5
sim3 <- generateRF(n=81, STmodel = fit)

## check
checkRF(RF = sim3, nfields = 9*9, method = "field")


## ---- warning=FALSE, message=FALSE--------------------------------------------
data("precip")
quickTSPlot(precip$value)

## ---- fig.height = 9, warning=FALSE, message=FALSE----------------------------
## CPU time: ~75s 
precip_ggamma <- analyzeTS(TS = precip, season = "month", dist = "ggamma",
                           acsID = "weibull", lag.max = 12)

reportTS(aTS = precip_ggamma, method = "dist") + theme_light()
reportTS(aTS = precip_ggamma, method = "acs") + theme_light()
reportTS(aTS = precip_ggamma, method = "stat")

## ---- fig.height = 9, warning=FALSE, message=FALSE----------------------------
precip_pareto <- analyzeTS(TS = precip, season = "month", dist = "paretoII", acsID = "fgn", lag.max = 12)

reportTS(aTS = precip_pareto, method = "dist")+ theme_light()
reportTS(aTS = precip_pareto, method = "acs") + theme_light()

## ---- warning=FALSE, message=FALSE--------------------------------------------
sim_precip <- simulateTS(aTS = precip_ggamma, from = as.POSIXct(x = "1978-12-01 00:00:00"),
                         to = as.POSIXct(x = "2008-12-01 00:00:00"))
dta <- precip
dta[, id := "observed"]
sim_precip[, id := "simulated"]
dta <- rbind(dta, sim_precip)

ggplot(data = dta) + geom_line(mapping = aes(x = date, y = value)) + 
  facet_wrap(facets = ~id, ncol = 1) + theme_light()

## ---- fig.height = 9, warning=FALSE, message=FALSE----------------------------
## CPU time: ~240s
data("disch")

str <- analyzeTS(TS = disch, dist = "lnorm", norm = "N2", acsID = "paretoII", 
                 lag.max = 20, constrain = TRUE, season = "month")

reportTS(aTS = str) + theme_light()
reportTS(aTS = str, method = "stat")

## ---- fig.height = 3.5, warning=FALSE, message=FALSE--------------------------
sim_str <- simulateTS(aTS = str)

dta <- disch
dta[, id := "observed"]
sim_str[, id := "simulated"]
dta <- rbind(dta, sim_str)

ggplot(data = dta) + geom_line(mapping = aes(x = date, y = value)) + 
  facet_wrap(facets = ~id, ncol = 1) + theme_light()