```{r mlr-setup, include=FALSE, purl=FALSE}
knitr::opts_knit$set(unnamed.chunk.label = "usssolns-")
knitr::opts_chunk$set(fig.align='center', out.width = '90%', cache=TRUE, warning=FALSE)
```

# Solutions Chapter 6
\chaptermark{Solutions Ch 6}

Here are the answers for the homework problems based on the material in [Univariate state-space models](https://nwfsc-timeseries.github.io/atsa-labs/chap-univariate-state-space.html).  You can download the Rmd file for this key from [here](HW_3_uss_Key.Rmd).

## Data Set Up

```{r}
library(MARSS)
data(grouse, package="MARSS")
dat <- log(grouse[,2])
```

## Problem 1

**Write the equations for each of these models: ARIMA(0,0,0), ARIMA(0,1,0), ARIMA(1,0,0), ARIMA(0,0,1), ARIMA(1,0,1).**

In all cases, $e_t \sim N(0,\sigma^2)$ and $e_t$ and $e_{t-1}$ are independent.  Thus the $e_t$ are i.i.d. (independent and identically distributed).

ARIMA(0,0,0) = white noise
$$x_t = e_t$$

ARIMA(0,1,0) = differenced data is white noise
$$x_t - x_{t-1} = e_t$$

ARIMA(1,0,0) = Autoregressive lag-1.  This is a mean-reverting random walk if $|b|<1$. If $b=1$, it is a simple random walk and the same as ARIMA(0,1,0).
$$x_t = bx_{t-1} + e_t$$

ARIMA(0,0,1) is a moving-average lag-1 model.
$$x_t = e_t + \theta e_{t-1}$$

ARIMA(1,0,1) is a moving-average lag-1 model with autoregression lag-1.
$$x_t = b x_{t-1} + e_t + \theta e_{t-1}$$


## Problem 2

a. **Plot the data.**

    ```{r}
    plot(grouse[,1], grouse[,2], type="l", ylab="log count", xlab="")
    ```

b. **Fit each model using `MARSS()`.**

    ```{r mod-lists, results='hide'}
    mod.list1 <- list(
      B=matrix(1), U=matrix(0), Q=matrix("q"),
      Z=matrix(1), A=matrix(0), R=matrix(0),
      x0=matrix("a"), tinitx=0)
    fit1.marss <- MARSS(dat, model=mod.list1)

    mod.list2 <- list(
      B=matrix(1), U=matrix("u"), Q=matrix("q"),
      Z=matrix(1), A=matrix(0), R=matrix(0),
      x0=matrix("a"), tinitx=0)
    fit2.marss <- MARSS(dat, model=mod.list2)
    ```

    Parameter estimates:
    ```{r}
    coef(fit1.marss, type="vector")
    coef(fit2.marss, type="vector")
    ```

c. **Which one appears better supported given AICc?**

    ```{r}
    c(fit1.marss$AICc, fit2.marss$AICc)
    ```
    
    Model 2 is better supported with a $\Delta AICc$ of `r fit2.marss$AICc-fit1.marss$AICc`.

c. **Load the forecast package. Use `auto.arima(dat)`  to fit the data.  Next run `auto.arima()` on the data with `trace=TRUE` to see all the ARIMA models it compared.**

    ```{r}
    library(forecast)
    auto.arima(dat)
    ```

    Let's look at all the models it tried using `trace=TRUE`.
    ```{r}
    auto.arima(dat, trace=TRUE)
    ```
    
    It picked model 2 as the best among those tested.  "ARIMA(0,1,0) with drift" is model 2.  

e. **Is the difference in the AICc values between a random walk with and without drift comparable between `MARSS()` and `auto.arima()`?**

    ```{r}
    fit1.arima <- Arima(dat, order=c(0,1,0))
    fit2.arima <- Arima(dat, order=c(0,1,0), include.drift=TRUE)
   
    fit2.arima$aicc-fit1.arima$aicc
    fit2.marss$AICc-fit1.marss$AICc
    ```
   
    Similar but not identical.  

## Problem 3

This code produces $x_t = x_{t-1} + u + w_t$ data with $u=0.1$ and $q=1$.
```{r}
dat <- cumsum(rnorm(100,0.1,1))
```

a. **Write out the equation for that random walk as a univariate state-space model.**

    \begin{equation}
    x_t = x_{t-1} + u + w_t, w_t \sim \text{N}(0,q) \\
    x_0 = \mu \text{ or } x_1 = y_1 \\
    y_t = x_t
    \end{equation}
    where $u=0.1$ and $q=1$.

b. **What is the order of the $\mathbf{x}$ part of the model written as ARIMA(p, d, q)?** 

    From Problem 1, you should be able to deduce it is ARIMA(0,1,0) but if you said ARIMA(1,0,0) with b=1, that's ok.  That's not how `Arima()` writes $x_t = x_{t-1} + u + w_t$ but it is correct.

c. **Fit that model using `Arima()` in the forecast package.  You'll need to specify the `order` and `include.drift` term.**

    ```{r hw2-arima, results='hide'}
    fit.arima <- Arima(dat, order=c(0,1,0), include.drift=TRUE)
    ```
    
d. **Fit that model with `MARSS()`.**

    ```{r hw2-marss, results='hide'}
    mod.list <- list(
      B=matrix(1), U=matrix("u"), Q=matrix("q"),
      Z=matrix(1), A=matrix(0), R=matrix(0),
      x0=matrix("mu"), tinitx=0)
    fit.marss <- MARSS(dat, model=mod.list)
    ```
    
    or since you know that $x_1 = y_1$ from the observation model, you could use:

    ```{r hw2-marss-alt, results='hide'}
    mod.list.alt <- list(
      B=matrix(1), U=matrix("u"), Q=matrix("q"),
      Z=matrix(1), A=matrix(0), R=matrix(0),
      x0=matrix(dat[1]), tinitx=1)
    fit.alt.marss <- MARSS(dat, model=mod.list.alt, method="BFGS")
    ```

    But hopefully you didn't try that because this makes the likelihood surface flat and you need to use many more iterations or try a Newton method which happens to help (sometimes it doesn't or makes thing worse).

