```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
set.seed(123)
```

## Topics

Univariate state-space models

1. state (process) model

2. observation model


## State-space models

Consist of 2 parts


## Part 1: State model

Describes the __true state of nature__ over time

```{r, fig.height=4, fig.width=8, out.height="100%", out.width="100%", fig.align='center'}
par(mai=c(0.8,0.8,0,0), omi=rep(0,4))
## boundaries
ss <- 5
nn <- 7
rr <- ss*3
cc <- ss*nn
## mid-points
xm <- ss/2 + seq(0,cc-ss,ss)
ymt <- rr - ss/2
ymb <- ss/2
## arrow locs
x0t <- seq(ss, by=2*ss, len=3)
x1t <- x0t + ss
## empty plot space
plot(c(0,cc), c(0,rr), type="n", xlab="", ylab="",
     xaxt="n", yaxt="n", bty="n")
## top row: state
symbols(x=xm[c(1,3,5,7)], y=rep(ymt,4), circles=rep(ss/2,4),
        lty="solid",  fg=NA, bg="#488fdf",
        inches=FALSE, add=TRUE, lwd=3)
text("Truth", x=-ss, y=ymt, adj=c(0,0.5), xpd=NA,
     cex=2, col="#488fdf")
arrows(x0=x0t,x1=x1t,y0=ymt, col="#488fdf", lwd=3, length=0.12)
## Time or space
arrows(x0=ss/2, x1=cc-ss/2, y0=-ss/3+ss*2,
       length=0.12, lwd=3, xpd=NA)
text("Time", x=cc/2, y=-ss/2+ss*2, xpd=NA, pos=1, cex=2)
```


## Part 1: State model

States of nature might be:

 * animal's location
 
 * density of organisms
 
 * reproductive state
 

## Changes in the state of nature

We can think of changes in the state of nature being driven by a combination of

 * intrinsic (eg, fecundity), and

 * extrinsic factors (eg, temperature)


## Process errors

Some of the extrinsic drivers may be unknown

In time series modeling, we often call these unknown extrinsic factors _process errors_


## Part 2: Observation model

__Data = Truth + Error__
 
```{r obs_diag, fig.height=4, fig.width=8, out.height="100%", out.width="100%", fig.align='center'}
par(mai=c(0.8,0.8,0,0), omi=rep(0,4))
## arrow locs
x0t <- seq(ss, by=2*ss, len=3)
x1t <- x0t + ss
y0b <- rr - ss
y1b <- ss
## empty plot space
plot(c(0,cc), c(0,rr), type="n", xlab="", ylab="",
     xaxt="n", yaxt="n", bty="n")
## top row: state
symbols(x=xm[c(1,3,5,7)], y=rep(ymt,4), circles=rep(ss/2,4),
        lty="solid",  fg=NA, bg="#488fdf",
        inches=FALSE, add=TRUE, lwd=3)
text("Truth", x=-ss, y=ymt, adj=c(0,0.5), xpd=NA,
     cex=2, col="#488fdf")
## arrows
arrows(x0=x0t,x1=x1t,y0=ymt, col="#488fdf", lwd=3, length=0.12)
## bottom row: obs
symbols(x=xm[c(1,3,5,7)], y=rep(ss/2,4), circles=rep(ss/2,4),
        lty="solid",  fg=NA, bg="#844870",
        inches=FALSE, add=TRUE, lwd=3)
text("Data", x=-ss, y=ss/2, adj=c(0,0.5), xpd=NA,
     cex=2, col="#844870")
## arrows
arrows(x0=xm[c(1,3,5,7)], y0=y0b, y1=y1b,
       col="#c10101", lwd=3, length=0.12)
## Time or space
arrows(x0=ss/2, x1=cc-ss/2, y0=-ss/3,
       length=0.12, lwd=3, xpd=NA)
text("Time", x=cc/2, y=-ss/2, xpd=NA, pos=1, cex=2)
```


## Observation (sampling) errors

Consider for a moment the very act of collecting data

Information gathered depends on many factors

* Environmental conditions (eg, turbidity, currents)

* Behavior (eg, vertical migration, threat avoidance)

* Demographics (age, sex, maturity)

* Sampling design/coverage


## Consider these data

```{r}
set.seed(211)

TT <- 30

uu <- -0.2
  
ww <- xx <- rnorm(TT)
for(t in 2:TT) {
  xx[t] <- xx[t-1] + uu + ww[t]
}

ee <- rnorm(TT)
yy <- xx + ee

par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(yy, ylim = range(xx,yy),
        lwd = 2, type = "o", pch = 16, cex = 1.5, col = "blue",
        ylab = expression(italic(y[t])))
```


