R/sim.r

Defines functions sim

Documented in sim

##' Simulate a random walk with time-varying autocorrelation
##'
##' @title Simulate a random walk with time-varying autocorrelation
##' @param Nt number of time steps to simulate
##' @param sigma_g sd of the AR(1) time-series for g
##' @param Sigma covar matrix for RW process
##' @param tstep increment, in days, in time-series
##' @return a dataframe with components:
##' \item{\code{dates}}{timestamps for each location}
##' \item{\code{lon}}{longitude without obs error}
##' \item{\code{lat}}{latitude without obs error}
##' \item{\code{gamma}}{autocorrelations}
##' @importFrom mvtnorm rmvnorm rmvt
##' @importFrom tibble as_data_frame data_frame
##' @export

sim <-
  function(Nt = 1000,
           sigma = c(0.2, 0.2),
           rho = 0,
           tstep = 0.5) {
    start.date <- ISOdate(2016, 5, 31, 00, 00, 00)

    x <- matrix(NA, Nt, 2)
    d <- matrix(NA, Nt - 1, 2)

    # simulate gamma as a slow cycle between low and high autocorrelation
    #   as per Auger-Methe et al. 2017
    g <- rep((cos(seq(-pi, pi, l = Nt / 3)) + 1.01) / 2.02, l = Nt)
    # lg <- -2
    # for(i in 2:Nt) {
    #   lg[i] <- rnorm(1, lg[i-1], 0.2)
    # }
    # g <- plogis(lg)

    Sigma <-
      matrix(c(sigma[1] ^ 2,
               sigma[1] * sigma[2] * rho,
               sigma[1] * sigma[2] * rho,
               sigma[2] ^ 2),
             2,
             2)

    ## Start as simple RW
    proc.err <- rmvnorm(Nt, c(0, 0), Sigma)
    x[1,] <- rmvnorm(1, c(0, 0), Sigma)
    x[2,] <- x[1, ] + proc.err[1, ]
    d[1,] <- x[2,] - x[1,]

    ## RW on displacements
    for (t in 2:(Nt - 1)) {
      d[t,] <- d[t - 1,] * g[t] + proc.err[t, ]
    }
    ## calculate locations from displacements
    for (t in 3:Nt) {
      x[t,] <- x[t - 1,] + d[t - 1,]
    }

    ## time interval is nominally 1 h
    dt <- seq(start.date, start.date + (Nt - 1) * tstep * 86400,
              by = tstep * 86400)

    data_frame(
      date = dt,
      lon = x[, 1],
      lat = x[, 2],
      g = g
    )

  }
ianjonsen/bssm documentation built on July 3, 2017, 10:33 p.m.