e. **How are the two estimates different?**

    ```{r hw2-compare-coef}
    coef(fit.marss, type="vector")
    coef(fit.alt.marss, type="vector")
    c(coef(fit.arima), s2=fit.arima$sigma2)
    ```
    
    `MARSS()` is estimating 3 parameters while `Arima()` is estimating 2.  The $u$ estimates are identical (or very similar) but the $q$ estimate is different.

    `coef()` is the standard function for getting estimates from fits.  Try `?coef` to find the help file for `coef` applied to MARSS objects.  For Arima objects, `coef()` doesn't return sigma2.  It is in `fit.arima$sigma2`.

    Now fit the first-differenced data:
    ```{r hw3.data}
    diff.dat <- diff(dat)
    ```

f. **If $x_t$ denotes a time series.  What is the first difference of $x$?  What is the second difference?**

    First difference `diff(x)` is $x_t - x_{t-1}$.

    Second difference is `diff(diff(x))` or $(x_t - x_{t-1}) - (x_{t-1} - x_{t-2})$.

g. **What is the $\mathbf{x}$ model for `diff.dat`?**

$$\text{diff}(x)=(x_t - x_{t-1}) = u + w_t$$

h. **Fit `diff.dat` using `Arima()`. You'll need to change `order` and `include.mean`.** _Note: I should have had you use `arima()` in the question; `Arima()` is reporting the unbiased variance estimate while `MARSS()` and `arima()` are reporting the straight maximum-likelihood estimate.  The difference is $(n-1)/n$ where $n$ is the length of the data, in this case the difference data.  So $(n-1)/n=98/99$._

    ```{r hw2-fit-diff-dat}
    fit.diff.Arima <- Arima(diff.dat, order=c(0,0,0), include.mean=TRUE)
    fit.diff.arima <- arima(diff.dat, order=c(0,0,0), include.mean=TRUE)
    ```
    