## Model options

__Linear regression__

State model

$$
x_t = \alpha + \beta t
$$

Observation model

$$
y_t = x_t + v_t ~ \text{with} ~ v_t \sim \text{N}(0,r)
$$

## Linear regression

```{r}
theta <- coef(lm(yy ~ seq(TT)))
y_obs <- theta[1] + theta[2] * seq(TT)

par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(yy, lwd = 2, type = "o", pch = 16,  cex = 1.5, col = "blue",
        ylab = expression(italic(x[t])~~or~~italic(y[t])))
lines(y_obs, lwd = 2, type = "o", pch = 16, col = "gray")
```

All of the variance in $y_t$ is due to observation error


## Observation errors

```{r}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(yy, lwd = 2, type = "o", pch = 16,  cex = 1.5, col = "blue",
        ylab = expression(italic(x[t])~~or~~italic(y[t])))
lines(y_obs, lwd = 2, type = "o", pch = 16, col = "gray")
segments(seq(TT), y0 = y_obs, y1 = yy, col = "red")
```

Errors $v_t = y_t - x_t$


## Model options

__Biased random walk__

State model

$$
x_t = x_{t-1} + u + w_t ~ \text{with} ~ w_t \sim \text{N}(0,q)
$$

Observation model

$$
y_t = x_t
$$


## Biased random walk

```{r}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(yy, lwd = 2, type = "o", pch = 16, cex = 1.5, col = "blue",
        ylab = expression(italic(y[t])))
lines(yy, lwd = 2, type = "o", pch = 16, col = "gray")
```

All of the variance in $y_t$ is due to process error


## Process errors

```{r}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(yy, lwd = 2, type = "o", pch = 16, cex = 1.5, col = "blue",
        ylab = expression(italic(y[t])))
lines(yy, lwd = 2, type = "o", pch = 16, col = "gray")
segments(seq(2,TT), y0 = yy[-TT], y1 = yy[-1], col = "red")
```

Errors: $w_t = y_t - y_{t-1}$


## Model options

__Biased random walk with observation error__

State model

$$
x_t = x_{t-1} + u + w_t ~ \text{with} ~ w_t \sim \text{N}(0,q)
$$

Observation model

$$
y_t = x_t + v_t ~ \text{with} ~ v_t \sim \text{N}(0,r)
$$


## Biased random walk

```{r}
xw <- rbind(cbind(seq(TT), xx), cbind(seq(TT)+0.5, xx+uu))
xw <- xw[order(xw[,1]),]
xw[,1] <- c(1, rep(seq(2, TT), ea = 2), TT+1)
  
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
## blank slate
plot.ts(yy, ylim = range(xx,yy), type = "n",
        ylab = expression(italic(x[t])))
lines(xw[,1], xw[,2], lwd = 2, col = "gray")
lines(xx, lwd = 2, type = "o", pch = 16, col = "black")
```


## Biased random walk with obs error

```{r}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
## blank slate
plot.ts(yy, ylim = range(xx,yy), type = "n",
        ylab = expression(italic(x[t])~~or~~italic(y[t])))
lines(xw[,1], xw[,2], lwd = 2, col = "gray")
lines(xx, lwd = 2, type = "o", pch = 16, col = "black")
segments(seq(TT), y0 = xx, y1 = yy, col = "red")
lines(yy, lwd = 2, type = "o", pch = 16, col = "blue", cex = 1.5)
```

The variance in $y_t$ is due to both process & observation error


## Biased random walk with obs error

```{r}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(yy, lwd = 2, type = "o", pch = 16, cex = 1.5, col = "blue",
        ylab = expression(italic(x[t])~~or~~italic(y[t])))
lines(xx, lwd = 2, type = "o", pch = 16, col = "gray")
```


## Separating process & obs errors

How is it possible to separate out both the process and observation errors?

They have different temporal patterns


## Separating process & obs errors

