# R/gbm.R In spatPomp: Inference for Spatiotemporal Partially Observed Markov Processes

#### Documented in gbm

```#' Geometric Brownian motion spatPomp simulator
#'
#' Generate a spatPomp object representing a \code{U}-dimensional
#' geometric Brownian motion with spatial correlation decaying geometrically with
#' distance around a circle. The model is defined in continuous time, but
#' an Euler approximation is used for this numerical implementation.
#'
#' @param U A length-one numeric signifying dimension of the process.
#' @param N A length-one numeric signifying the number of time steps to evolve the process.
#' @param delta_t process simulations are performed every \code{delta_t} time units
#' @param delta_obs observations occur every \code{delta_obs} time units
#' @param IVP_values initial value parameters for the latent states
#' @return An object of class \sQuote{spatPomp} representing a simulation from a \code{U}-dimensional
#' geometric Brownian motion
#' @examples
#' g <- gbm(U=4, N=20)
#' # See all the model specifications of the object
#' spy(g)
#' @export
gbm <- function(U=5, N=100, delta_t=0.1, IVP_values = 1, delta_obs = 1){
dist <- function(u,v,n=U) min(abs(u-v),abs(u-v+U),abs(u-v-U))
dmat <- matrix(0,U,U)
for(u in 1:U) {
for(v in 1:U) {
dmat[u,v] <- dist(u,v)
}
}
to_C_array <- function(v)paste0("{",paste0(v,collapse=","),"}")
dist_C_rows <- apply(dmat,1,to_C_array)
dist_C_array <- to_C_array(dist_C_rows)
dist_C <- paste0("const double dist[",U,"][",U,"] = ",dist_C_array,"; ")
gbm_globals <- Csnippet(paste0("#define U ", U, " \n ", dist_C))

obs_names <- paste0("Y",(1:U))
gbm_data <- data.frame(time=rep((1:N)*delta_obs,U),unit=rep(obs_names,each=N),Y=rep(NA,U*N),stringsAsFactors=F)

gbm_unit_statenames <- c("X")
gbm_statenames <- paste0(gbm_unit_statenames,1:U)

gbm_IVPnames <- paste0(gbm_statenames,"_0")
gbm_RPnames <- c("rho","sigma","tau")
gbm_paramnames <- c(gbm_RPnames,gbm_IVPnames)

gbm_rprocess <- Csnippet("
double *X = &X1;
double Xbm[U];
double dW[U];
int u,v;
for (u = 0 ; u < U ; u++) {
Xbm[u] = log(X[u]);
dW[u] = rnorm(0,sigma*sqrt(dt));
}
for (u = 0 ; u < U ; u++) {
for (v=0; v < U ; v++) {
Xbm[u] += dW[v]*pow(rho,dist[u][v]);
}
X[u] = exp(Xbm[u]);
}
")

gbm_skel <- Csnippet("
double *DX = &DX1;
double *X = &X1;
double cumsigsq[U];
int u,v;
for (u = 0; u < U; u++){
cumsigsq[u] = 0;
}
for (u = 0 ; u < U ; u++) {
for (v = 0 ; v < U ; v++) {
cumsigsq[u] += pow(sigma*(pow(rho, dist[u][v])), 2);
}
DX[u] = X[u]*((cumsigsq[u])/2);
}
")

gbm_rinit <- Csnippet("
double *X = &X1;
const double *X_0=&X1_0;
int u;
for (u = 0; u < U; u++) {
X[u]=X_0[u];
}
")

gbm_dmeasure <- Csnippet("
const double *X = &X1;
const double *Y = &Y1;
double tol = pow(1.0e-18,U);
int u;
lik=0;
// Jacobian to get Y|logX since we know logY|logX is Normal(0,tau-squred)
for (u=0; u<U; u++)lik += (dnorm(log(Y[u]),log(X[u]),tau,1) + log(1/Y[u]));
if(!give_log) lik = exp(lik) + tol;
")

gbm_rmeasure <- Csnippet("
const double *X = &X1;
double *Y = &Y1;
double tol = pow(1.0e-18,U);
int u;
// Y|X = X*exp(eta) where eta is Normal(0, tau-squared)
for (u=0; u<U; u++) Y[u] = X[u]*exp(rnorm(0,tau+tol));
")

gbm_dunit_measure <- Csnippet("
lik = dnorm(log(Y),log(X),tau,1) + log(1/Y);
if(!give_log) lik = exp(lik);
")

gbm_runit_measure <- Csnippet("
double tol = pow(1.0e-18,U);
double Y;
Y = X*exp(rnorm(0,tau+tol));
")

gbm_eunit_measure <- Csnippet("
ey = X*exp(tau*tau/2);
")

gbm_munit_measure <- Csnippet("
M_tau = sqrt(log(0.5 + 0.5*sqrt(1 + (4*vc/(X*X)))));
")

gbm_vunit_measure <- Csnippet("
vc = X*X*(exp(2*tau*tau) - exp(tau*tau));
")

gbm_spatPomp <- spatPomp(gbm_data,
times="time",
t0=0,
units="unit",
unit_statenames = gbm_unit_statenames,
rprocess=euler(gbm_rprocess,delta.t = delta_t),
skeleton=vectorfield(gbm_skel),
paramnames=gbm_paramnames,
globals=gbm_globals,
rmeasure=gbm_rmeasure,
dmeasure=gbm_dmeasure,
dunit_measure=gbm_dunit_measure,
runit_measure=gbm_runit_measure,
eunit_measure=gbm_eunit_measure,
munit_measure=gbm_munit_measure,
vunit_measure=gbm_vunit_measure,
partrans = parameter_trans(logit = c("rho"), log = c("sigma", "tau")),
rinit=gbm_rinit
)
## We need a parameter vector. For now, we initialize the Brownian process at zero.
test_ivps <- rep(IVP_values,U)
names(test_ivps) <- gbm_IVPnames
test_params <- c(rho=0.4, sigma=1, tau=1, test_ivps)
simulate(gbm_spatPomp,params=test_params)

}
```

## Try the spatPomp package in your browser

Any scripts or data that you put into this service are public.

spatPomp documentation built on Jan. 15, 2022, 5:06 p.m.