i. **Fit the model in (h) with `MARSS()`.**

    `data` ($y$) is now `diff.dat` and state-space model is below. It doesn't matter what $x_0$ is as it does not appear in the model, but it is important to use $x_0$ instead of $x_1$ to match `arima()`.

    \begin{equation}
    x_t = u + w_t, w_t \sim \mathbf{N}(0,q) \\
    x_0 = 0 \\
    y_t = x_t
    \end{equation}

    ```{r hw2-fit-diff-marss, results='hide'}
    mod.list.diff.1 <- list(
      B=matrix(0), U=matrix("u"), Q=matrix("q"),
      Z=matrix(1), A=matrix(0), R=matrix(0),
      x0=matrix(0), tinitx=0)
    fit.diff.marss.1 <- MARSS(diff.dat, model=mod.list.diff.1)
    ```
    
    Or we could have written it like so
    \begin{equation}
    x_t = 0, w_t \sim \mathbf{N}(0,0) \\
    x_0 = 0 \\
    y_t = a + v_t, v_t \sim \mathbf{N}(0,r)
    \end{equation}
    
    In this case, we fit the model as:
    
    ```{r hw2-fit-diff-marss-2, results='hide'}
    mod.list.diff.2 <- list(
      B=matrix(0), U=matrix(0), Q=matrix(0),
      Z=matrix(0), A=matrix("u"), R=matrix("r"),
      x0=matrix(0), tinitx=0)
    fit.diff.marss.2 <- MARSS(diff.dat, model=mod.list.diff.2)
    ```

    Note, we can also fit with lm():
    ```{r hw2-fit-diff-marss-lm, results='hide'}
    fit.diff.lm <- lm(diff.dat~1)
    ```

    Here are the parameter estimates. 
    ```{r hw2-fit-diff-coef}
    rbind(
      marss.diff.1 <- coef(fit.diff.marss.1, type="vector"),
      marss.diff.2 <- coef(fit.diff.marss.2, type="vector"),
      arima.diff <- c(coef(fit.diff.arima), s2=fit.diff.arima$sigma2),
      Arima.diff <- c(coef(fit.diff.Arima), s2=(98/99)*fit.diff.Arima$sigma2),
      lm.diff <- c(coef(fit.diff.lm), s2=(98/99)*summary(fit.diff.lm)$sigma^2)
      )
    ```
    
    They are all the same except the variances reported by `Arima()` and `lm()` have to be multiplied by $98/99$ to be the same as `MARSS()` and `arima()` because the former are reporting the unbiased estimates and the latter are reporting the straight (biased) maximum-likelihood estimates.

## Problem 4

\begin{equation}
x_t = b x_{t-1}+u+w_t \text{ where } w_t \sim \mathbf{N}(0,q)  
\label{eq:gompertz}\end{equation}

a. **Write R code to simulate the equation above.  Make $b$ less than 1 and greater than 0.  Set $u$ and $x_0$ to whatever you want.  You can use a `for` loop.**

    ```{r hw3-sim}
    #set up my parameter values
    b <- .8; u <- 2; x0 <- 10; q <- 0.1
    nsim <- 1000
    #set up my holder for x
    x <- rep(NA, nsim)
    x[1] <- b*x0+u+rnorm(1,0,sqrt(q))
    for(t in 2:nsim) x[t] <- b*x[t-1]+u+rnorm(1,0,sqrt(q))
    ```

b. **Plot the trajectories and show that this model does not ``drift'' upward or downward.  It fluctuates about a mean value.**

    ```{r hw3-fig-plot}
    par(mar=c(2, 4, 2, 2))
    plot(x, type="l",xlab="", ylab="x", ylim=c(-6+u/(1-b),6+u/(1-b)))
    ```

c. **Hold $b$ constant and change $u$.  How do the trajectories change?**

    ```{r}
    #set up my parameter values
    u2 <- u+1
    x2 <- rep(NA, nsim)
    x2[1] <- b*x0+u2+rnorm(1,0,sqrt(q))
    for(t in 2:nsim) x2[t] <- b*x2[t-1]+u2+rnorm(1,0,sqrt(q))
    #second u
    u3 <- u-1
    x3 <- rep(NA, nsim)
    x3[1] <- b*x0+u3+rnorm(1,0,sqrt(q))
    for(t in 2:nsim) x3[t] <- b*x3[t-1]+u3+rnorm(1,0,sqrt(q))
    ```
    
    Plot the simulations:
    
    ```{r hw3-fig-plot2}
    par(mar=c(2, 4, 2, 2))
    plot(x, type="l",xlab="", ylab="x", ylim=c(-6+u/(1-b),6+u/(1-b)))
    lines(x2, col="blue")
    lines(x3, col="red")
    legend("bottomright", c("u+1","u-1"), col=c("blue", "red"), lty=1, bg="white")
    ```
    
    $u$ moves the mean of the trajectories up or down.


