Forecasting-4-1---Exp-Smoothing.html

<!DOCTYPE html>
<html lang="" xml:lang="">
  <head>
    <title>Forecasting-4-1---Exp-Smoothing.utf8</title>
    <meta charset="utf-8" />
    <link rel="stylesheet" href="my-theme.css" type="text/css" />
  </head>
  <body>
    <textarea id="source">


layout: true

.hheader[&lt;a href="index.html"&gt;&lt;svg style="height:0.8em;top:.04em;position:relative;fill:steelblue;" viewBox="0 0 576 512"&gt;&lt;path d="M280.37 148.26L96 300.11V464a16 16 0 0 0 16 16l112.06-.29a16 16 0 0 0 15.92-16V368a16 16 0 0 1 16-16h64a16 16 0 0 1 16 16v95.64a16 16 0 0 0 16 16.05L464 480a16 16 0 0 0 16-16V300L295.67 148.26a12.19 12.19 0 0 0-15.3 0zM571.6 251.47L488 182.56V44.05a12 12 0 0 0-12-12h-56a12 12 0 0 0-12 12v72.61L318.47 43a48 48 0 0 0-61 0L4.34 251.47a12 12 0 0 0-1.6 16.9l25.5 31A12 12 0 0 0 45.15 301l235.22-193.74a12.19 12.19 0 0 1 15.3 0L530.9 301a12 12 0 0 0 16.9-1.6l25.5-31a12 12 0 0 0-1.7-16.93z"/&gt;&lt;/svg&gt;&lt;/a&gt;]

---

class: center, middle, inverse
# Forecasting Time Series
## Simple Exponential Smoothing

.futnote[Eli Holmes, UW SAFS]

.citation[eeholmes@uw.edu]

---

## Naive forecast

For a naive forecast of the anchovy catch in 1988, we just use the 1987 catch.

`$$\hat{x}_{1988} = x_{1987}$$`
Which is the same as saying that we put 100% of the 'weight' on the most recent value and no weight on any value prior to that.

`$$\hat{x}_{1988} = 1 \times x_{1987} + 0 \times x_{1986} + 0 \times x_{1985} + \dots$$`
---

Past values in the time series have information about the current state, but only the most recent past value.

&lt;img src="Forecasting-4-1---Exp-Smoothing_files/figure-html/unnamed-chunk-1-1.png" style="display: block; margin: auto;" /&gt;


---

That's a bit extreme.  Often the values prior to the last value also have some information about future states.  But the 'information content' should decrease the farther in the past that we go.

&lt;img src="Forecasting-4-1---Exp-Smoothing_files/figure-html/unnamed-chunk-2-1.png" style="display: block; margin: auto;" /&gt;

---

Simple exponential smoothing uses this type of weighting that falls off exponentially and the objective to estimate the best weighting ( `\(\alpha\)` ):

&lt;img src="Forecasting-4-1---Exp-Smoothing_files/figure-html/unnamed-chunk-3-1.png" style="display: block; margin: auto;" /&gt;

---

## Fitting exponential smoothing models

The forecast package will fit a wide variety of exponential smoothing models.  The main fitting function is `ets()`:

```
ets(y, model = "ZZZ", &lt; + many other arguments &gt;)
```

* y : your data.  A time series of responses.

* model: what type of exponential smoothing model.

---

We are going to `ets()` to fit three simple types of exponential smoothing models:

model  | "ZZZ" | alternate function |
------------- | ------------- | --------- |
exponential smoothing no trend | "ANN" | `ses()` |
exponential smoothing with trend  | "AAN" | `holt()` |
exponential smoothing choose trend  | "AZN" | NA |

The alternate function does exactly the same fitting.  It is just a 'shortcut'.

---

## Fit exponential smoothing with no trend

This is like the naive model that just uses the last value to make the forecast, but instead of only using the last value it will use values farther in the past also. The weighting fall off exponentially.

Load the data and **forecast** package.


```r
data(greeklandings, package="FishForecast")
anchovy &lt;- subset(greeklandings, 
                  Species=="Anchovy" &amp; Year&lt;=1987)$log.metric.tons
anchovy &lt;- ts(anchovy, start=1964)
library(forecast)
```


