In mijn vorige blog heb ik uitgelegd hoe je met hulp van de forecast package een voorspelling kan maken met een bestaande tijdreeks. Hier werd gebruik gemaakt van auto.arima, een algoritme dat automatisch het best verklarende Arima model selecteerd. Hoewel dit algoritme enorm handig is, vooral voor de beginnende ‘forecaster’, is het ook van belang om het proces te begrijpen waar auto.arima gebruik van maakt. Dit gaan we in deze blog doen door zelf een tijdreeks te analyseren met behulp van (handmatige)Arima.

Het doel van deze blog is om meer kennis te krijgen over Arima door begrijpbare voorbeelden en uitleg. Je gaat bekend worden met :

In deze blog wordt gebruik gemaakt van de package Forecast van Rob J. Hyndman. Rob Hyndman is Professor of Statistics en schrijver van het online textboek Forecasting: principles and practice.

ARMA model

In mijn vorige blog is meer informatie te vinden over het maken van een tijdreeks, het ontleden hiervan en informatie over Arima.

ARMA modellen zijn ontwikkeld voor het analyseren van tijdreeksen. Hierin wordt rekening gehouden met afhankelijkheid van tijdpunten. Punten zijn niet onafhankelijk van elkaar. - Box en Jenkins (1976)

ARMA staat voor AutoRegressive Moving Average.

‘AutoRegressive’(p) is het procesgedeelte waar met behulp van een gewogen som van waarden geschat wordt wat de toekomstige waarden zullen zijn. Met andere woorden, dit gedeelte van het model geeft aan in welke mate vorige waarden een effect hebben op de huidige waarden.

Een moving average(q), in het Nederlands voortschrijdend gemiddelde, heeft als belangrijkste eigenschap een onderliggende trend te herkennen. Dit onderdeel meet de in tijd verschoven waarden van de residuals . Dit gedeelte neemt dus een gedeelte van de foutenmarge in de vorige voorspelling mee in de voorspelling van de volgende waarde.

Stationary series

Een tijdreeks - ARMAmodel kan alleen gebouwd worden indien de data stationaroftwel stabiel zijn. Met stationair wordt bedoeld dat de eigenschappen van de data onafhankelijk zijn van de tijd. Met andere woorden: er bevindt zich geen trend en seasonality in de tijd. Daarnaast dient het gemmidelde, variantie en autocorrelatie gelijkmatig te zijn over de tijd heen.

Definitie volgens Hyndman; If {yt} is a stationary time series, then for all s, the distribution of (yt , . . . , yt+s) does not depend on t.

In een stationaire tijdreeks ga je er van uit de statistische kenmerken voor de toekomst het zelfde zullen zijn en dit maakt het voorspellen makkelijker.

Waarom moet een dataset stationair zijn ?

Indien de tijdreeks niet stationair is kan het zijn dat het gemiddelde niet gelijkmatig over de tijd heen en dit kan er voor zorgen dat de forecast niet betrouwbaar is. Stel dat een tijdreekstrend blijft groeien, het gemiddelde wordt dan ook meer naarmate de tijd ‘vooruitgaat’. Bij het forecasten wordt het gemiddelde onderschat aangezien deze alleen naar de geschiedenis kijkt en dit brengt waarschijnlijk te lage voorspellingen voor de toekomst. Met andere woorden, hier wordt geen rekening gehouden met een groeiende gemiddelde. Daarnaast kan er bijvoorbeeld een relatie gesuggereerd worden tussen twee metingen, die niet aanwezig is.

Wanneer is een tijdreeks Stationar?

Een stationary tijdreeks lijkt op het eerste gezicht ‘horizontaal’
Er is sprake van een constante variantie
Geen patronen te herkennen voor in de toekomst.

Gelukkig zijn er verschillende testen die aangeven of een tijdreeks stationair is. En er zijn ook manieren om een tijdreeks stationair te maken.

Differentieren

De eerste stap van een ARMA model maken is het stationair maken van de data. Dit heet differentieren, hier wordt in het Arima model naar verwezen door het ‘Integrated’ gedeelte.