d. **Hold $u$ constant and change $b$.  Make sure to use a $b$ close to 1 and another close to 0. How do the trajectories change?**

    ```{r hw3-sim2}
    #set up my parameter values
    b1 <- 0.9
    x0 <- u/(1-b1)
    x1 <- rep(NA, nsim)
    x1[1] <- b1*x0+u+rnorm(1,0,sqrt(q))
    for(t in 2:nsim) x1[t] <- b1*x1[t-1]+u+rnorm(1,0,sqrt(q))
    # second b
    b2 <- 0.1
    x0 <- u/(1-b2)
    x2 <- rep(NA, nsim)
    x2[1] <- b2*x0+u+rnorm(1,0,sqrt(q))
    for(t in 2:nsim) x2[t] <- b2*x2[t-1]+u+rnorm(1,0,sqrt(q))
    ```
    
    Plot the simulations:
    
    ```{r hw3-fig-plot3}
    par(mfrow=c(1,2), mar=c(2,2,2,1))
    plot(x1, type="l",xlab="", ylab="x", main="b=0.9")
    plot(x2, type="l",xlab="", ylab="x", main="b=0.1")
    ```

    The one with smaller $b$ has less auto-regression and is `tighter' (explores less of a range of the y axis).

e. **Do 2 simulations each with the same $w_t$.  In one simulation, set $u=1$ and in the other $u=2$.  For both simulations, set $x_1 = u/(1-b)$.  You can set $b$ to whatever you want as long as $0<b<1$.  Plot the 2 trajectories on the same plot.  What is different?**

    ```{r hw3-sim3}
    #set up my parameter values
    b <- 0.9
    u <- 1
    x0 <- u/(1-b)
    err <- rnorm(nsim,0,sqrt(q))
    x1 <- rep(NA, nsim)
    x1[1] <- b*x0+u+err[1]
    for(t in 2:nsim) x1[t] <- b*x1[t-1]+u+err[t]
    # second u
    u <- 2
    x0 <- u/(1-b)
    x2 <- rep(NA, nsim)
    x2[1] <- b*x0+u+err[1]
    for(t in 2:nsim) x2[t] <- b*x2[t-1]+u+err[t]
    ```
    
    Plot the simulations:
    ```{r hw3-fig-plot4}
    par(mfrow=c(1,2), mar=c(2,2,2,1))
    plot(x1, type="l",xlab="", ylab="x", main="u=1")
    plot(x2, type="l",xlab="", ylab="x", main="u=2")
    ```
    
    They are exactly the same except that the mean has changed from $1/(1-b)$ to $2/(1-b)$.  The mean level in the AR-1 model $x_t = b x_{t-1} + u + w_t$ is $u/(1-b)$.  For a given $b$, $u$ just changes the level.


## Problem 5

The MARSS package includes a data set of gray whales.  Load the data to use as follows:
```{r}
library(MARSS)
data(graywhales, package="MARSS")
dat <- log(graywhales[,2])
```

Fit a random walk with drift model observed with error to the data:

\begin{equation}
x_t = x_{t-1}+u+w_t \text{ where } w_t \sim \mathbf{N}(0,q) \\
y_t = x_t+v_t \text{ where } v_t \sim \mathbf{N}(0,r) \\
x_0 = a 
\end{equation}

$y$ is the whale count in year $t$. $x$ is interpreted as the 'true' unknown population size that we are trying to estimate.

a. **Fit this model with `MARSS()`**

    ```{r hw4-fit}
    mod.list=list(
      B=matrix(1), U=matrix("u"), Q=matrix("q"),
      Z=matrix(1), A=matrix(0), R=matrix("r"),
      x0=matrix("mu"), tinitx=0)
    fit.marss = MARSS(dat, model=mod.list)
    ```

b. **Plot the estimated $x$ as a line with the actual counts added as points.**

    ```{r hw4-fig-plot1}
    par(mar=c(2,2,2,2))
    plot(graywhales[,1], fit.marss$states[1,], type="l",xlab="", ylab="log count")
    points(graywhales[,1], dat)
    ```

c. **Simulate 1000 sample trajectories using the estimated $u$ and $q$ starting at the estimated $x$ in 1997.  You can do this with a couple `for` loops or write something terse with `cumsum()` and `apply()`.**

    ```{r hw4-sim}
    #1997 is the 39th (last) data point
    x0 <- fit.marss$states[1,39]
    q <- coef(fit.marss)$Q
    u <- coef(fit.marss)$U
    #next question asks for pop size in 2007 so nforeward=10
    nsim <- 1000
    nforeward <- 10
    #each row holds a simulation
    x <- matrix(NA, nsim, nforeward)
    x[,1] <- x0+u+rnorm(nsim,0,sqrt(q))
    for(t in 2:nforeward) x[,t] <- x[,t-1]+u+rnorm(nsim,0,sqrt(q))
    ```
    
d. **Using this what is your estimated probability of reaching 50,000 graywhales in 2007.**

    The question was phrased a big vaguely.  It does not specify if this means ``in 2007, x>=log(50000)'', ``at some point by or before 2007, x reaches log(50000) at least once'', or ``in 2007, the population is at least 50000 whales''.  As long as you stated what you were trying to estimate, you were fine.  Let's use a simulation to compute the last option.
    
    ```{r hw4-prob}
    xthresh <- log(50000)
    sum(x[,10]>=xthresh)/nsim
    ```
    
e. **What kind of uncertainty does that estimate NOT include?**

    By using the point estimates of $u$, $q$ and $x_0$, we are not including the uncertainty in those estimates in our forecasts.


## Problem 6

**Fit the following 3 models to the graywhales data using `MARSS()`:**

1. Process error only model with drift
2. Process error only model without drift
3. Process error with drift and observation error with observation error variance fixed = 0.05. 
4. Process error with drift and observation error with observation error variance estimated. 

Process error only with drift. $x_t = x_{t-1} + u + w_t$ with $y_t = x_t$.

```{r hw5-fit1, results='hide'}
mod.list <- list(
  B=matrix(1), U=matrix("u"), Q=matrix("q"),
  Z=matrix(1), A=matrix(0), R=matrix(0),
  x0=matrix("mu"), tinitx=0)