```{r, fig.align="center"}
xp <- matrix(NA, TT, 3)
xp[1,] <- c(0,0,0)

for(t in 2:TT) {
  xp[t,] <- xp[t-1,] + uu + rnorm(3)
}

yp <- xp[,1] + matrix(rnorm(TT*3), TT, 3)

par(mfcol = c(2,3), mai = c(0.6,0.6,0.1,0), omi = c(0,0,0.2,0))
for(i in 1:3) {
  plot.ts(xp[,i], xlab = "", ylab = expression(italic(x[t])))
  if(i == 2) {
    mtext(side = 3, " Different realizations of same process",
          line = 0.5, xpd = NA)
    }
  plot.ts(yp[,i], ylab = expression(italic(y[t])))
  if(i == 2) {
    mtext(side = 3, " Different observations of same process",
          line = 0.5, xpd = NA)
    }
}
```


## Density-dependent population growth

Stochastic Gompertz model

$$
N_t = N_{t-1} \text{e}^{(r_{max} + (b - 1) \log(N_{t-1}))} \text{e}^{w_t}  ~ \text{with} ~ w_t \sim \text{N}(0,q)
$$


## Density-dependent population growth

Stochastic Gompertz model

$$
N_t = N_{t-1} \text{e}^{(r_{max} + (b - 1) \log(N_{t-1}))} \text{e}^{w_t}  ~ \text{with} ~ w_t \sim \text{N}(0,r)
$$

Taking the log of both sides and substituting $x_t$ for $\log(N_t)$

$$
\begin{align}
\log(N_t) & = \log(N_{t-1}) + r_{max} + (b - 1) \log(N_{t-1}) + w_t \\
& \Downarrow \\
x_t &= x_{t-1} + r_{max} + (b - 1) x_{t-1} + w_t \\
    &= x_{t-1} + r_{max} + b x_{t-1} - x_{t-1} + w_t \\
    &= r_{max} + b x_{t-1} + w_t
\end{align}
$$

## Density-dependent population growth

What is the mean of $\{x_t\}$?

$$
x_t = r_{max} + b x_{t-1} + w_t
$$


## Density-dependent population growth

What is the mean of $\{x_t\}$?

$$
x_t = r_{max} + b x_{t-1} + w_t
$$

$$
\text{E}(x_t) = \frac{r_{max}}{1 - b}
$$


## Density-dependent population growth

$$
x_t = r_{max} + b x_{t-1} + w_t
$$

It turns out that $r_{max}$ and $b$ are confounded and create a likelihood surface with a strong ridge

This means we need a _long_ time series or mulitple observations of the process 

In practice, we will de-mean our data and instead use a familiar AR(1) model

$$
x_t = b x_{t-1} + w_t
$$


## Density-dependent population growth

If our population censuses contain observation or sampling errors, we can use a state-space version

$$
x_t = b x_{t-1} + w_t \\
y_t = x_t + v_t
$$


## Simulate a Gompertz model

Simulate the AR(1) state model

```{r, echo=TRUE}
## number of time steps
TT <- 50
## strength of density-dependence (0 < b < 1)
bb <- 0.5
## time series of process errors with SD = 1
ww <- rnorm(TT, 0, sqrt(1))
## initialize state & set x0 = w0
xx <- ww
## loop over time steps
for(t in 2:TT) {
  xx[t] <- bb * xx[t-1] + ww[t]
}
```


## Simulate a Gompertz model

Add some observation error

```{r, echo=TRUE}
## obs errors with var = 0.5
vv <- rnorm(TT, 0, sqrt(0.5))
## obs data
yy <- xx + vv
```


## Plot a random walk

Plot the state and observation

```{r, echo=TRUE, eval=FALSE}
plot.ts(xx, lwd = 2, type = "o", pch = 16, col = "gray",
        ylim = c(min(xx,yy), max(xx,yy)),
        ylab = expression(italic(x[t])~~or~~italic(y[t])))
lines(yy, lwd = 2, type = "o", pch = 16, col = "blue")
```


## Plot a random walk

```{r, align.fig="center"}
par(mai = c(0.8,0.8,0,0), omi = c(0,0,0,0))
plot.ts(xx, lwd = 2, type = "o", pch = 16, col = "gray",
        ylim = c(min(xx,yy), max(xx,yy)),
        ylab = expression(italic(x[t])~~or~~italic(y[t])))
lines(yy, lwd = 2, type = "o", pch = 16, col = "blue")
```


## Fitting models with `MARSS()`

We can use the __MARSS__ pkg for R to fit state-space models

To do so, we must write the model out just as we would on a white board


## Fitting models with `MARSS()`

Here are the main function arguments

```{r, echo=TRUE, eval=FALSE}
MARSS(y, model = NULL, inits = NULL, control = NULL, ...)
```

* `y` is an $n \times T$ `matrix` of data (observations)

* `model` is a `list` that defines the state-space model

