Lab 13: Forecasting med ARIMA

I de neste to labbene skal vi ta i bruk noen av tidsrekke-verktøyene vi har lært, og anvende dem til å skape en prognose, eller forecast. Labbene er basert på boken Forecasting: Principles and Practice og bruker den som referanse. Det kan være særlig nyttig å lese gjennom k.8 og k.9.

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.0.6     ✓ dplyr   1.0.4
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(fpp2)

## Loading required package: forecast

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Loading required package: fma

## Loading required package: expsmooth

Prognoser basert på ARIMA

AR modellene som vi har brukt er en del av en større klasse modeller som heter ARIMA (AR: Autoregressive, MA: Moving Average. Vi kommer tilbake til hva ‘I’ betyr.)

An autoregressive model har vi allerede sett.

En Moving Average model er også dynamisk, men her er det feilene i modellen (eller sjokkene, hvis vi tenker økonomisk) som skaper dynamikken, og ikke y verdien.

Vi kan skrive en MA med orden 1 (MA(1)) modell som:

\(y_t = c + \epsilon_t + \epsilon_{t-1}\)

En MA med orden 2 (MA(2)) modell som:

\(y_t = c + \epsilon_t + \epsilon_{t-1} + \epsilon_{t-2}\)

Økonomiske intuisjonøen for en sånn modell er at en sjokk i en periode kan påvirke ikke bare nåværende periode, men også senere perioder.

En liten interessant fakta: alle AR modeller kan også skrives som MA modeller. Og alle MA modeller kan også skrives som AR modeller. Men som en praktisk sak, velger de fleste å bruke en blanding av AR og MA komponenter til å best modellere en tidsrekke.

Vi bruker en eksempel fra Forecasting boka. Her bruker de en datasett over bestillinger av utstry for kraftindustrien. Datasettet er inkludert i fpp2 pakken, og finnes med navnet elecequip.

str(elecequip)

##  Time-Series [1:195] from 1996 to 2012: 79.3 75.8 86.3 72.6 74.9 ...

Vi bruker autoplot funksjonen til fpp2 pakken til å plotte vår dataserie.

autoplot(elecequip)

Hva ser vi i denne dataserien?

Det første vi kanskje ser er at det virker som det er noen sesonger, eller faste opp- og ned-mønsteret.

Så først ønsker vi å sesongjustere datasettet, så at vi kun har den mer langsiktige elementet.

stl deler opp tidsrekken i flere komponenter: en sesong komponent, en trend komponent og en feil-komponent (seasonal, trend and irregular)

decomp = elecequip %>% stl(s.window="periodic")

head(decomp$time.series)

##           seasonal    trend  remainder
## Jan 1996 -5.444213 80.35590  4.4383098
## Feb 1996 -5.967189 80.27515  1.4720416
## Mar 1996  8.039247 80.19439 -1.9136391
## Apr 1996 -6.385339 80.16226 -1.1769171
## May 1996 -4.839102 80.13012 -0.4310188
## Jun 1996  7.732144 80.13135 -4.0534979

Vi sender deretter vår objekt med alle komponentene til funksjonen seasadj, som justerer den opprinnelige datasettet for sesongvariasjon (trekker fra sesongene).

eeadj = decomp %>% seasadj()

autoplot(eeadj)

Etter sesongene er tatt bort, ser vi at det er store svingninger i bestillinger, som beveger seg med makroøkonomien. (Kraftustyr er en form av investering, så dette er i tråd med teorien.)
Vi ser også at dataserien ser ikke stasjonær ut.
Vi kan bekrefte dette ved å kjøre en dickey-fuller test

library(tseries)
adf.test(eeadj)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  eeadj
## Dickey-Fuller = -2.4033, Lag order = 5, p-value = 0.4073
## alternative hypothesis: stationary

Vi kan isteden ta en first-difference å se om det løser vårt problem

d_eeadj = eeadj %>% diff()
adf.test(d_eeadj)

## Warning in adf.test(d_eeadj): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  d_eeadj
## Dickey-Fuller = -4.2461, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Nå ser det bedre ut!

Vi kan deretter se på ACF og pACF for å få litt innsikt i dynamikken

acf(d_eeadj)

pacf(d_eeadj)