AutoRegressive integrated Moving Average.

I ‘Integrated’ (d) Het ‘integrated’ gedeelte van een ARIMA model zorgt ervoor dat de data stationary wordt, door te differenti?ren. Bij differenti?ren neem je het verschil van achtereenvolgende waarden. Op deze manier kunnen trends en seasonality verwijderd worden uit data. Indien de data nog steeds niet stationair zijn, kan dit nog een keer worden gedaan.

Het non-seasonal Arima model ziet er in notaties als volgt uit: ARIMA(p,d,q) .

= AR AutoRegressive
= I Integrated
= MA Moving Average

Na het vaststellen van het aantal keer dat er gedifferentieerd (d) zal worden (meestal niet meer dan 2 keer) worden de waarden voor AR(p) en MA(q) stapsgewijs bepaald.

Belangrijk om te onthouden is dat de dataset eerst gedifferentieerd moet worden indien nodig, daarna kunnen we pas verder gaan met de volgende stap.

Voorbeeld differentieren

Verschillende testen zijn er aanwezig om te controleren of een tijdreeks stationary is. Voorbeelden zijn de Augmented Dickey-Fuller test (ADF) of de Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test.

Installeer en laadt, indien je dat nog niet hebt gedaan, de fpp package van Hyndman.

#install.packages("fpp")
suppressMessages(require(fpp))

Unit Root test 1

Test of een dataset stationair is

Augmented Dickey-Fuller test (ADF). In deze test is de alternatieve hypothese dat de data stationary is. Een voorbeeld wordt gegeven aan de hand van de DowJones dataset. Een grote p-value geeft aan dat de data gedifferentieerd dient te worden.

Als grenswaarde wordt standaard p=0.05 gebruikt.

data(dj)
adf.test(dj,alternative = "stationary",k=12)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dj
## Dickey-Fuller = -2.1488, Lag order = 12, p-value = 0.5135
## alternative hypothesis: stationary

Unit Root test 2

Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test. Deze test heeft een andere null hypothese, namelijk dat de data stationary is. Een grote p-value (>0.05) geeft aan dat de data NIET gedifferentieerd dient te worden.

data(dj)
kpss.test(dj)

## Warning in kpss.test(dj): p-value smaller than printed p-value

## 
##  KPSS Test for Level Stationarity
## 
## data:  dj
## KPSS Level = 1.9812, Truncation lag parameter = 3, p-value = 0.01

In dit voorbeeld dient er gedifferentieerd te worden.

Hoe vaak differentieren?

Deze code rekent voor je uit hoevaak er gedifferentieerd dient te worden.

ndiffs(dj)

## [1] 1

Differentieren kan in R door diff te plaatsen voor de data

Is de data nu wel stationary?

adf.test(diff(dj),alternative = "stationary")

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(dj)
## Dickey-Fuller = -6.7623, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

kpss.test(diff(dj))

## 
##  KPSS Test for Level Stationarity
## 
## data:  diff(dj)
## KPSS Level = 0.050693, Truncation lag parameter = 3, p-value = 0.1

ARIMA(p,d,q)

Na het vaststellen van het aantal keer dat er gedifferentieerd (d) zal worden (meestal niet meer dan 2 keer) worden de waarden voor AR(p) en MA(q) stapsgewijs bepaald. Maar.. hoe doen we dat?

Autocorrelatie

Door te kijken naar de autocorrelatie met eerdere waarden (die we ‘lags’ noemen). Autocorrelation meet de ‘lagged values’ van een timeserie. Hierin wordt gekeken naar de correlatie met zichzelf in een eerder moment. In timeseries wordt er gesproken van ‘lags’. Lag 1 meet bijvoorbeeld de relatie tussen Yt en yt-1 en lag 2 meet de relatie Yt en Yt-2

Als een ‘spike’ boven de stippellijn (=CI) uit komt is er sprake van een AR of MA model.

par(mfrow=c(1,2))
acf(diff(dj))
pacf(diff(dj))

par(mfrow=c(1,1))

Tsdisplay geeft de originele data en de ACF+PACF plot weer

tsdisplay(diff(dj))