* `inits` is an optional `list` of initial values

* `control` is an optional `list` of control parameters


## Fitting models with `MARSS()`

`MARSS()` uses a slightly expanded form of state-space model, such that

$$
\begin{align}
x_t &= b x_{t-1} + w_t \\
y_t &= x_t + v_t \\
& \Downarrow \\
x_t &= b x_{t-1} + u + w_t \\
y_t &= Z x_t + a + v_t
\end{align}
$$

with $u = 0$, $Z = 1$, and $a = 0$


## Fitting models with `MARSS()`

$$
x_t = b x_{t-1} + u + w_t ~ \text{with} ~ w_t \sim \text{N}(0,q)  \\
y_t = Z x_t + a + v_t ~ \text{with} ~ v_t \sim \text{N}(0,r)
$$

`MARSS()` works on matrices, so we have to define any scalar as a $1 \times 1$ matrix

```{r, echo=TRUE}
mod_list <- list(
  ## state model
  B = matrix("b"), U = matrix(0), Q = matrix("q"),
  ## obs model
  Z = matrix(1), A = matrix(0), R = matrix("r")
  )
```


## Fitting models with `MARSS()`

What about the initial state $x_0$?

$$
x_0 \sim \text{N}(\pi,\lambda)
$$

Estimating both $\pi$ and $\lambda$ is nearly impossible, so the default is to treat $x_0$ as fixed, but unknown, effect rather than a random effect, such that

$$
x_0 \sim \text{N}(\pi,0)
$$

`MARSS()` will do this for us


## Fitting models with `MARSS()`

We can ignore the `inits` and `control` arguments for now and fit the model

```{r, echo=TRUE, eval=FALSE}
## load MARSS package
library(MARSS)
## define the data as an N (rows) x T (cols) matrix
dat <- matrix(yy, nrow = 1, ncol = TT)
## fit the model
mod_fit <- MARSS(y = dat, model = mod_list)
```


## Fitting models with `MARSS()`

```{r, echo=FALSE, eval=TRUE}
## load MARSS package
library(MARSS)
## define the data as an N (rows) x T (cols) matrix
dat <- matrix(yy, nrow = 1, ncol = TT)
## fit the model
mod_fit <- MARSS(y = dat, model = mod_list)
```


## Covariates in state-space models

We can include covariates (explanatory variables) in our models as well

Covariates in the state model _affect the underlying process_

$$
x_t = b x_{t-1} + C c_t + w_t
$$


## Covariates in state-space models

We can include covariates (explanatory variables) in our models as well

Covariates in the state model _affect the underlying process_

$$
x_t = b x_{t-1} + C c_t + w_t
$$

Covariates in the observation model are _offsets to the underlying process_

$$
y_t = Z x_t + D d_t + v_t
$$


## Covariates in state-space models

```{r}
ww <- xx <- xy <- rnorm(TT)

bb <- 0.7

CC <- DD <- 2
cc <- dd <- sin(2*pi*seq(TT)/12)
  
for(t in 2:TT) {
  xx[t] <- bb * xx[t-1] + CC * cc[t] + ww[t]
  xy[t] <- bb * xy[t-1] + ww[t]
}

ee <- rnorm(TT)
yy <- xy + DD * dd # + ee

par(mai = c(0.9,0.9,0.5,0.1), omi = c(0,0,0,0))
plot.ts(xx, ylim = range(xx,yy),
        lwd = 2, type = "o", pch = 16, col = "gray",
        ylab = expression(italic(x[t])~~or~~italic(y[t])))
mtext(side = 3, expression(italic(x[t])==italic(b)~italic(x[t-1])~+~italic(C)~italic(c[t])~+~italic(w[t])),
      line = 0.5, adj = 0)
```


## Covariates in state-space models

```{r}
par(mai = c(0.9,0.9,0.5,0.1), omi = c(0,0,0,0))
plot.ts(xx, ylim = range(xx,yy),
        lwd = 2, type = "o", pch = 16, col = "gray",
        ylab = expression(italic(x[t])~~or~~italic(y[t])))
lines(yy, lwd = 2, type = "o", pch = 16, cex = 1.5, col = "blue")
mtext(side = 3, expression(italic(x[t])==italic(b)~italic(x[t-1])~+~italic(w[t])),
      line = 1.5, adj = 0)
mtext(side = 3, expression(italic(y[t])==italic(x[t])~+~italic(D)~italic(d[t])),
      line = 0.5, adj = 0)
```



