910-Appendices.Rmd

# Appendices

## Random Walk Model with dlm {#randomWalkWithDLM}

The [Random Walk](#randomWalk) chapter shows how to
create a random walk model using the `StructTS` function.
You can also create a random walk model by using the dlm package,
and this appendix shows how.

### Estimation
The code for estimating parameters is very similar to
the code for estimating the local level model.
The difference is that we force $V$, the variance of the observations,
to be zero.

We define a function, `buildRandomWalk`, that builds a `dlm` model object
from two input parameters, `dW` and `m0`.
The parameters are packed into a single, 2-element vector.

```{r, eval=FALSE}
buildRandomWalk <- function(v) {
  dW <- exp(v[1])
  m0 <- v[2]
  dlmModPoly(order=1, dV=0, dW=dW, m0=m0)
}
```

The function calls the `dlmModPoly` function from `dlm` to create the model object.

We need initial guesses for the model parameters.

```{r, eval=FALSE}
varGuess <- var(diff(y), na.rm=TRUE)
mu0Guess <- as.numeric(y[1])
```

Next we call the `dlmMLE` function to estimate the MLE parameters
using numerical optimzation.
Always check for convergence.

```{r, eval=FALSE}
parm <- c(log(varGuess), mu0Guess)
mle <- dlmMLE(y, parm=parm, buildRandomWalk)
if (mle$convergence != 0) stop(mle$message)
```

From the MLE estimates, we can build the final `dlm` model.

```{r, eval=FALSE}
model <- buildRandomWalk(mle$par)
```

We can extract the estimated parameters from `model`, the returned object.

---------- -----------------------------------
`model$W`  Variance of the random walk errors
`model$m0` Initial level
---------- -----------------------------------

### Example {-}

```{r, eval=TRUE}
library(dlm)

y <- datasets::Nile

buildRandomWalk <- function(v) {
  dW <- exp(v[1])
  m0 <- v[2]
  dlmModPoly(order=1, dV=0, dW=dW, m0=m0)
}

varGuess <- var(diff(y), na.rm=TRUE)
mu0Guess <- as.numeric(y[1])

parm <- c(log(varGuess), mu0Guess)
mle <- dlmMLE(y, parm=parm, buildRandomWalk)
if (mle$convergence != 0) stop("Optimizer did not converge")

model <- buildRandomWalk(mle$par)

cat("Transitional variance:", model$W, "\n",
    "Initial level:", model$m0, "\n")
```

## Local Level Model with dlm {#localLevelModelWithDLM}

The chapter on the [local level model](#localLevelModel)
uses the `StructTS` function because that's the easiest way to estimate the model parameters.
Sometimes, however, you might want to use the `dlm` package instead,
even though it's a bit more work.
Why would one do that?
The local level model might be your first step in model building,
leading to more complicate models.
Or you might want to bootstrap your model, which is more easily done using `dlm`.
Or you might want to combine a local level model with another model
using the model "addition" feature of `dlm`.

The `dlm` authors refer to the local level model as
the _random walk with noise_ model:
the underlying level follows a random walk,
and our observation of it is polluted by noise.

Mathematically, the local level models used by the `StructTS` function and the `dlm` package
are the same,
but they use different variable names
and slightly different notational conventions.

\begin{eqnarray*}
  Y_{t} & = & \mu_{t} + v_{t}, \qquad v_{t} \sim N(0, V) \\
  \mu_{t} & = & \mu_{t-1} + w_{t}, \qquad w_{t} \sim N(0, W)
\end{eqnarray*}

Under these conventions,
we observe $Y_{t}$ (not $y_{t}$),
and the variances of the error terms are generalized to be matrices $V$ and $W$.

(*Move to footnote:* Generalizing $V$ and $W$ to matrices will open the door to the multivariate case.)

Following those conventions, the model has these three parameters.

------ ------------------------------------
 `dV`  Variance of the observation errors
 `dW`  Variance of the transition errors
 `m0`  The initial value ($\mu_{0}$)
------ ------------------------------------

### Estimation {-}

The R code for parameter estimation begins by defining the `buildModPoly1` function
which can create the needed `dlm` model object from three parameters.

```{r, eval=FALSE}
buildModPoly1 <- function(v) {
  dV <- exp(v[1])
  dW <- exp(v[2])
  m0 <- v[3]
  dlmModPoly(1, dV=dV, dW=dW, m0=m0)
}
```

The R function itself takes one parameter, a 3-element vector,
into which the model parameters are packed.
The first two parameters are log-variance, not variance,
to prevent the optimizer from exploring negative values for variance.

To start, we need some reasonable guesses at the parameters.
They don't need to be perfect, but being in the right ballpark is useful.
```{r, eval=FALSE}
varGuess <- var(diff(y), na.rm=TRUE)
mu0Guess <- as.numeric(y[1])
```

The `dlmMLE` function finds the maximum likelihood estimate of the parameters,
starting with our reasonable guesses and
repeatedly calling our `buildModPoly1` until it converges on the MLE solution.
Always check for convergence.

```{r, eval=FALSE}
parm <- c(log(varGuess), log(varGuess), mu0Guess)
mle <- dlmMLE(y, parm=parm, buildModPoly1)

if (mle$convergence != 0) stop(mle$message)
```

From the MLE parameter estimates, we can build the final model.
```{r, eval=FALSE}
model <- buildModPoly1(mle$par)
```

### Example {-}

```{r, eval=TRUE}
library(dlm)

y <- datasets::Nile

buildModPoly1 <- function(v) {
  dV <- exp(v[1])
  dW <- exp(v[2])
  m0 <- v[3]
  dlmModPoly(1, dV=dV, dW=dW, m0=m0)
}

varGuess <- var(diff(y), na.rm=TRUE)
mu0Guess <- as.numeric(y[1])

parm <- c(log(varGuess), log(varGuess), mu0Guess)

mle <- dlmMLE(y, parm=parm, buildModPoly1)
if (mle$convergence != 0) stop(mle$message)

model <- buildModPoly1(mle$par)

cat("Observational variance:", model$V, "\n",
    "Transitional variance:", model$W, "\n",
    "Initial level:", model$m0, "\n")
```