fit.whales1 = MARSS(dat, model=mod.list)
```

Process error only without drift. $x_t = x_{t-1} + w_t$ with $y_t = x_t$.
```{r hw5-fit2, results='hide'}
mod.list <- list(
  B=matrix(1), U=matrix(0), Q=matrix("q"),
  Z=matrix(1), A=matrix(0), R=matrix(0),
  x0=matrix("mu"), tinitx=0)
fit.whales2 = MARSS(dat, model=mod.list)
```

Process error only with drift. $x_t = x_{t-1} + w_t$ with $y_t = x_t+v_t, v_t \sim N(0,0.05)$.
```{r hw5-fit3, results='hide'}
mod.list <- list(
  B=matrix(1), U=matrix("u"), Q=matrix("q"),
  Z=matrix(1), A=matrix(0), R=matrix(0.05),
  x0=matrix("mu"), tinitx=0)
fit.whales3 = MARSS(dat, model=mod.list)
```

Process error only with drift. $x_t = x_{t-1} + w_t$ with $y_t = x_t+v_t, v_t \sim N(0,r)$.
```{r hw5-fit4, results='hide'}
mod.list <- list(
  B=matrix(1), U=matrix("u"), Q=matrix("q"),
  Z=matrix(1), A=matrix(0), R=matrix("r"),
  x0=matrix("mu"), tinitx=0)
fit.whales4 = MARSS(dat, model=mod.list)
```

**Compute the AICc's for each model and likelihood or deviance (-2 * log likelihood)**

```{r hw5-aicc}
c(fit.whales1$AICc, fit.whales2$AICc, 
  fit.whales3$AICc, fit.whales4$AICc)