```r
fit &lt;- ets(anchovy, model="ANN")
fr &lt;- forecast(fit, h=5)
```

---


```r
plot(fr)
```

&lt;img src="Forecasting-4-1---Exp-Smoothing_files/figure-html/unnamed-chunk-4-1.png" style="display: block; margin: auto;" /&gt;

---

## Look at the estimates


```r
fit
```

```
## ETS(A,N,N) 
## 
## Call:
##  ets(y = anchovy, model = "ANN") 
## 
##   Smoothing parameters:
##     alpha = 0.7065 
## 
##   Initial states:
##     l = 8.5553 
## 
##   sigma:  0.2166
## 
##       AIC      AICc       BIC 
##  6.764613  7.964613 10.298775
```

---

## The weighting function

&lt;img src="Forecasting-4-1---Exp-Smoothing_files/figure-html/ann.weighting-1.png" style="display: block; margin: auto;" /&gt;

---

## Produce forecast with ets from a previous fit

Say you want to estimate a forecasting model from one dataset and use that model to forecast another dataset or another area.  Here is how to do that.

This is the fit to the 1964-1987 data:


```r
fit1 &lt;- ets(anchovy, model="ANN")
```

Use that model with the 2000-2007 data and produce a forecast:


```r
dat &lt;- subset(landings, Species=="Anchovy" &amp; Year&gt;=2000 &amp; Year&lt;=2007)
dat &lt;- ts(dat$log.metric.tons, start=2000)
fit2 &lt;- ets(dat, model=fit1)
```

```
## Model is being refit with current smoothing parameters but initial states are being re-estimated.
## Set 'use.initial.values=TRUE' if you want to re-use existing initial values.
```

```r
fr2 &lt;- forecast(fit2, h=5)
```

---


```r
plot(fr2)
```

&lt;img src="Forecasting-4-1---Exp-Smoothing_files/figure-html/unnamed-chunk-7-1.png" style="display: block; margin: auto;" /&gt;

---

## Naive model with drift

Fit a model that uses the last observation as the forecast but includes a trend estimated from ALL the data.  This is what the naive model with drift does.


```r
fit.rwf &lt;- Arima(anchovy, order=c(0,1,0), include.drift=TRUE)
fr.rwf &lt;- forecast(fit.rwf, h=5)
plot(fr.rwf)
```

&lt;img src="Forecasting-4-1---Exp-Smoothing_files/figure-html/unnamed-chunk-8-1.png" style="display: block; margin: auto;" /&gt;

---

The trend seen in the blue line is estimated from the overall trend in ALL the data.


```r
coef(fit.rwf)
```

```
##      drift 
## 0.06577281
```

The trend from all the data is (last-first)/(number of steps).


```r
mean(diff(anchovy))
```

```
## [1] 0.06577281
```

So we only use the latest data to choose the level for our forecast but use all the data to choose the trend?  It would make more sense to weight the more recent trends more heavily.

---

## Exponential smoothing model with trend

The exponential smoothing model with trend does this.  The one-year trend is 
`$$x_t - x_{t-1}$$`
That is how much the data increased or decreased.  

---


```r
plot(diff(anchovy),ylim=c(-0.3,.3))
abline(h=0, col="blue")
abline(h=mean(diff(anchovy)),col="red")
title("0 means no change")
```

&lt;img src="Forecasting-4-1---Exp-Smoothing_files/figure-html/unnamed-chunk-11-1.png" style="display: block; margin: auto;" /&gt;

---

If we take the average of all `\(x_t - x_{t-1}\)` we are using the average trend like the naive model with drift.  We put an equal weighting on all trends.

But we could use a weighting that falls off exponentially so that we more recent trends affect the forecast more than trends in the distant past.  That is what an exponential smoothing model with trend does.

&lt;img src="Forecasting-4-1---Exp-Smoothing_files/figure-html/unnamed-chunk-12-1.png" style="display: block; margin: auto;" /&gt;

---

## Naive model with trend

If your training data are length `\(T\)`, then a forecast for `\(T+h\)` is

`$$\hat{x}_{T+h} = l_T + h \bar{b}$$`

where `\(\hat{b}\)` is the mean of the the yearly changes in `\(x\)`, so the mean of `\(x_t - x_{t-1}\)`.

