# cross correlation function besteedbaar inkomen (CBS) en boekverkoop (KVB)
# Frank Huysmans | WareKennis onderzoek en advies, Den Haag
# februari 2018
# https://wp.me/p7dgyC-Fz


# benodigde pakketten laden
library(tseries)
library(forecast)

# data entry
# reeksen zijn procentuele daling/stijging tov jaar ervoor voor de jaren 2007-2016
# inkomen_1 is lag inkomen (-1)
afzet   <- c(2.101,4.733,-2.947,0.202,-5.253,-4.478,-8.482,-9.024,4.558,4.359)
omzet   <- c(5.369,2.070,0.156,-2.991,-4.448,-6.402,-5.745,-7.619,2.680,4.618)
inkomen <- c(3.100,1.400,1.900,-0.500,-0.700,-1.100,-1.100,1.900,1.300,2.700)

# aanmaak lag(inkomen, -1) voor gebruik in regressie (zie verder)
inkomen_1 <- c(NA,3.100,1.400,1.900,-0.500,-0.700,-1.100,-1.100,1.900,1.300)

# tijdreeksen
ts_afzet   <- ts(afzet, frequency=1, start=2007)
ts_omzet   <- ts(omzet, frequency=1, start=2007)
ts_inkomen <- ts(inkomen, frequency=1, start=2007)

# inspectie
plot.ts(ts_afzet)
plot.ts(ts_omzet)
plot.ts(ts_inkomen)

# cross correlation functions
adf.test(ts_afzet)
kpss.test(ts_afzet)

adf.test(ts_omzet)
kpss.test(ts_omzet)

adf.test(ts_inkomen)
kpss.test(ts_inkomen)

ccf(ts_inkomen, ts_afzet, ylab = "cross-correlation", plot=TRUE)
ccf(ts_inkomen, ts_afzet, ylab = "cross-correlation", plot=FALSE)
ccf(ts_inkomen, ts_omzet, ylab = "cross-correlation", plot=TRUE)
ccf(ts_inkomen, ts_omzet, ylab = "cross-correlation", plot=FALSE)

# interpretatie cf. documentatie aif: The lag k value returned by ccf(x, y) estimates 
# the correlation between x[t+k] and y[t]. 
# dus (hoogste) correlatie bij -1 betekent dat inkomen 2007 en afzet/omzet 2008
# het sterkst zijn gecorreleerd, gevolgd door correlatie bij 0 (in hetzelfde jaar)

# dynamische regressie (zie https://www.otexts.org/fpp/9/1 example 9.3)
fit_afzet <- auto.arima(afzet, xreg=inkomen, stationary=TRUE, seasonal=FALSE)
fit_afzet 

fit_omzet <- auto.arima(omzet, xreg=inkomen, stationary=TRUE, seasonal=FALSE)
fit_omzet 

# auto.arima geeft in beide gevallen een (0,0,0)-model als best passend (AIC)
# oftewel: modellering van autocorrelaties is niet nodig en reeksen zijn stationair
# (zie ook tests hierboven)

# regressie afzet op lag(inkomen, -1): lag-variabele geeft sterkste Rˆ2 e
# en acceptabele residuen (m.n op observatie 8 = 2014)
fit_afzet <- lm(afzet ~ inkomen_1)
summary(fit_afzet)
coefficients(fit_afzet)
confint(fit_afzet, level=0.95)
anova(fit_afzet)
fitted(fit_afzet)
influence(fit_afzet)

# regressie omzet op lag(inkomen, -1)
fit_omzet <- lm(omzet ~ inkomen_1)
summary(fit_omzet)
coefficients(fit_omzet)
confint(fit_omzet, level=0.95)
anova(fit_omzet)
fitted(fit_omzet)
influence(fit_omzet)