```{r jags-setup, include=FALSE, purl=FALSE}
removefiles <- c()
knitr::opts_knit$set(unnamed.chunk.label = "jags-")
```


# JAGS for Bayesian time series analysis  {#chap-jags}
\chaptermark{JAGS}

In this lab, we will illustrate how to use JAGS to fit time series models with Bayesian methods. The purpose of this chapter is to teach you some basic JAGS models. To go beyond these basics, study the wide variety of software tools to do time series analysis using Bayesian methods, e.g. packages listed on the R Cran [TimeSeries](http://cran.r-project.org/web/views/TimeSeries.html) task view.   

A script with all the R code in the chapter can be downloaded  [here](./Rcode/intro-to-jags.R). The Rmd for this chapter can be downloaded [here](./Rmds/intro-to-jags.Rmd).


### Data and packages {-}

For data for this lab, we will use a dataset on air quality in New York. For the majority of our models, we are going to treat wind speed as the response variable for our time series models. 

```{r jags-loaddata, echo=TRUE, results='hide', eval=TRUE}
data(airquality, package = "datasets")
Wind <- airquality$Wind # wind speed
Temp <- airquality$Temp # air temperature
N <- dim(airquality)[1] # number of data points
```

To run this code, you will need to install JAGS for your operating platform using the instructions [here](http://sourceforge.net/projects/mcmc-jags/files/).  Click on JAGS, then the most recent folder, then the platform of your machine.  You will also need the **coda**, **rjags** and **R2jags** packages.
```{r jags-loadpackages, results='hide', message=FALSE, warnings=FALSE}
library(coda)
library(rjags)
library(R2jags)
```

## Overview {#sec-jags-overview}

In this chapter, we will be working up to simple univariate state-space JAGS models. We will write each of the models with the same univariate state-space form.

\begin{equation}
\begin{gathered}
x_t = b x_{t-1} + u + w_t, \, w_t \sim N(0,q)\\
y_t = x_t + a + v_t, \, v_t \sim \,\text{N}(0, r)
\end{gathered}
(\#eq:jags-uniss)
\end{equation}

We will be fitting linear regressions with this form, and this will mean the JAGS code is more verbose than necessary, but the goal is to build up to our univariate state-space code by building off simpler models.

## Univariatate response models {#sec-jags-univariate}

### Linear regression with no covariates {#sec-jags-lr-no-covariates}

We will start with a linear regression with only an intercept. We will write the model in the form of Equation \@ref(eq:jags-uniss). Our model is
\begin{equation}
\begin{gathered}
x_t = u \\
y_t = x_t + v_t, v_t \sim \,\text{N}(0, r)
\end{gathered}
(\#eq:jags-lr1)
\end{equation}
An equivalent way to think about this model is 
\begin{equation}
Y_t \sim \,\text{N}(E[Y_t], r)
\end{equation}
$E[Y_{t}] = x_t$ where $x_t = u$. 
In this linear regression model, we will treat the residual error as independent and identically distributed Gaussian observation error.

To run the JAGS model, we will need to start by writing the model in JAGS notation.  We can construct the model in Equation \@ref(eq:jags-lr1) as

```{r jags-lr1, results='hide', cache=TRUE}
# LINEAR REGRESSION intercept only.

model.loc <- "lm_intercept.txt" # name of the txt file
jagsscript <- cat("
model {  
   # priors on parameters
   u ~ dnorm(0, 0.01);
   inv.r ~ dgamma(0.001,0.001); # This is inverse gamma
   r <- 1/inv.r; # derived value
   
   # likelihood
   for(i in 1:N) {
      X[i] <- u
      EY[i] <- X[i]; # derived value
      Y[i] ~ dnorm(EY[i], inv.r); 
   }
}  
", file = model.loc)
```
```{r clean1, echo=FALSE, purl=FALSE}
removefiles <- c(removefiles, model.loc)
```

The JAGS code has three parts: our parameter priors, our data model and derived parameters.

**Parameter priors** There are two parameters in the model ($u$, the mean, and $r$, the variance of the observation error). We need to set a prior on both of these. We will set a vague prior of a Gaussian with varianc 10 on $u$. In JAGS instead of specifying the normal distribution with the variance, $N(0, 10)$, you specify it with the precision (1/variance), so our prior on $u$ is  `dnorm(0, 0.01)`. For $r$, we need to set a prior on the precision $1/r$, which we call `inv.r` in the code. The precision receives a gamma prior, which is equivalent to the variance receiving an inverse gamma prior (fairly common for standard Bayesian regression models). 

**Likelihood** Our data distribution is $Y_t \sim \,\text{N}(E[Y_t], r)$. We use the `dnorm()` distribution with the precision ($1/r$) instead of $r$. So our data model is `Y[t] = dnorm(EY[t], inv.r)`. JAGS is not vectorized so we need to use for loops (instead of matrix multiplication) and use the for loop to specify the distribution for each `Y[t]`. For, this model we didn't actually need `X[t]` but we use it because we are building up to a state-space model which has both $x_t$ and $y_t$.

**Derived values** Derived values are things we want output so we can track them. In this example, our derived values are a bit useless but in more complex models they will be quite handy. Also they can make your code easier to understand.


To run the model, we need to create several new objects, representing (1) a list of data that we will pass to JAGS `jags.data`, (2) a vector of parameters that we want to monitor and have returned back to R `jags.params`, and (3) the name of our text file that contains the JAGS model we wrote above. With those three things, we can call the `jags()` function. 

```{r jags-call-fun1, results='hide', cache=TRUE}
jags.data <- list("Y" = Wind, "N" = N) 
jags.params <- c("r", "u") # parameters to be monitored
mod_lm_intercept <- R2jags::jags(jags.data,
  parameters.to.save = jags.params,
  model.file = model.loc, n.chains = 3, n.burnin = 5000,
  n.thin = 1, n.iter = 10000, DIC = TRUE
)
```

The function from the **R2jags** package that we use to run the model is `jags()`. There is a parallel version of the function called `jags.parallel()` which is useful for larger, more complex models. The details of both can be found with `?jags` or `?jags.parallel`.

Notice that the `jags()` function contains a number of other important arguments. In general, larger is better for all arguments: we want to run multiple MCMC chains (maybe 3 or more), and have a burn-in of at least 5000. The total number of samples after the burn-in period is n.iter-n.burnin, which in this case is 5000 samples. Because we are doing this with 3 MCMC chains, and the thinning rate equals 1 (meaning we are saving every sample), we will retain a total of 1500 posterior samples for each parameter.  

The saved object storing our model diagnostics can be accessed directly, and includes some useful summary output.
```{r jags-lm1-mod}
mod_lm_intercept
```

The last two columns in the summary contain `Rhat` (which we want to be close to 1.0), and `neff` (the effective sample size of each set of posterior draws). To examine the output more closely, we can pull all of the results directly into R,

```{r jags-lm1-attach, eval=FALSE}
R2jags::attach.jags(mod_lm_intercept)
```
Attaching the **R2jags** object loads the posteriors for the parameters and we can call them directly, e.g. `u`. If we don't want to attach them to our workspace, we can find the posteriors within the model object.
```{r jags-lm1-attach2}
post.params <- mod_lm_intercept$BUGSoutput$sims.list
```

We make a histogram of the posterior distributions of the parameters `u` and `r` with the following code,

```{r jags-plot-lm1, echo=TRUE, eval=TRUE, fig.show='hide'}
# Now we can make plots of posterior values
par(mfrow = c(2, 1))
hist(post.params$u, 40, col = "grey", xlab = "u", main = "")
hist(post.params$r, 40, col = "grey", xlab = "r", main = "")
```

(ref:jags-plot-hist-post) Plot of the posteriors for the linear regression model.

```{r jags-plot-hist-post, fig=TRUE, echo=FALSE, fig.width=6, fig.height=6, fig.cap='(ref:jags-plot-hist-post)'}
par(mfrow = c(2, 1))
hist(post.params$u, 40, col = "grey", xlab = "u", main = "")
hist(post.params$r, 40, col = "grey", xlab = "r", main = "")
```
 
We can run some useful diagnostics from the **coda** package on this model output. We have written a small function to make the creation of a MCMC list (an argument required for many of the diagnostics). The function is

```{r jags-lm1-mcmclist-func, cache=TRUE}
createMcmcList <- function(jagsmodel) {
  McmcArray <- as.array(jagsmodel$BUGSoutput$sims.array)
  McmcList <- vector("list", length = dim(McmcArray)[2])
  for (i in 1:length(McmcList)) McmcList[[i]] <- as.mcmc(McmcArray[, i, ])
  McmcList <- mcmc.list(McmcList)
  return(McmcList)
}
```

Creating the  MCMC list preserves the random samples generated from each chain and allows you to extract the samples for a given parameter (such as $\mu$) from any chain you want. To extract $\mu$ from the first chain, for example, you could use the following code. Because `createMcmcList()` returns a list of **mcmc** objects, we can summarize and plot these directly. Figure \@ref(fig:jags-plot-myList) shows the plot from `plot(myList[[1]])`.

```{r jags-make-myList, fig.show='hide'}
myList <- createMcmcList(mod_lm_intercept)
summary(myList[[1]])
plot(myList[[1]])
```

(ref:jags-plot-myList) Plot of an object output from $\texttt{creatMcmcList}$.

```{r jags-plot-myList,fig=TRUE, echo=FALSE, fig.width=6, fig.height=6, fig.cap='(ref:jags-plot-myList)'}
plot(myList[[1]])
```
 
For more quantitative diagnostics of MCMC convergence, we can rely on the **coda** package in R. There
are several useful statistics available, including the Gelman-Rubin diagnostic (for one or several chains), autocorrelation diagnostics (similar to the ACF you calculated above), the Geweke diagnostic, and Heidelberger-Welch test of stationarity. 

```{r jags-coda, results='hide', message=FALSE, warning=FALSE}
library(coda)
gelmanDiags <- coda::gelman.diag(createMcmcList(mod_lm_intercept), multivariate = FALSE)
autocorDiags <- coda::autocorr.diag(createMcmcList(mod_lm_intercept))
gewekeDiags <- coda::geweke.diag(createMcmcList(mod_lm_intercept))
heidelDiags <- coda::heidel.diag(createMcmcList(mod_lm_intercept))
```


### Linear regression with covariates {#sec-jags-covariates}

We can introduce `Temp` as the covariate explaining our response variable `Wind`. Our new equation is

\begin{equation}
\begin{gathered}
x_t = u + C\,c_t\\
y_t = x_t + v_t, v_t \sim \,\text{N}(0, r)
\end{gathered}
(\#eq:jags-obs-eqn)
\end{equation}

To create JAGS code for this model, we (1) add a prior for our new parameter  `C`, (2) update `X[i]` equation to include the new covariate, and (3) we include the new covariate in our named data list. 

```{r jags-cov, results='hide', cache=TRUE}
# 1. LINEAR REGRESSION with covariates

model.loc <- ("lm_covariate.txt")
jagsscript <- cat("
model {  
   # priors on parameters
   u ~ dnorm(0, 0.01); 
   C ~ dnorm(0,0.01);
   inv.r ~ dgamma(0.001,0.001); 
   r <- 1/inv.r; 
   
   # likelihood
   for(i in 1:N) {
      X[i] <- u + C*c[i];
      EY[i] <- X[i]
      Y[i] ~ dnorm(EY[i], inv.r);
   }
}  
", file = model.loc)

jags.data <- list("Y" = Wind, "N" = N, "c" = Temp)
jags.params <- c("r", "EY", "u", "C")
mod_lm <- R2jags::jags(jags.data,
  parameters.to.save = jags.params,
  model.file = model.loc, n.chains = 3, n.burnin = 5000,
  n.thin = 1, n.iter = 10000, DIC = TRUE
)
```
```{r clean6, echo=FALSE, purl=FALSE}
removefiles <- c(removefiles, model.loc)
```


We can show the the posterior fits (the model fits) to the data. Here is a simple function whose arguments are one of our fitted models and the raw data. The function is: 

```{r jags-lm1ar-plot-func, results='hide'}
plotModelOutput <- function(jagsmodel, Y) {
  # attach the model
  EY <- jagsmodel$BUGSoutput$sims.list$EY
  x <- seq(1, length(Y))
  summaryPredictions <- cbind(
    apply(EY, 2, quantile, 0.025), apply(EY, 2, mean),
    apply(EY, 2, quantile, 0.975)
  )
  plot(Y,
    col = "white", ylim = c(min(c(Y, summaryPredictions)), max(c(Y, summaryPredictions))),
    xlab = "", ylab = "95% CIs of predictions and data", main = paste(
      "JAGS results:",
      jagsmodel$model.file
    )
  )
  polygon(c(x, rev(x)), c(summaryPredictions[, 1], rev(summaryPredictions[, 3])),
    col = "grey70", border = NA
  )
  lines(summaryPredictions[, 2])
  points(Y)
}
```

We can use the function to plot the predicted posterior mean with 95\% CIs, as well as the raw data. Note that the shading is for the CIs on the expected value of $y_t$ so will look narrow relative to the data. For example, try
```{r jags-lm1ar-plot, results='hide', eval=FALSE, fig.show='hide'}
plotModelOutput(mod_lm, Wind)
```


(ref:jags-lm1ar-plot1) Predicted posterior mean with 95\% CIs 

```{r jags-lm1ar-plot1, fig=TRUE, echo=FALSE, fig.width=6, fig.height=6, fig.cap='(ref:jags-lm1ar-plot1)', cache=FALSE}
plotModelOutput(mod_lm, Wind)
```


### Random walk with drift {#sec-jags-rw}

The previous models were observation error only models. Switching gears, we can create process error models. We will start with a random walk model. In this model, the assumption is that the underlying state $x_t$ is measured perfectly. All stochasticity is originating from process variation: variation in $x_t$ to $x_{t+1}$. 

For this simple model, we will assume that wind behaves as a random walk. We will call this process $x$ to prepare for the state-space model to come. We have no $y_t$ part of the equation in this model.

\begin{equation}
x_t = x_{t-1} + u + w_t, \text{ where }w_t \sim \,\text{N}(0,q)
(\#eq:jags-obs-eqn-rw)
\end{equation}
Now $x_t$ is stochastic and $E[X_t] = x_{t-1} + u$ and $X_t \sim \,\text{N}(E[X_t],q)$.

We are going to need to put a prior on $x_0$, which appears in $E[X_1]$. We could start with $t=2$ and skip this but we will start at $t=1$ since we will need to do that for later problems. The question is what prior should we put on $x_0$? This is not a stationary process. We will just put a vague prior on $x_0$.

The JAGS random walk model is:
```{r jags-rw, results='hide', cache=TRUE}
# RANDOM WALK with drift

model.loc <- ("rw_intercept.txt")
jagsscript <- cat("
model {  
   # priors on parameters
   u ~ dnorm(0, 0.01); 
   inv.q ~ dgamma(0.001,0.001); 
   q <- 1/inv.q;
   X0 ~ dnorm(0, 0.001);
   
   # likelihood
   X[1] ~ dnorm(X0 + u, inv.q);
   for(i in 2:N) {
      X[i] ~ dnorm(X[i-1] + u, inv.q);
   }
}  
", file = model.loc)
```
```{r clean3, echo=FALSE, purl=FALSE}
removefiles <- c(removefiles, model.loc)
```

To fit this model, we need to change `jags.data` to pass in `X = Wind` instead of `Y = Wind`. Obvioously we could have written the JAGS code with `Y` in place of `X` and kept our `jags.data` code the same as before, but we are working up to a state-space model where we have a hidden random walk called `X` and an observation of that called `Y`. 

```{r jags-rw-fit, results='hide', cache=TRUE}
jags.data <- list("X" = Wind, "N" = N)
jags.params <- c("q", "u")
mod_rw_intercept <- R2jags::jags(jags.data,
  parameters.to.save = jags.params, model.file = model.loc,
  n.chains = 3, n.burnin = 5000, n.thin = 1, n.iter = 10000, DIC = TRUE
)
```

### Autoregressive AR(1) time series models {#sec-jags-ar1}

A variation of the random walk model is the autoregressive time series model of order 1, AR(1). This model introduces a coefficient, which we will call $b$. The parameter $b$ controls the degree to which the random walk reverts to the mean. When $b = 1$, the model is identical to the random walk, but at smaller values, the model will revert back to the mean (which in this case is zero). Also, $b$ can take on negative values.

\begin{equation}
x_t = b \, x_{t-1} + u + w_t, \text{ where }w_t \sim \,\text{N}(0,q)
(\#eq:jags-obs-eqn-ar1)
\end{equation}
Now $E[X_t] = b \, x_{t-1} + u$.

Once again we need to put a prior on $x_0$, which appears in $E[X_1]$. An AR(1) with $|b|<1$ is a stationary process and the variance of the stationary distribution of $x_t$ is $q/(1-b^2)$. If you think that $x_0$ has the stationary distribution (does your data look stationary?) then you can use the variance of the stationary distribution of $x_t$ for your prior. We specify priors with the precision (1 over the variance) instead of the variance. Thus the precision of the stationary distribution of $x_0$ is $(1/q)(1-b^2)$. In the code, `inv.q` is $1/q$ and the precision is `inv.q * (1-b*b)`.

```{r jags-ar1est, echo=TRUE, results='hide', cache=TRUE}
# AR(1) MODEL WITH AND ESTIMATED AR COEFFICIENT

model.loc <- ("ar1_intercept.txt")
jagsscript <- cat("
model {  
   # priors on parameters
   u ~ dnorm(0, 0.01); 
   inv.q ~ dgamma(0.001,0.001); 
   q <- 1/inv.q; 
   b ~ dunif(-1,1);
   X0 ~ dnorm(0, inv.q * (1 - b * b));
   
   # likelihood
   X[1] ~ dnorm(b * X0 + u, inv.q);
   for(t in 2:N) {
      X[t] ~ dnorm(b * X[t-1] + u, inv.q);
   }
}  
", file = model.loc)

jags.data <- list("X" = Wind, "N" = N)
jags.params <- c("q", "u", "b")
mod_ar1_intercept <- R2jags::jags(jags.data,
  parameters.to.save = jags.params,
  model.file = model.loc, n.chains = 3, n.burnin = 5000, n.thin = 1,
  n.iter = 10000, DIC = TRUE
)
```
```{r clean4, echo=FALSE, purl=FALSE}
removefiles <- c(removefiles, model.loc)
```

### Regression with AR(1) errors {#sec-jags-regression-with-corr-errors2}

The AR(1) model in the previous section suggests a way that we could include correlated errors in our linear regression. We could use the $x_t$ AR(1) process as our errors for $y_t$. Here is an example of modifying the intercept only linear regression model. We will set $u$ to 0 so that our AR(1) errors have a mean of 0.

\begin{equation}
\begin{gathered}
x_t = b \, x_{t-1} + w_t, \text{ where }w_t \sim \,\text{N}(0,q) \\
y_t = a + x_t
\end{gathered}
(\#eq:jags-obs-eqn-lmar2)
\end{equation}

The problem with this is that we need a distribution for $y_t$. We cannot use `Y[t] <- a + X[t]` in our JAGS code ($Y_t$ is a random varible with a distribution; you cannot assign it a value). We need to re-write this as $Y_t \sim N(a + b \, x_{t-1}, q)$.

\begin{equation}
\begin{gathered}
Y_t \sim N(a + b \, x_{t-1}, q) \\
x_t = y_t - a
\end{gathered}
(\#eq:jags-obs-eqn-lmar2b)
\end{equation}

We will create the variable `EY` so we can keep track of the expected value of $Y_t$, conditioned on $t-1$.

```{r jags-lm1ar2, results='hide', cache=TRUE}
# LINEAR REGRESSION with autocorrelated errors
# no covariates, intercept only.

model.loc <- ("lm_intercept_ar1b.txt")
jagsscript <- cat("
model {  
   # priors on parameters
   a ~ dnorm(0, 0.01); 
   inv.q ~ dgamma(0.001,0.001); 
   q <- 1/inv.q; 
   b ~ dunif(-1,1);
   X0 ~ dnorm(0, inv.q * (1 - b * b));
   
   # likelihood
   EY[1] <- a + b * X0;
   Y[1] ~ dnorm(EY[1], inv.q);
   X[1] <- Y[1] - a;
   for(t in 2:N) {
      EY[t] <- a + b * X[t-1];
      Y[t] ~ dnorm(EY[1], inv.q);
      X[t] <- Y[t]-a;
   }
}  
", file = model.loc)

jags.data <- list("Y" = Wind, "N" = N)
jags.params <- c("q", "EY", "a", "b")
mod_ar1_intercept <- R2jags::jags(jags.data,
  parameters.to.save = jags.params,
  model.file = model.loc, n.chains = 3, n.burnin = 5000, n.thin = 1,
  n.iter = 10000, DIC = TRUE
)
```
```{r clean2, include=FALSE, purl=FALSE}
removefiles <- c(removefiles, model.loc)
```



### Univariate state space model {#sec-jags-uss}

Now we will combine the process and observation models to create a univariate state-space model. This is the classic stochastic level model.

\begin{equation}
\begin{gathered}
x_t = x_{t-1} + u + w_t, \, w_t \sim N(0,q)\\
y_t = x_t + v_t, \, v_t \sim \,\text{N}(0, r)
\end{gathered}
(\#eq:jags-proc-eqn-uss)
\end{equation}

Because $x$ is a random walk model not a stationary AR(1), we will place a vague weakly informative prior on $x_0$: $x_0 \sim \,\text{N}(y_1, 1000)$. We had to pass in `Y1` as data because JAGS would complain if we used `Y[1]` in our prior (because have `X0` in our model for $Y[1]$). `EY` is added so that we can track the model fits for $y$. In this case it is just `X` but in more complex models it will involve more parameters.

```{r jags-ss1, echo=TRUE, results='hide', cache=TRUE}
model.loc <- ("ss_model.txt")
jagsscript <- cat("
model {  
   # priors on parameters
   u ~ dnorm(0, 0.01); 
   inv.q ~ dgamma(0.001,0.001); 
   q <- 1/inv.q;
   inv.r ~ dgamma(0.001,0.001);
   r <- 1/inv.r; 
   X0 ~ dnorm(Y1, 0.001);
   
   # likelihood
   X[1] ~ dnorm(X0 + u, inv.q);
   EY[1] <- X[1];
   Y[1] ~ dnorm(EY[1], inv.r);
   for(t in 2:N) {
      X[t] ~ dnorm(X[t-1] + u, inv.q);
      EY[t] <- X[t];
      Y[t] ~ dnorm(EY[t], inv.r); 
   }
}  
", file = model.loc)
```

We fit as usual with the addition of `Y1` in `jags.data`.

```{r jags-ss1-fit, echo=TRUE, results='hide', cache=TRUE}
jags.data <- list("Y" = Wind, "N" = N, Y1 = Wind[1])
jags.params <- c("q", "r", "EY", "u")
mod_ss <- jags(jags.data,
  parameters.to.save = jags.params, model.file = model.loc, n.chains = 3,
  n.burnin = 5000, n.thin = 1, n.iter = 10000, DIC = TRUE
)
```

```{r clean5, echo=FALSE, purl=FALSE}
removefiles <- c(removefiles, model.loc)
```

## Multivariate state-space models {#sec-jags-marss}

In the multivariate state-space model, our observations and hidden states can be multivariate along with all the parameters:
\begin{equation}
\begin{gathered}
\mathbf{x}_t = \mathbf{B} \mathbf{x}_{t-1}+\mathbf{u}+\mathbf{w}_t \text{ where } \mathbf{w}_t \sim \,\text{N}(0,\mathbf{Q}) \\
\mathbf{y}_t = \mathbf{Z}\mathbf{x}_t+\mathbf{a}+\mathbf{v}_t \text{ where } \mathbf{v}_t \sim \,\text{N}(0,\mathbf{R}) \\
\mathbf{x}_0 = \boldsymbol{\mu}
\end{gathered}   
(\#eq:jags-marss)
\end{equation}

### One hidden state

Let's start with a very simple MARSS model with JAGS: two observation time-series and one hidden state. Our $\mathbf{x}_t$ model is $x_t = x_{t-1} + u + w_t$ and our $\mathbf{y}_t$ model is 
\begin{equation}
 \begin{bmatrix}
    y_{1} \\
    y_{2}\end{bmatrix}_t = 
    \begin{bmatrix}
    1\\
    1\end{bmatrix} x_t +  
    \begin{bmatrix}
    0 \\
    a_2\end{bmatrix} + 
    \begin{bmatrix}
    v_{1} \\
    v_{2}\end{bmatrix}_t, \,
    \begin{bmatrix}
    v_{1} \\
    v_{2}\end{bmatrix}_t \sim 
    \,\text{MVN}\left(0, \begin{bmatrix}
    r_1&0 \\
    0&r_2\end{bmatrix}\right)
(\#eq:jags-marss-y)
\end{equation}

We need to put a prior on our $x_0$ (initial $x$). Since $b=1$, we have a random walk rather than a stationary process and we will put a vague prior on the $x_0$. We need to deal with the $\mathbf{a}$ so that our code doesn't run in circles by trying to match $x$ up with different $y_t$ time series. We force $x_t$ to track the mean of $y_{1,t}$ and then use $a_2$ to scale the other $y_t$ relative to that. The problem is that a random walk is very flexible and if we tried to estimate $a_1$ then we would have infinite solutions.

To keep our JAGS code organized, let's separate the $\mathbf{x}$ and $\mathbf{y}$ parts of the code.
```{r jags-jagsscript-marss1}
jagsscript <- cat("
model {  
   # process model priors
   u ~ dnorm(0, 0.01); # one u
   inv.q~dgamma(0.001,0.001);
   q <- 1/inv.q; # one q
   X0 ~ dnorm(Y1,0.001); # initial state
   # process model likelihood
   EX[1] <- X0 + u;
   X[1] ~ dnorm(EX[1], inv.q);
   for(t in 2:N) {
         EX[t] <- X[t-1] + u;
         X[t] ~ dnorm(EX[t], inv.q);
   }

   # observation model priors
   for(i in 1:n) { # r's differ by site
     inv.r[i]~dgamma(0.001,0.001);
     r[i] <- 1/inv.r[i];
   }
   a[1] <- 0; # first a is 0, rest estimated
   for(i in 2:n) {
     a[i]~dnorm(0,0.001);
   }   
   # observation model likelihood
   for(t in 1:N) {
     for(i in 1:n) {
       EY[i,t] <- X[t]+a[i]
       Y[i,t] ~ dnorm(EY[i,t], inv.r[i]);
     }
   }
}  

",file="marss-jags1.txt")
```
```{r clean-marss1, echo=FALSE, purl=FALSE}
removefiles <- c(removefiles, model.loc)
```

To fit the model, we write the data list, parameter list, and pass the model to the `jags()` function.
```{r jags-marss1-fit, results='hide', message=FALSE, cache=TRUE}
data(harborSealWA, package="MARSS")
dat <- t(harborSealWA[,2:3])
jags.data <- list("Y" = dat, n = nrow(dat), N = ncol(dat), Y1 = dat[1,1]) 
jags.params <- c("EY", "u", "q", "r")
model.loc <- "marss-jags1.txt" # name of the txt file
mod_marss1 <- R2jags::jags(jags.data,
  parameters.to.save = jags.params,
  model.file = model.loc, n.chains = 3,
  n.burnin = 5000, n.thin = 1, n.iter = 10000, DIC = TRUE
)
```

We can make a plot of our estimated parameters:
```{r jags-marss1-hist}
post.params <- mod_marss1$BUGSoutput$sims.list
par(mfrow=c(2,2))
hist(log(post.params$q), main="log(q)", xlab="")
hist(post.params$u, main="u", xlab="")
hist(log(post.params$r[,1]), main="log(r_1)", xlab="")
hist(log(post.params$r[,2]), main="log(r_2)", xlab="")
```

We can make a plot of the model fitted $y_t$ with 50% credible intervals and the data. Note that the credible intervals are for the expected value of $y_{i,t}$ so will be narrower than the data.
```{r jags-marss1-plot-fun, warning=FALSE, message=FALSE}
make.ey.plot <- function(mod, dat){
   library(ggplot2)
EY <- mod$BUGSoutput$sims.list$EY
n <- nrow(dat); N <- ncol(dat)
df <- c()
for(i in 1:n){
tmp <- data.frame(n = paste0("Y",i),
                  x = 1:N, 
                  ey=apply(EY[,i,, drop=FALSE],3,median),
                  ey.low=apply(EY[,i,, drop=FALSE],3,quantile,probs=0.25),
                  ey.up=apply(EY[,i,, drop=FALSE],3,quantile,probs=0.75),
                  y=dat[i,]
                  )
df <- rbind(df, tmp)
}
ggplot(df, aes(x=x, y=ey)) + geom_line() +
   geom_ribbon(aes(ymin=ey.low, ymax=ey.up), alpha=0.25) +
   geom_point(data=df, aes(x=x, y=y)) +
   facet_wrap(~n) + theme_bw()
}
```

```{r jags-marss1-plot, warning=FALSE, message=FALSE}
make.ey.plot(mod_marss1, dat)
```


### $m$ hidden states

Let's add multiple hidden states. We'll say that each $y_t$ is observing its own $x_t$ but the $x_t$ share the same $q$ but not $u$. Our $\mathbf{x}_t$ model is \begin{equation}
 \begin{bmatrix}
    x_{1} \\
    x_{2}\end{bmatrix}_t = 
    \begin{bmatrix}
    1&0\\
    0&1\end{bmatrix}
    \begin{bmatrix}
    x_{1} \\
    x_{2}\end{bmatrix}_{t-1} +  
    \begin{bmatrix}
    u_1 \\
    u_2\end{bmatrix} + 
    \begin{bmatrix}
    w_{1} \\
    w_{2}\end{bmatrix}_t, \,
    \begin{bmatrix}
    w_{1} \\
    w_{2}\end{bmatrix}_t \sim 
    \,\text{MVN}\left(0, \begin{bmatrix}
    q&0 \\
    0&q\end{bmatrix}\right)
(\#eq:jags-marss2)
\end{equation}

Here is the JAGS model. Note that $a_i$ is 0 for all $i$ because each $y_t$ is associated with its own $x_t$.
```{r jags-jagsscript-marss2}
jagsscript <- cat("
model {  
   # process model priors
   inv.q~dgamma(0.001,0.001);
   q <- 1/inv.q; # one q
   for(i in 1:n) {
      u[i] ~ dnorm(0, 0.01); 
      X0[i] ~ dnorm(Y1[i],0.001); # initial states
   }
   # process model likelihood
   for(i in 1:n) {
     EX[i,1] <- X0[i] + u[i];
     X[i,1] ~ dnorm(EX[i,1], inv.q);
   }
   for(t in 2:N) {
      for(i in 1:n) {
         EX[i,t] <- X[i,t-1] + u[i];
         X[i,t] ~ dnorm(EX[i,t], inv.q);
      }
   }

   # observation model priors
   for(i in 1:n) { # The r's are different by site
     inv.r[i]~dgamma(0.001,0.001);
     r[i] <- 1/inv.r[i];
   }
   # observation model likelihood
   for(t in 1:N) {
     for(i in 1:n) {
       EY[i,t] <- X[i,t]
       Y[i,t] ~ dnorm(EY[i,t], inv.r[i]);
     }
   }
}  

",file="marss-jags2.txt")
```
```{r clean-marss2, echo=FALSE, purl=FALSE}
removefiles <- c(removefiles, model.loc)
```

Our code to fit the model changes a little.
```{r jags-marss2-fit, results='hide', message=FALSE, cache=TRUE}
data(harborSealWA, package="MARSS")
dat <- t(harborSealWA[,2:3])
jags.data <- list("Y" = dat, n = nrow(dat), N = ncol(dat), Y1 = dat[,1]) 
jags.params <- c("EY", "u", "q", "r")
model.loc <- "marss-jags2.txt" # name of the txt file
mod_marss1 <- R2jags::jags(jags.data,
  parameters.to.save = jags.params,
  model.file = model.loc, n.chains = 3,
  n.burnin = 5000, n.thin = 1, n.iter = 10000, DIC = TRUE
)
```

```{r jags-marss2-plot, warning=FALSE, message=FALSE}
make.ey.plot(mod_marss1, dat)
```

## Non-Gaussian observation errors {#sec-jags-non-gaussian}

### Poisson observation errors {#sec-jags-poisson}

So far we have used the following observation model $y_t \sim \,\text{N}(x_t, r)$. 
We can change this to a Poisson observation error model:
$Y_t \sim \text{Pois}(\lambda_t)$ where $E[Y_t] = \lambda_t$. $\text{log}(\lambda_t) = x_t$ where $x_t$ is our process model. 

All we need to change to allow Poisson errors is to change the `Y[t]` part to
```
log(EY[t]) <- X[i]
Y[t] ~ dpois(EY[t])
```
We also need to ensure that our data are 
integers and we remove the `r` part from our model code since the Poisson does not have that. 

Our univariate state-space code with Poisson observation errors is the following:

```{r jags-ss1-pois, echo=TRUE, results='hide'}
# SS MODEL with Poisson errors

model.loc <- ("ss_model_pois.txt")
jagsscript <- cat("
model {  
   # priors on parameters
   u ~ dnorm(0, 0.01); 
   inv.q ~ dgamma(0.001,0.001); 
   q <- 1/inv.q;
   X0 ~ dnorm(0, 0.001);
   
   # likelihood
   X[1] ~ dnorm(X0 + u, inv.q);
   log(EY[1]) <- X[1]
   Y[1] ~ dpois(EY[1])
   for(t in 2:N) {
      X[t] ~ dnorm(X[t-1] + u, inv.q);
      log(EY[t]) <- X[t]
      Y[t] ~ dpois(EY[t]); 
   }
}  
", file = model.loc)
```

We will fit this to the wild dogs data in the **MARSS** package. 

```{r jags-ss1-fit-pois, echo=TRUE, results='hide'}
data(wilddogs, package="MARSS")
jags.data <- list("Y" = wilddogs[,2], "N" = nrow(wilddogs))
jags.params <- c("q", "EY", "u")
mod_ss <- jags(jags.data,
  parameters.to.save = jags.params, model.file = model.loc, n.chains = 3,
  n.burnin = 5000, n.thin = 1, n.iter = 10000, DIC = TRUE
)
```

```{r clean7, echo=FALSE, purl=FALSE}
removefiles <- c(removefiles, model.loc)
```

When we use this univariate state-space model with population data, like the wild dogs, we would log the data$^\dagger$, and our $y_t$ in our JAGS code is really $log(y_t)$. In that case, $E[log(Y_t)]) = f(x_t)$. So there is a log-link that we are not really explicit about when we pass in the log of our data. In the Poisson model, that log relationship is explicit, aka we specify $log(E[Y_t]) = x_t$ and we pass in the raw count data not the log of the data. 

$\dagger$ Why would we typically log population data in this case? Because we would typically think of population processes as multiplicative. Population size at time $t$ is growth **times** population size at time $t-1$. By logging the data, we convert to an additive process. Log population size at time $t$ is log growth **plus** log population size at time $t-1$.

### Negative binomial observation errors {#sec-jags-negbin}

In the Poisson distribution, the mean and variance are the same. Using the negative binomial distribution, we can relax that assumption and allow the mean and variance to be different. The negative binomial distribution has two parameters, $r$ and $p$. $r$ is the dispersion parameter. As $r \rightarrow \infty$, the distribution becomes the Poisson distribution and when $r$ is small, the distribution is overdispersed (higher variance) relative to the Poisson. In practice, $r > 30$ is going to be very close to the Poisson. $p$ is the success parameter, $p = r/(r+E[Y_t])$. As for the Poisson, $log E[Y_{t}] = x_t$---for the univariate state-space model in this example with one state, $z=1$ and $a=0$.

To allow negative binomial errors we change the `Y[t]` part to
```
log(EY[t]) <- X[t]
p[t] <- r/(r + EY[t])
Y[t] ~ dnegbin(p[t], r)
```
Now that we have $r$ again in the model, we will need to put a prior on it. $r$ is positive and 50 is close to infinity. The following is a sufficiently vague prior.
```
r ~ dunif(0,50)
```

Our univariate state-space code with negative binomial observation errors is the following:

```{r jags-ss1-negbin, echo=TRUE, results='hide'}
# SS MODEL with negative binomial errors

model.loc <- ("ss_model_negbin.txt")
jagsscript <- cat("
model {  
   # priors on parameters
   u ~ dnorm(0, 0.01); 
   inv.q ~ dgamma(0.001,0.001); 
   q <- 1/inv.q;
   r ~ dunif(0,50);
   X0 ~ dnorm(0, 0.001);
   
   # likelihood
   X[1] ~ dnorm(X0 + u, inv.q);
   log(EY[1]) <- X[1]
   p[1] <- r/(r + EY[1])
   Y[1] ~ dnegbin(p[1], r)
   for(t in 2:N) {
      X[t] ~ dnorm(X[t-1] + u, inv.q);
      log(EY[t]) <- X[t]
      p[t] <- r/(r + EY[t])
      Y[t] ~ dnegbin(p[t], r)
   }
}  
", file = model.loc)
```

We will fit this to the wild dogs data in the **MARSS** package. 

```{r jags-ss1-fit-negbin, echo=TRUE, results='hide'}
data(wilddogs, package="MARSS")
jags.data <- list("Y" = wilddogs[,2], "N" = nrow(wilddogs))
jags.params <- c("q", "EY", "u", "r")
mod_ss <- jags(jags.data,
  parameters.to.save = jags.params, model.file = model.loc, n.chains = 3,
  n.burnin = 5000, n.thin = 1, n.iter = 10000, DIC = TRUE
)
```

```{r clean-negbin, echo=FALSE, purl=FALSE}
removefiles <- c(removefiles, model.loc)
```

## Forecasting with JAGS models {#sec-jags-forecast}

There are a number of different approaches to using Bayesian time series models to perform forecasting. One approach might be to fit a model, and use those posterior distributions to forecast as a secondary step (say within R). A more streamlined approach is to do this within the JAGS code itself. We can take advantage of the fact that JAGS allows you to include NAs in the response variable (but never in the predictors). Let's use the same Wind dataset, and the univariate state-space model described above to forecast three time steps into the future. We can do this by including 3 more NAs in the dataset, and incrementing the variable ```N``` by 3.

```{r jags-cov-forecast, results='hide', cache=TRUE}
jags.data <- list("Y" = c(Wind, NA, NA, NA), "N" = (N + 3), Y1 = Wind[1])
jags.params <- c("q", "r", "EY", "u")
model.loc <- ("ss_model.txt")
mod_ss_forecast <- jags(jags.data,
  parameters.to.save = jags.params,
  model.file = model.loc, n.chains = 3, n.burnin = 5000, n.thin = 1,
  n.iter = 10000, DIC = TRUE
)
```

We can inspect the fitted model object, and see that ```EY``` contains the 3 new predictions for the forecasts from this model. 

\clearpage

## Problems {#sec-jags-problems}


1. Fit the intercept only model from section \@ref(sec-jags-lr-no-covariates). Set the burn-in to 3, and when the  model completes, plot the time series of the parameter ```u``` for the first MCMC chain.

    a. Based on your visual inspection, has the MCMC chain convered?
    
    b. What is the ACF of the first MCMC chain?

2. Increase the MCMC burn-in for the model in question 1 to a value that you think is reasonable. After the model has converged, calculate the Gelman-Rubin diagnostic for the fitted model object.

3. Compare the results of the ```plotModelOutput()``` function for the intercept only model from section \@ref(sec-jags-lr-no-covariates). You will to add "predY" to your JAGS model and to the list of parameters to monitor, and re-run the model.

5. Plot the posterior distribution of $b$ for the AR(1) model in section \@ref(sec-jags-ar1).  Can this parameter be well estimated for this dataset?

6. Plot the posteriors for the process and observation variances (not standard deviation) for the univariate state-space model in section \@ref(sec-jags-uss). Which is larger for this dataset?

7. Add the effect of temperature to the AR(1) model in section \@ref(sec-jags-ar1).  Plot the posterior for ```C``` and compare to the posterior for ```C``` from the model in section \@ref(sec-jags-covariates).

8. Plot the fitted values from the model in section \@ref(sec-jags-forecast), including the forecasts, with the 95\% credible intervals for each data point.

9. The following is a dataset from the Upper Skagit River (Puget Sound, 1952-2005) on salmon spawners and recruits:
    ```{r jags-hwdata, echo=TRUE}
Spawners <- c(2662, 1806, 1707, 1339, 1686, 2220, 3121, 5028, 9263, 4567, 1850, 3353, 2836, 3961, 4624, 3262, 3898, 3039, 5966, 5931, 7346, 4911, 3116, 3185, 5590, 2485, 2987, 3829, 4921, 2348, 1932, 3151, 2306, 1686, 4584, 2635, 2339, 1454, 3705, 1510, 1331, 942, 884, 666, 1521, 409, 2388, 1043, 3262, 2606, 4866, 1161, 3070, 3320)
Recruits <- c(12741, 15618, 23675, 37710, 62260, 32725, 8659, 28101, 17054, 29885, 33047, 20059, 35192, 11006, 48154, 35829, 46231, 32405, 20782, 21340, 58392, 21553, 27528, 28246, 35163, 15419, 16276, 32946, 11075, 16909, 22359, 8022, 16445, 2912, 17642, 2929, 7554, 3047, 3488, 577, 4511, 1478, 3283, 1633, 8536, 7019, 3947, 2789, 4606, 3545, 4421, 1289, 6416, 3647)
logRS <- log(Recruits / Spawners)
    ```

    a. Fit the following Ricker model to these data using the following linear form of this model with normally distributed errors:   
\begin{equation*}
log(R_t/S_t) = a + b \times S_t + e_t,\text{ where } e_t \sim \,\text{N}(0,\sigma^2)
\end{equation*}
    You will recognize that this form is exactly the same as linear regression, with independent errors (very similar to the intercept only model of Wind we fit in section \@ref(sec-jags-lr-no-covariates)).

    b. Within the constraints of the Ricker model, think about other ways you might want to treat the errors. The basic model described above has independent errors that are not correlated in time. Approaches to analyzing this dataset might involve

        *  modeling the errors as independent (as described above)
        
        *  modeling the errors as autocorrelated
        
        *  fitting a state-space model, with independent or correlated process errors

    Fit each of these models, and compare their performance (either using their predictive ability, or forecasting ability).


```{r jags-reset, include=FALSE, purl=FALSE}
file.remove(removefiles)
```