`$$\hat{b} = \sum_{t=2}^T (x_t - x_{t-1})$$`

## Exponential smoothing model with trend

`$$\hat{x}_{T+h} = l_T + h b_T$$`

where `\(b_T\)` is a weighted average with the more recent trends given more weight.

`$$b_t = \sum_{t=2}^T \beta (1-\beta)^{t-2}(x_t - x_{t-1})$$`
---
 
## Fit exponential smoothing with trend


```r
fit &lt;- ets(anchovy, model="AAN")
fr &lt;- forecast(fit, h=5)
plot(fr)
```

&lt;img src="Forecasting-4-1---Exp-Smoothing_files/figure-html/unnamed-chunk-13-1.png" style="display: block; margin: auto;" /&gt;

---

## Decomposing your model fit

Sometimes you would like to see the smoothed level and smoothed trend that the model estimated. You can see that with `plot(fit)` or `autoplot(fit)`.


```r
autoplot(fit)
```

&lt;img src="Forecasting-4-1---Exp-Smoothing_files/figure-html/unnamed-chunk-14-1.png" style="display: block; margin: auto;" /&gt;

---

## Forecast performance

We can evaluate the forecast performance with forecasts of our test data or we can use all the data and use time-series cross-validation.

Let's start with the former.


---

## Test forecast performance against a test data set

We will fit an an exponential smoothing model with trend to the training data and make a forecast for the years that we 'held out'. 


```r
fit1 &lt;- ets(traindat, model="AAN")
h=length(testdat)
fr &lt;- forecast(fit1, h=h)
plot(fr)
points(end(traindat)[1]+1:h, testdat, pch=2, col="red")
legend("topleft", c("forecast","actual"), pch=c(20,2), col=c("blue","red"))
```

&lt;img src="Forecasting-4-1---Exp-Smoothing_files/figure-html/unnamed-chunk-16-1.png" style="display: block; margin: auto;" /&gt;

---

We can calculate a variety of forecast error metrics with


```r
accuracy(fr, testdat)
```

```
##                      ME      RMSE       MAE        MPE     MAPE      MASE
## Training set  0.0155561 0.1788989 0.1442712  0.1272938 1.600532 0.7720807
## Test set     -0.5001701 0.5384355 0.5001701 -5.1678506 5.167851 2.6767060
##                      ACF1 Theil's U
## Training set -0.008371542        NA
## Test set     -0.500000000  2.690911
```

We would now repeat this for all the models in our candidate set and choose the model with the best forecast performance.

---

## Forecast performance with time-series cross-validation

Another approach is to use all the data and test a series of forecasts made by fitting the model to different lengths of the data.

In this approach, we don't have test data.  Instead we will use all the data for fitting and for forecast testing.

We will redefine `traindat` as all our Anchovy data.


---

## tsCV() function

We will use the `tsCV()` function.  We need to define a function that returns a forecast.


```r
far2 &lt;- function(x, h, model){
  fit &lt;- ets(x, model=model)
  forecast(fit, h=h)
  }
```

---

Now we can use `tsCV()` to run our `far2()` function to a series of training data sets.  We will specify that a NEW ets model be estimated for each training set.  We are not using the weighting estimated for the whole data set but estimating the weighting new for each set.

The `e` are our forecast errors for all the forecasts that we did with the data.


```r
e &lt;- tsCV(traindat, far2, h=1, model="AAN")
e
```

```
## Time Series:
## Start = 1 
## End = 26 
## Frequency = 1 
##  [1] -0.245378390  0.366852341  0.419678595 -0.414861770 -0.152727933
##  [6] -0.183775208 -0.013799590  0.308433377 -0.017680471 -0.329690537
## [11] -0.353441463  0.266143346 -0.110848616 -0.005227309  0.157821831
## [16]  0.196184446  0.008135667  0.326024067  0.085160559  0.312668447
## [21]  0.246437781  0.117274740  0.292601670 -0.300814605 -0.406118961
## [26]           NA
```

---

Let's look at the first few `e` so we see exactly with `tsCV()` is doing.


```r
e[2]
```

```
## [1] 0.3668523
```

This uses training data from `\(t=1\)` to `\(t=2\)` so fits an ets to the first two data points alone.  Then it creates a forecast for `\(t=3\)` and compares that forecast to the actual value observed for `\(t=3\)`.