c(fit.whales1$logLik, fit.whales2$logLik, 
  fit.whales3$logLik, fit.whales4$logLik)
```

**Calculate a table of delta-AICc values and AICc weights.**

```{r hw5-aicc-table}
AICc <- c(fit.whales1$AICc, fit.whales2$AICc, 
       fit.whales3$AICc, fit.whales4$AICc)
delAIC <- AICc-min(AICc)
relLik <- exp(-0.5*delAIC)
aic.table <- data.frame(
  AICc = AICc,
  delAICc = delAIC,
  relLik = relLik,
  weight = relLik/sum(relLik)
)
rownames(aic.table) <- c(
  "proc only with drift", 
  "proc only no drift", 
  "proc with drift and obs error fixed",
  "proc with drift and obs error estimated")
knitr::kable(round(aic.table, digits=3))
```

There is not much data support for including observation error with $r=0.05$.  But that is because $r=0.05$ is too big.  If we estimate $r$, the process error with drift and observation error model would is best.

**Show the acf of the model and state residuals for the best model.  Do they suggest any problems?**

We need to specify `na.action` since the state residuals have an NA at the last value.

```{r}
best.fit <- fit.whales4
sres <- residuals(best.fit)$state.residuals[1,]
mres <- residuals(best.fit)$model.residuals[1,]
par(mfrow=c(1,2), mar=c(1,2,2,1))
acf(sres, na.action=na.omit)
acf(mres, na.action=na.omit)
```

The state residuals have some autocorrelation suggesting that the model with uncorrelated errors does not fully capture the temporal correlation in the data.  The model residuals look 'ok' but on plotting the residuals, we see something odd.  

```{r}
par(mfrow=c(1,2), mar=c(1,2,2,1))
plot(sres, main="state residuals")
plot(mres, main="model residuals")
```

The model residuals have a 0 when the data are NA.  In the MARSS output, model residual = expected value of (y|all data) - fitted value of (y|all data).  Fitted is $x_{t|y_{1:T}}+a$.  When the data are missing, expected value of y is fitted value, so model residual is 0.  However, when evaluating model residuals, we should set the model residual to 0 when the data are missing.  The state residuals are also 0 when data are missing.  We might want to set those to NA also so they do not skew our auto-correlation estimates.  Since this creates internal NAs, we need to set `na.action=na.pass`.  Now the state residuals look better.

```{r}
par(mfrow=c(1,2), mar=c(1,2,2,1))
sres[is.na(dat)]=NA
mres[is.na(dat)]=NA
acf(sres, na.action=na.pass)
acf(mres, na.action=na.pass)
```

We can plot the residuals nicely with `plot()` also.

```{r}
plot(best.fit, plot.type="state.residuals")
plot(best.fit, plot.type="model.residuals")
```

## Problem 7

**Load the data to use as follows and set up so you can use the last 3 data points to validate your fits.**

```{r hw5-data}
library(forecast)
dat <- log(airmiles)
n <- length(dat)
training.dat <- dat[1:(n-3)]
test.dat <- dat[(n-2):n]
```

a. **Fit the following four models using `Arima()`: ARIMA(0,0,0), ARIMA(1,0,0), ARIMA(0,0,1), ARIMA(1,0,1).**

    ```{r hw6-fit}
    fit.1 <- Arima(training.dat, order =c(0,0,0))
    fit.2 <- Arima(training.dat, order =c(1,0,0))
    fit.3 <- Arima(training.dat, order =c(0,0,1))
    fit.4 <- Arima(training.dat, order =c(1,0,1))
    ```

b. **Use `forecast()` to make 3 step ahead forecasts from each.**

    ```{r hw6-forecast}
    forecast.1 <- forecast(fit.1, h=3)
    forecast.2 <- forecast(fit.2, h=3)
    forecast.3 <- forecast(fit.3, h=3)
    forecast.4 <- forecast(fit.4, h=3)
    ```

c. **Calculate the MASE statistic for each using the `accuracy()` function in the forecast package.** 
 
    ```{r hw-accuracy-1}
    accuracy(forecast.1, test.dat)
    ```
    
    The MASE statistic we want is in the Test set row and MASE column. 
    
    ```{r hw6-mase}
    MASEs <- c(
      accuracy(forecast.1, test.dat)["Test set","MASE"],
      accuracy(forecast.2, test.dat)["Test set","MASE"],
      accuracy(forecast.3, test.dat)["Test set","MASE"],
      accuracy(forecast.4, test.dat)["Test set","MASE"]
      )
    ```

d. **Present the results in a table.**

    ```{r hw6-mase-table}
    data.frame(
      name=paste("Arima",c("(0,0,0)","(1,0,0)","(0,0,1)","(1,0,1)"),sep=""),
      MASE=MASEs
      )
    ```

e. **Which model is best supported based on the MASE statistic?**

    What this table shows is that the ARMA(1,0,1) is the best, and the AR component strongly improves predictions

## Problem 8

**Set up the data**
```{r hw8-data}
turtlename <- "MaryLee"
dat <- loggerheadNoisy[which(loggerheadNoisy$turtle==turtlename),5:6]
dat <- t(dat) 
```

a. **Plot MaryLee's locations (as a line not dots).  Put the latitude locations on the y-axis and the longitude on the y-axis.**

    ```{r hw8-plot-code}
    plot(dat[1,],dat[2,], type="l")
    ```

b. **Analyze the data with a state-space model (movement observed with error) using**

    ```{r hw-movement-fit1, results='hide'}
    fit0 <- MARSS(dat)
    ```

    $U_{lon}$ is the average velocity in N-S direction. $U_{lat}$ is the average velocity in E-W direction. $R_{diag}$ is the observation error variance. $Q$'s are the movement error variances. $x_0$'s are the estimated positions (lat/lon) at $t=0$.

c. **What assumption did the default MARSS model make about observation error and process error?**

    The observation errors in the lat and lon direction are independent but have identical variance.  The movement errors are independent (not correlated) and allowed to have different variances.  So the model doesn't allow a average NE movement; that would require correlation in the movement errors.  It allows that turtles tend to move faster N-S (along the coast) than E-W (out to sea).

d. **Does MaryLee move faster in the latitude direction versus longitude direction?**

    No. The estimated $u$'s in the lat and lon direction are similar.

e. **Add MaryLee's estimated "true" positions to your plot of her locations.** 

    ```{r hw8-plot2-code}
    plot(dat[1,],dat[2,], type="l")
    lines(fit0$states[1,], fit0$states[2,], col="red")
    ```

f. **Compare the following models for these data.  Movement in the lat/lon direction is  (1) independent but the variance is the same, (2) is correlated and lat/lon variances are different, and (3) is correlated and the lat/lon variances are the same.**

    ```{r hw8-movement-model-compare, results='hide'}
    fit1 <- MARSS(dat, model=list(Q="diagonal and equal"))
    fit2 <- MARSS(dat, model=list(Q="unconstrained"))
    fit3 <- MARSS(dat, model=list(Q="equalvarcov"))
    ```

    ```{r hw8-movement-model-compare-aics}
    c(fit0=fit0$AICc, fit1=fit1$AICc, 
      fit2=fit2$AICc, fit3=fit3$AICc)
    ```

    The model with correlated movement but equal movement error variances is best supported.  This suggests a tendency to move in a particular direction (probably up down the coast).  However, actually this is caused by strong directional movement in the middle of the movement track.

g. **Plot your state residuals (true location residuals).  What are the problems? Discuss in reference to your plot of the location data.**  

    Here is how to get state residuals from MARSS().

    ```{r hw8-resids-plot-code}
    par(mfrow=c(2,2),mar=c(3,5,3,5))
    resids.lon <- residuals(fit3)$state.residuals[1,]
    plot(resids.lon, xlab=""); abline(h=0)
    acf(resids.lon, na.action=na.pass)
    resids.lat <- residuals(fit3)$state.residuals[2,]
    plot(resids.lat, xlab=""); abline(h=0)
    acf(resids.lat, na.action=na.pass)
    ```

    There is a period in the middle of the track where the model does not describe the movement well.  We can see in the plot that the turtle has a long northward movement in the middle of the track.