Det ser ut som det er dynamikk i serien, og på pacf så ser det ut som det er kanskje 3 signifikante lag.

Så en første modell vi kan prøve er AR(3) med en differansse (I(1)). Det vil si en ARIMA(3,1,0).

diff y_t = diff y_t-1 + diff y_t-2 + diff y_t-3

fit1 = Arima(eeadj, order=c(3,1,0))

Vi kan sjekke residualene for å se om det er noe autokorrelasjon igjen?

acf(residuals(fit1))

pacf(residuals(fit1))

Det ser bra ut!

Nå kan vi snu modellen rundt og skape en forecast:

forecast1 = forecast(fit1, h=20)

autoplot(forecast1)

Oppgaver:

Hva er noen av problemene ved denne prognosen, fra et praktisk standspunkt?
Kjør en ARIMA modell uten differanse I(0). Hvordan ser prognosen ut da?
I boken bruker de en ARIMA(3,1,1) isteden. Hvorfor har de valgt denne modellen? Blir resultatene annerledes?

Du kan lese mer om nøyaktig hvordan prognosen er laget her

Prognoser med sesong.

I eksemepelen over, lagde vi en prognose ut i fra en serie uten sesonger. Nå skal vi se på en prognose med sesonger.

Vi skal se på data om detaljhandel fra europa. Dette kommer også med pakken fpp2 og heter euretail.

str(euretail)

##  Time-Series [1:64] from 1996 to 2012: 89.1 89.5 89.9 90.1 89.2 ...

Jeg bruker en funksjon i fpp2 som automatisk viser både datasettet pluss ACF og PACF

euretail %>%  ggtsdisplay()

Her har vi data som er både ikke-stasjonær og som ser ut som de har noen sesonger i det. Dette er kvartalvis data, og vi kan ta differansen årlig (hver 4. kvartal)

euretail %>% diff(lag=4) %>% ggtsdisplay()

Datasettet ser fortsatt ikke-stasjonær ut, så vi gjør en ekstra differanse:

euretail %>% diff(lag=4) %>% diff() %>% ggtsdisplay()

Nå begynner det å se bra ut.

Vi kjører vår ARIMA model. Hvis vi ser på vår PACF, så ser det ut som vi har autokorrelasjon i den første laggen. Vi har også sesong-korrelasjon på den 4. laggen.

Vi kan modellere det på følgende måte: 1 AR ledd, 1 sesong (4 kvartaler) ledd og både en sesong differanse og en vanlig differanse.

fit2 = euretail %>%
  Arima(order=c(1,1,0), seasonal=c(1,1,0))
residuals(fit2) %>% ggtsdisplay()

Nå ser vi ikke mange signifikante autokorrelasjoner

Nå er vi klar til å skape en prognose 3 år (12 kvartaler) i framtiden:

forecast2 = fit2 %>% forecast(h=12)
autoplot(forecast2)

Oppgave

Fra labb 12 lastet vi inn kronekurs data fra Norges Bank:

kronekurs = read_csv("https://jmaurit.github.io/anvendt_macro/lab/data/kroneKurs.csv")

## Parsed with column specification:
## cols(
##   date = col_date(format = ""),
##   NOK_USD = col_double()
## )

tail(kronekurs)

## # A tibble: 6 x 2
##   date       NOK_USD
##   <date>       <dbl>
## 1 2020-04-01   10.4 
## 2 2020-05-01   10.1 
## 3 2020-06-01    9.53
## 4 2020-07-01    9.30
## 5 2020-08-01    8.94
## 6 2020-09-01    9.14

kk_ts = ts(kronekurs$NOK_USD, start=c(1960,1), frequency=12)
autoplot(kk_ts)

Du skal lage en prognosemodell for kronen basert på en ARIMA modell.

Ved å se på tidsrekken, ville det vært fornuftig å avgrense perioden som du bruker til å lage din prognose? Hvorfor?
Skap en passende ARIMA modell og bruk den til å skape en prognose 12 måneder fram for valutakurset.
Tror du at du kunne skape en bedre prognose ved å bruke mer informasjon. Hva slaks informasjon ville du ønsket å bruke?

Lab 13: Forecasting med ARIMA

Johannes Mauritzen

Høsten 2020

Prognoser basert på ARIMA

Prognoser med sesong.

Oppgave