```r
TT &lt;- 2 # end of the temp training data
temp &lt;- traindat[1:TT]
fit.temp &lt;- ets(temp, model="AAN")
fr.temp &lt;- forecast(fit.temp, h=1)
traindat[TT+1] - fr.temp$mean
```

```
## Time Series:
## Start = 3 
## End = 3 
## Frequency = 1 
## [1] 0.3668523
```

---


```r
e[3]
```

```
## [1] 0.4196786
```

This uses training data from `\(t=1\)` to `\(t=2\)` so fits an ets to the first two data points alone.  Then it creates a forecast for `\(t=3\)` and compares that forecast to the actual value observed for `\(t=3\)`.


```r
TT &lt;- 3 # end of the temp training data
temp &lt;- traindat[1:TT]
fit.temp &lt;- ets(temp, model="AAN")
fr.temp &lt;- forecast(fit.temp, h=1)
traindat[TT+1] - fr.temp$mean
```

```
## Time Series:
## Start = 4 
## End = 4 
## Frequency = 1 
## [1] 0.4196786
```

---

## Testing a specific ets model

By specifying `model="AAN"`, we estimated a new ets model (meaning new weighting) for each training set used.  We might want to specify that we use only the weighting we estimated for the full data set.

We do this by passing in a fit to `model`.

The `e` are our forecast errors for all the forecasts that we did with the data. `fit1` below is the ets estimated from all the data 1964 to 1989.  Note, the code will produce a warning that it is estimating the initial value and just using the weighting.  That is what we want.


```r
fit1 &lt;- ets(traindat, model="AAN")
e &lt;- tsCV(traindat, far2, h=1, model=fit1)
e
```

```
## Time Series:
## Start = 1 
## End = 26 
## Frequency = 1 
##  [1]           NA  0.576663901  1.031385937  0.897828249  1.033164616
##  [6]  0.935274283  0.958914499  1.265427119 -0.017241938 -0.332751184
## [11] -0.330473144  0.255886314 -0.103926617  0.031206730  0.154727479
## [16]  0.198328366 -0.020605522  0.297475742  0.005297401  0.264939892
## [21]  0.196256334  0.129798648  0.335887872 -0.074017535 -0.373267163
## [26]           NA
```

---

## Forecast accuracy metrics

Now we can compute forecast accuracy metrics from the forecast errors (`e`).

RMSE: root mean squared error

```r
rmse &lt;- sqrt(mean(e^2, na.rm=TRUE))
```

MAE: mean absolute error

```r
mae &lt;- mean(abs(e), na.rm=TRUE)
```

We would repeat this process for all models in our candidate set and compare forecast fits to choose our forecast model.

---

## Next steps

Now you are ready to compare forecasts from a variety of models and choose a forecast model.  Note when you compare models, you can use both 'training data/test data' and use time-series cross-validation, but report the metrics in separate columns.  Example, 'RMSE from tsCV' and 'RMSE from test data'.

---

## Example candidate set for anchovy

* Exponential smoothing model with trend
```
fit &lt;- ets(traindat, model="AAN")
fr &lt;- forecast(fit, h=1)
```
* Exponential smoothing model no trend
```
fit &lt;- ets(traindat, model="ANN")
fr &lt;- forecast(fit, h=1)
```
* ARIMA(0,1,1) with drift (best)
```
fit &lt;- Arima(traindat, order(0,1,1), include.drift=TRUE)
fr &lt;- forecast(fit, h=1)
```
* ARIMA(2,1,0) with drift (within 2 AIC of best)
```
fit &lt;- Arima(traindat, order(2,1,0), include.drift=TRUE)
fr &lt;- forecast(fr)
```

---

## Candidata models continued

* Time-varying regression with linear time
```
traindat$t &lt;- 1:24
fit &lt;- lm(log.metric.tons ~ t, data=traindat)
fr &lt;- forecast(fit, newdata=data.frame(t=25))
```