In de dj dataset is geen correlatie of andere duidelijk verband te zien met de ‘lags’ in de data. Dit wijst op een AR = 0 en MA = 0 model. Dit brengt ons op een Arima(0,1,0) model

Door visuale inspectie van de ACF en PACF plot kunnen we de waarde voor p en q invullen.

ACF vs PACF

In ACF wordt er gekeken naar de correlatie tussen yt-1 en yt-2, yt-3 en yt-4 etc. Echter, op deze manier correleert yt-3 ook met yt-1. Daarom bekijken we ook de PACF plot . De partial autocorrelations geven de correlatie weer tussen lag 1 en lag k terwijl het corrigeert voor de correlaties tussen de overige lags.

De ACF plot geeft ons informatie over het MA model terwijl de PACF ons informatie geeft over het AR model.

Hyndman :

De data volgt een ARIMA(p,d,0) model als de ACF- en PACF-plots van de gediffentieerde data het volgende patroon laten zien:

-De ACF is exponentieel aan het afnemen;

-Er is een significante spike bij lag p in PACF, maar vervolgens geen verdere spikes waar te nemen.

De data volgt een ARIMA(0,d,q) model als de ACF- en PACF-plots van de gediffentieerde data het volgende patroon laten zien:

-De PACF is exponentieel aan het afnemen;

-Er is een significante spike bij lag p in ACF, maar vervolgens geen verdere spikes waar te nemen.

Ardj<-Arima(dj, order=c(0,1,0))

Residuals analyse

Controleer of de residuals ongecontroleerd zijn met elkaar. Indien dit niet het geval is, is er meer uit de data te halen.

Autocorrelatie tussen residuals wijst op verborgen informatie.
Residuals plot moet eruit zien als ‘white noise’. Dit betekent dat er geen correlatie is met vorige waarden.
(Het gemiddelde van de residuals moet nul zijn, anders is er een bias) Indien het gemiddelde van de residuals hoger ( of lager) is dan 0, dan kan deze waarde bij de puntvoorspelling opgeteld worden.

Checkresiduals in Forecast

In de nieuwe forecast 8.0 package van Hyndman;

is een nieuwe functie ‘checkresiduals’, hiermee kan je de assumpties controleren zoals ongecorreleerdheid en een normale distributie.

Daarnaast voert deze ook een Ljung-Box test uit met de juiste lag en vrijheidsgraden.Een hoge p-waarde voor deze test suggereert dat de residuals white noise zijn.

Installeer minimaal forecastpackage version 8.0

packageDescription("forecast")$Version
install.packages("forecast")
require(forecast)

checkresiduals(Ardj)

## 
##  Ljung-Box test
## 
## data:  residuals
## Q* = 14.462, df = 10, p-value = 0.153
## 
## Model df: 0.   Total lags used: 10

NB Hoewel de residuals er uit zien als white noise, hoeft het niet altijd te betekenen dat het beste model is gevonden. Het model lijkt goed te zijn om te forecasten, maar het is de moeite waard om altijd nog even omliggende modellen te bekijken.

Overzicht Modelselectie

Een manier voor het selecteren van het beste model is verschillende ARIMA volgorden te proberen. Het model met de laagste AICc maat is het beste model!

Voor meer informatie over de AIC waarde zie “https://www.otexts.org/fpp/8/6”

model

Voorbeeld Non-seasonal Arima

Als voorbeeld maken we een voorbeeld arimamodel van de dataset consumption US. Deze dataset ziet er als volgt uit:

data(usconsumption)
plot(usconsumption[,1])

Stap 1: differentieren

Stap 2: plot de Acf en Pacf

ndiffs(usconsumption[,1])

## [1] 0

par(mfrow=c(1,2))
Acf(usconsumption[,1],main="")
Pacf(usconsumption[,1],main="")

par(mfrow=c(1,1))

In het consumptionmodel is er sprake van een Arima model met order (0,0,3). Dit is terug te zien in de afnemende groei van de spikes in de PACF plot.

arimaus<-Arima(usconsumption[,1],order=c(0,0,3))
arimaus

## Series: usconsumption[, 1] 
## ARIMA(0,0,3) with non-zero mean 
## 
## Coefficients:
##          ma1     ma2     ma3    mean
##       0.2542  0.2260  0.2695  0.7562
## s.e.  0.0767  0.0779  0.0692  0.0844
## 
## sigma^2 estimated as 0.3953:  log likelihood=-154.73
## AIC=319.46   AICc=319.84   BIC=334.96

Voorbeeld Seasonal ARIMA model

Voorbeeld SeizoenArima model met Car sales

Door een tsdisplay te plotten kan er gezien worden of er sprake is van seasonality. Een andere mogelijkheid is decompose.

#Laad data in R
data.orig<-read.csv(file="../monthly-car-sales-in-quebec-1960.csv")
#Gebruik alleen eerste 7 jaar 
data<-data.orig[1:96,]
myts <- ts(data$Monthly.car.sales.in.Quebec.1960.1968 ,start=c(1960, 1), end=c(1967,12), frequency=12)

Wat valt je op aan de ACF plot?

tsdisplay(myts)

Spikes na 12 maanden wijst op seizoenseffecten.

Hyndman;

Het seizoensgedeelte van het AR of MA model kan gezien worden in de seizoenslags in de PACF en ACF. Zo zal een ARIMA(0,0,0)(0,0,1)12 model een spike laten zien in lag 12 en zal de PACF een afnamen laten zien in lag 12, 24, 36, ..

Ook voor een seizoensArima kan het zijn dat er gedifferentieerd dient te worden. Dit kan met diff(x,12)

nsdiffs(myts)

## [1] 1

mytsdiff<-diff(myts,12) ## Diff 12 maanden (seizoen)
tsdisplay(mytsdiff)

kpss.test(mytsdiff)

## 
##  KPSS Test for Level Stationarity
## 
## data:  mytsdiff
## KPSS Level = 0.43458, Truncation lag parameter = 2, p-value =
## 0.06225

ndiffs(mytsdiff)

## [1] 0

In de ACF en PACF plot is te zien dat de jaarlijkste trend nu verwijderd is uit de data.

In de ACF plot is nu alleen een spike te zien bij lag 12 en niet meer op lag 24. Dit is een indicator van een seasonal MA model.

Verder is er nog een AR of MA 2 model. Heel duidelijk is niet te zien of het om een MA of AR gaat, dus we proberen beide en kiezen degene met de laagste AICc waarde. Nadat we de waarden voor het seizoens-arima model hebben geschat proberen we naast verschillende combinaties ook het toevoegen van een constante. De constante is bijvoorbeeld het toevoegen van een gemiddelde, of drift (mee’groeiende’ gemiddelde)

model1<-Arima(myts,order = c(2,0,0), seasonal= c(0,1,1))
model2<-Arima(myts,order = c(0,0,2), seasonal= c(0,1,1))
model1<-Arima(myts,order = c(2,0,0), seasonal= c(0,1,1),include.constant = TRUE)
model2<-Arima(myts,order = c(0,0,2), seasonal= c(0,1,1),include.constant = TRUE)
model1$aicc

## [1] 1480.302

Het laatste model lijkt de laagste AICc waarden te hebben. Vervolgens kunnen we deze forecasten.

plot(forecast(model1))

Conclusie

In deze blog ben ik dieper ingegaan op de onderliggende processen van het Arima model. Is dit toch nog te lastig? Bekijk dan nog een keer vorige blog , waar gebruik wordt gemaakt van de auto.arima.

Andere blogs;

In mijn volgende blog leg ik uit hoe je naast Arima modellen ook op andere manieren voorspellingen kan maken. Deze methoden zijn eenvoudiger maar niet altijd slechter.

Daarnaast ga ik uitleggen hoe je het forecast model kan testen en hoe je uit verschillende modellen de bestpassende kan kiezen.

Verder wordt er ook verteld welke tijdaggregatie je het beste kunt kiezen voor de dataset en hoe je verschillende tijdaggregaties kan combineren.

Tijdreeksanalyse ARIMA in R

Handleiding modelselectie in R

Mandy oud