* Naive no trend
```
fit &lt;- Arima(traindat, order(0,1,0))
fr &lt;- forecast(fit, h=1)
# or simply
fr &lt;- rwf(traindat)
```
* Naive with trend
```
fit &lt;- Arima(traindat, order(0,1,0), include.drift=TRUE)
fr &lt;- forecast(fit)
# or simply
fr &lt;- rwf(traindat, drift=TRUE)
```
* Average or mean
```
fit &lt;- Arima(traindat, order(0,0,0))
fr &lt;- forecast(fit)
```
    </textarea>
<style data-target="print-only">@media screen {.remark-slide-container{display:block;}.remark-slide-scaler{box-shadow:none;}}</style>
<script src="https://remarkjs.com/downloads/remark-latest.min.js"></script>
<script>var slideshow = remark.create({
"highlightStyle": "github",
"highlightLines": true
});
if (window.HTMLWidgets) slideshow.on('afterShowSlide', function (slide) {
  window.dispatchEvent(new Event('resize'));
});
(function(d) {
  var s = d.createElement("style"), r = d.querySelector(".remark-slide-scaler");
  if (!r) return;
  s.type = "text/css"; s.innerHTML = "@page {size: " + r.style.width + " " + r.style.height +"; }";
  d.head.appendChild(s);
})(document);

(function(d) {
  var el = d.getElementsByClassName("remark-slides-area");
  if (!el) return;
  var slide, slides = slideshow.getSlides(), els = el[0].children;
  for (var i = 1; i < slides.length; i++) {
    slide = slides[i];
    if (slide.properties.continued === "true" || slide.properties.count === "false") {
      els[i - 1].className += ' has-continuation';
    }
  }
  var s = d.createElement("style");
  s.type = "text/css"; s.innerHTML = "@media print { .has-continuation { display: none; } }";
  d.head.appendChild(s);
})(document);
// delete the temporary CSS (for displaying all slides initially) when the user
// starts to view slides
(function() {
  var deleted = false;
  slideshow.on('beforeShowSlide', function(slide) {
    if (deleted) return;
    var sheets = document.styleSheets, node;
    for (var i = 0; i < sheets.length; i++) {
      node = sheets[i].ownerNode;
      if (node.dataset["target"] !== "print-only") continue;
      node.parentNode.removeChild(node);
    }
    deleted = true;
  });
})();
// adds .remark-code-has-line-highlighted class to <pre> parent elements
// of code chunks containing highlighted lines with class .remark-code-line-highlighted
(function(d) {
  const hlines = d.querySelectorAll('.remark-code-line-highlighted');
  const preParents = [];
  const findPreParent = function(line, p = 0) {
    if (p > 1) return null; // traverse up no further than grandparent
    const el = line.parentElement;
    return el.tagName === "PRE" ? el : findPreParent(el, ++p);
  };

  for (let line of hlines) {
    let pre = findPreParent(line);
    if (pre && !preParents.includes(pre)) preParents.push(pre);
  }
  preParents.forEach(p => p.classList.add("remark-code-has-line-highlighted"));
})(document);</script>

<script>
(function() {
  var links = document.getElementsByTagName('a');
  for (var i = 0; i < links.length; i++) {
    if (/^(https?:)?\/\//.test(links[i].getAttribute('href'))) {
      links[i].target = '_blank';
    }
  }
})();
</script>

<script>
slideshow._releaseMath = function(el) {
  var i, text, code, codes = el.getElementsByTagName('code');
  for (i = 0; i < codes.length;) {
    code = codes[i];
    if (code.parentNode.tagName !== 'PRE' && code.childElementCount === 0) {
      text = code.textContent;
      if (/^\\\((.|\s)+\\\)$/.test(text) || /^\\\[(.|\s)+\\\]$/.test(text) ||
          /^\$\$(.|\s)+\$\$$/.test(text) ||
          /^\\begin\{([^}]+)\}(.|\s)+\\end\{[^}]+\}$/.test(text)) {
        code.outerHTML = code.innerHTML;  // remove <code></code>
        continue;
      }
    }
    i++;
  }
};
slideshow._releaseMath(document);
</script>
<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
(function () {
  var script = document.createElement('script');
  script.type = 'text/javascript';
  script.src  = 'https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML';
  if (location.protocol !== 'file:' && /^https?:/.test(script.src))
    script.src  = script.src.replace(/^https?:/, '');
  document.getElementsByTagName('head')[0].appendChild(script);
})();
</script>
  </body>
</html>