La funzione di autocorrelazione parziale calcola la funzione di autocorrelazione escludendo l’effetto regressivo dei termini precedenti.
pacf(1:100)
pacf(sin(1:100 * pi/10)) # periodo 20
pacf(sin(1:100 * pi/5)) # periodo 10
pacf(sin(1:100 * pi/3)) # periodo 6
pacf(rnorm(1:100))Esaminiamo un esempio più elaborato. Tipicamente, una serie storica con carattere autoregressivo che dipende dai primi \(p\) istanti precedenti avrà una pacf nulla ovunque tranne che in corrispondenza delle dipendenze, e una acf che decade a zero rapidamente.
U1 = rep(0, 200)
U2a = rep(0, 200)
U2b = rep(0, 200)
U3 = rep(0, 200)
for (i in 4:200) {
U1[i] = 0.5 * U1[i - 2] + rnorm(1) # AR(1)
U2a[i] = 0.2 * U2a[i - 1] + 0.5 * U2a[i - 2] + rnorm(1) # AR(2)
U2b[i] = -0.2 * U2b[i - 1] + 0.5 * U2b[i - 2] + rnorm(1) # AR(2)
U3[i] = -0.2 * U3[i - 1] + 0.5 * U3[i - 2] + 0.55 * U3[i - 3] + rnorm(1) # AR(3)
}
layout(matrix(c(1, 2, 1, 3), 2, 2))
# AR(1)
ts.plot(U1)
acf(U1, 60)
pacf(U1)
# AR(2)
ts.plot(U2a)
acf(U2a, 60)
pacf(U2a)
# AR(2) (alternato)
ts.plot(U2b)
acf(U2b, 60)
pacf(U2b)
# AR(3)
ts.plot(U3)
acf(U3, 60)
pacf(U3)
layout(1)Consideriamo il dataset relativo alle temperature di Nottingham. Valutiamo l’autocorrelazione e l’autocorrelazione al netto dell’effetto regressivo. Il risultato sembra suggerire una dipendenza che arriva fino a \(7\) lag, forse \(11\).
Proviamo il modello regressivo fatto a mano con dipendenze fino a \(13\) lag.
L = length(nottem)
l = 13 # numero di lag in ingresso
mnt = matrix(nrow = L - l, ncol = l + 1)
for (i in 1:(l + 1)) {
mnt[, i] = nottem[i:(L - l - 1 + i)]
}
mnt <- data.frame(mnt)
nt.lm <- lm(X14 ~ ., data = mnt) # X14 perché 13 lag in ingresso
summary(nt.lm)Eliminando ad uno ad uno i fattori con \(\wp\)-value alto (nell’ordine: 9,4,5,8,12,7,11,2), c’è un guadagno in termini di \(R^2_{adj}\) (o alla peggio una perdita trascurabile) e un \(R^2\) pressoché invariato. Sottrarre i fattori 1 e 10, seppur di \(\wp\)-value basso, non altera in maniera significativa la varianza spiegata. Sottrarre invece il terzo o il tredicesimo fattore comporta un calo un po’ più significativo di varianza spiegata.
nt.lm <- lm(X14 ~ . - X9 - X4 - X5 - X8 - X12 - X7 - X11 - X2 - X6 - X1 -
X10, data = mnt)
summary(nt.lm)Per confrontare il risultato e la predizione a questo livello naive non possiamo utilizzare predict, perché alcuni dei fattori in input sono predizioni a tempi precedenti.
anni = 2
L = length(nottem)
pt = rep(0, L + 12 * anni)
pt[1:L] = nottem
for (i in 1:(12 * anni)) {
pt[L + i] = coef(nt.lm) %*% c(1, pt[L + i - 11], pt[L + i - 1])
}
nt.lm.pt = ts(pt, frequency = 12, start = c(1920, 1))
nt.lm.a = window(nottem, c(1921, 2)) - resid(nt.lm)
ts.plot(nottem, nt.lm.a, window(nt.lm.pt, c(1939, 12)), col = c("black",
"blue", "red"))La predizione precedente è stata ottenuta con la seguente modalità. Il modello regressivo fornisce i coefficienti \((b,a_1,a_2,\dots,a_\ell)\), dove \(\ell\) è il numero di lag in ingresso, che a parte l’intercetta corrispondono alle colonne X1, X2, …, \(X\ell\). La previsione è \[
x_{L+i}
= b + a_1 x_{L+i-\ell} + a_2 x_{L+i-\ell-1} + \dots + a_\ell x_{L+i-1},
\] ed il generico termine corrispondente al coefficiente \(a_m\) (\(m=1,2,\dots,\ell\)) è \(x_{L+i-(l-m+1)}\). Nella previsione dunque devono apparire i fattori corrispondenti alle colonne che sono state incluse nel modello. Ad esempio nella predizione precedente abbiamo incluso le colonne X13 e X3. Il termine per X3 (\(\ell=13\) e \(m=3\))è pt[L+i-11], il termine per X13 (\(\ell=13\), \(m=13\)) è pt[L+i-1].
Proviamo con il modello completo.
nt.lm <- lm(X14 ~ ., data = mnt)
anni = 2
pt = rep(0, L + 12 * anni)
pt[1:L] = nottem
for (i in 1:(12 * anni)) {
pt[L + i] = coef(nt.lm) %*% c(1, rev(pt[L + i - 1:l]))
}
nt.lm.pt = ts(pt, frequency = 12, start = c(1920, 1))
nt.lm.a = window(nottem, c(1921, 2)) - resid(nt.lm)
ts.plot(nottem, nt.lm.a, window(nt.lm.pt, c(1939, 12)), col = c("black",
"blue", "red"))Confrontiamo con la previsione di Holt-Winters.
nt.hw = HoltWinters(nottem)
nt.lm.pt = window(nt.lm.pt, c(1939, 12))
nt.hw.pt = predict(nt.hw, 24)
ts.plot(nottem, nt.lm.pt, nt.hw.pt, col = c("black", "red", "blue"))Proviamo a confrontare i modelli ottenuti con il metodo autoregressivo e con Holt-Winters. Analizziamo i residui.
layout(t(1:2))
# estrazione dei residui
nt.lm.r = resid(nt.lm)
nt.hw.r = resid(nt.hw)
# confronto grafico rispetto al tempo
plot(as.numeric(nt.hw.r), pch = 20)
plot(nt.lm.r, pch = 20)
# confronto grafico rispetto ai valori
plot(nt.hw$fitted[, 1], nt.hw.r, pch = 20)
plot(nt.lm.a, nt.lm.r, pch = 20)
# varianze spiegate
var(nt.hw.r)/var(window(nottem, 1921))
var(nt.lm.r)/var(window(nottem, 1921))
# acf
acf(nt.hw.r)
acf(nt.lm.r)
# pacf
pacf(nt.hw.r)
pacf(nt.lm.r)
# frequenze
hist(nt.hw.r, 20, freq = F)
lines(density(nt.hw.r), col = "blue")
lines(sort(nt.hw.r), dnorm(sort(nt.hw.r), mean(nt.hw.r), sd(nt.hw.r)), col = "red")
hist(nt.lm.r, 20, freq = F)
lines(density(nt.lm.r), col = "blue")
lines(sort(nt.lm.r), dnorm(sort(nt.lm.r), mean(nt.lm.r), sd(nt.lm.r)), col = "red")
# quantili
qqnorm(nt.hw.r, pch = 20)
qqline(nt.hw.r)
qqnorm(nt.lm.r, pch = 20)
qqline(nt.lm.r)
# test
shapiro.test(nt.hw.r)
shapiro.test(nt.lm.r)
layout(1)Eseguiamo una autovalidazione dei due metodi, usando i dati del 2010 come insieme test.
# training e test set
train = window(nottem, end = c(1938, 12))
test = window(nottem, c(1939, 1))
# previsione di Holt-Winters
train.hw = HoltWinters(train)
train.hw.pt = predict(train.hw, 12)
# previsione del modello autoregressivo
train.lm = lm(X14 ~ ., data = mnt)
pt = rep(0, L)
pt[1:(L - l + 1)] = nottem[1:(L - l + 1)]
for (i in 1:12) {
pt[L - l + 1 + i] = coef(train.lm) %*% c(1, rev(pt[L - l + 1 + i - 1:l]))
}
train.lm.pt = window(ts(pt, frequency = 12, start = c(1920, 1)), 1939)
# errore nelle predizioni
sqrt(mean((train.hw.pt - test)^2))
sqrt(mean((train.lm.pt - test)^2))
# confronto grafico
ts.plot(test, train.hw.pt, train.lm.pt, col = c("black", "red", "blue"))Consideriamo il dataset sugli autisti vittime di incidenti stradali in UK. Studiamo la stagionalità e le autodipendenze lineari. Emerge una dipendenza fino al lag \(14\).
Invece di fare la regressione a mano, usiamo il comando ar.
Rappresentiamo graficamente quanto ottenuto. Il calcolo delle previsioni viene dal modello \[ x_t - m = a_1(x_{t-1}-m)+\dots+a_p(x_{t-p}-m), \] ed \(m\) è la media dei dati. Per ottenere i valori stimati implementiamo la formula precedente.
uk.ar.an = uk
o = uk.ar$order
a = uk.ar$ar
for (i in (o + 1):length(uk)) {
uk.ar.an[i] = sum(rev(a) * uk[(i - o):(i - 1)]) + mean(uk) * (1 - sum(a))
}
uk.ar.an <- ts(uk.ar.an, frequency = 12, start = 1969)
ts.plot(uk, uk.ar.an, col = c("black", "red"))Il modello può essere ottenuto in maniera più semplice osservando che ar restituisce i residui.
Predizione con il modello regressivo.
Valutiamo l’incertezza (non parametrica) nella previsione.
uk.ar.r = uk.ar$resid[15:129]
y.max = max(uk.ar.pt + quantile(uk.ar.r, 0.975))
y.min = min(uk.ar.pt + quantile(uk.ar.r, 0.025))
ts.plot(uk.ar.pt, ylim = c(y.min, y.max))
lines(uk.ar.pt + quantile(uk.ar.r, 0.975), col = "red")
lines(uk.ar.pt + quantile(uk.ar.r, 0.025), col = "red")Il metodo fornisce anche le incertezze sulle predizioni, che ricordiamo però sono sensate solo nelle ipotesi di distribuzione Gaussiana e indipendenza dei residui
uk.ar.pt = predict(uk.ar, n.ahead = 12)
# non parametrico
ts.plot(uk.ar.pt$pred, ylim = c(y.min, y.max))
lines(uk.ar.pt$pred + quantile(uk.ar.r, 0.975), col = "red")
lines(uk.ar.pt$pred + quantile(uk.ar.r, 0.025), col = "red")
# parametrico
lines(uk.ar.pt$pred - uk.ar.pt$se, col = "green3")
lines(uk.ar.pt$pred + uk.ar.pt$se, col = "green3")Confrontiamo con Holt-Winters.
Analisi dei residui del modello regressivo.
Confrontiamo la predizione del metodo autoregressivo con la predizione di Holt-Winters.
Confrontiamo i due metodi mediante autovalutazione.
train = window(uk, end = c(1983, 12))
test = window(uk, 1984)
ukcv.hw.p = predict(HoltWinters(train), 12)
ukcv.ar.p = predict(ar(train), n.ahead = 12, se.fit = FALSE)
ts.plot(test, ukcv.hw.p, ukcv.ar.p, col = c("black", "red", "blue"))
sqrt(mean((test - ukcv.hw.p)^2))
sqrt(mean((test - ukcv.ar.p)^2))Una autovalutazione più robusta (su 11 mesi per semplicità di codice).
res.hw = rep(0, 11)
res.ar = rep(0, 11)
for (i in 1:11) {
train = window(uk, end = c(1984, i))
test = window(uk, start = c(1984, i + 1))
res.hw[i] = predict(HoltWinters(train), 1)
res.ar[i] = predict(ar(train), n.ahead = 1, se.fit = F)
}
test = window(uk, start = c(1984, 2))
sqrt(mean((test - res.hw)^2))
sqrt(mean((test - res.ar)^2))
ts.plot(test, res.hw, res.ar, col = c("black", "blue", "red"))Questo metodo considera anche serie non stazionarie. Esaminiamo la serie sull’andamento dei passeggeri di voli internazionali
Rappresentiamo graficamente quanto ottenuto. Ricordiamo che il calcolo delle previsioni viene dal modello \[ x_t - m = a_1(x_{t-1}-m)+\dots+a_p(x_{t-p}-m) + b, \] ed \(m\) è la media dei dati. L’implementazione diretta è la seguente.
o = ap.ls$order
a = ap.ls$ar
b = ap.ls$x.intercept
ap.ls.an = ap
for (i in (o + 1):length(ap)) {
ap.ls.an[i] = sum(rev(a) * ap[(i - o):(i - 1)]) + mean(ap) * (1 - sum(a)) +
b
}
plot(ap)
lines(ap.ls.an, col = "blue", lwd = 2)Avendo a disposizione i residui, il grafico precedente può essere ottenuto più rapidamente nel modo seguente.
Verifica dei residui.
layout(t(1:2))
ap.ls.r = as.double(na.omit(ap.ls$resid))
ap.ls.fitted = as.double(na.omit(ap - ap.ls$resid))
plot(ap.ls.r, pch = 20)
plot(ap.ls.fitted, ap.ls.r, pch = 20)
var(ap.ls.r)/var(ap.ls.r + ap.ls.fitted)
acf(ap.ls.r)
pacf(ap.ls.r)
hist(ap.ls.r, 20, freq = F)
lines(density(ap.ls.r), col = "red")
lines(sort(ap.ls.r), dnorm(sort(ap.ls.r), mean(ap.ls.r), sd(ap.ls.r)), col = "blue")
qqnorm(ap.ls.r, pch = 20)
qqline(ap.ls.r)
shapiro.test(ap.ls.r)
layout(1)Ricaviamo la previsione con il metodo appena calcolato. Corrediamo la previsione con una stima delle incertezze.
ap.ls.pt = predict(ap.ls, n.ahead = 12, se.fit = FALSE)
y.max = max(ap.ls.pt + quantile(ap.ls.r, 0.975))
y.min = min(ap.ls.pt + quantile(ap.ls.r, 0.025))
ts.plot(window(ap, 1959), window(ap - ap.ls$resid, 1959), ap.ls.pt, col = c("black",
"blue", "red"), lwd = c(1, 1, 2), ylim = c(y.min, y.max))
# stima empirica dell'incertezza
lines(ap.ls.pt + quantile(ap.ls.r, 0.975), col = "green3")
lines(ap.ls.pt + quantile(ap.ls.r, 0.025), col = "green3")
# stima parametrica dell'incertezza lines(ap.ls.pt +
# qnorm(0.975,mean(ap.ls.r),sd(ap.ls.r)), col = 'blue') lines(ap.ls.pt +
# qnorm(0.025,mean(ap.ls.r),sd(ap.ls.r)), col = 'blue')Confrontiamo analisi e previsioni con Holt-Winters.
# analisi
ap.hw = HoltWinters(ap, seasonal = "m")
ts.plot(ap, ap - ap.ls$resid, ap.hw$fitted[, 1], col = c("black", "blue",
"green3"))
# previsioni
ap.hw.pt = predict(ap.hw, 12)
ts.plot(ap.ls.pt, ap.hw.pt, col = c("blue", "green3"))
lines(ap.ls.pt + quantile(ap.ls.r, 0.975), col = "lightblue")
lines(ap.ls.pt + quantile(ap.ls.r, 0.025), col = "lightblue")Confrontiamo i due metodi con l’autovalutazione.
train = window(ap, end = c(1959, 12))
test = window(ap, 1960)
apcv.hw.p = predict(HoltWinters(train, seasonal = "m"), 12)
apcv.ls.p = predict(ar(train, method = "ols"), n.ahead = 12, se.fit = FALSE)
ts.plot(test, apcv.hw.p, apcv.ls.p, col = c("black", "red", "blue"))
sqrt(mean((apcv.hw.p - test)^2))
sqrt(mean((apcv.ls.p - test)^2))La funzione di cross-correlazione rivela le dipendenze tra traslazioni temporali di due serie storiche.
L’ordine delle due serie è importante
EuStockMarketsRiprendiamo il dataset degli indici delle borse europeee.
library(datasets)
dax = as.vector(EuStockMarkets[, "DAX"])
cac = as.vector(EuStockMarkets[, "CAC"])
smi = as.vector(EuStockMarkets[, "SMI"])
ftse = as.vector(EuStockMarkets[, "FTSE"])Proviamo la previsione usando ftse come fattore esogeno.
l = length(smi)
Op = smi[2:l]
Ip1 = smi[1:(l - 1)]
Ip2 = ftse[1:(l - 1)]
summary(lm(Op ~ Ip1 + Ip2))Esaminiamo le cross-correlazioni. Essendo dati economici, conviene esaminare i logaritmi.
Le serie hanno un trend comune (che spiega le cross-correlazioni con lag altissimi). Per eliminare questa influenza studiamo la cross-correlazione tra la serie delle differenze.
Il comando diff applicato a una serie storica \(x_1,x_2,\dots,x_n\) restituisce la nuova serie \(x_2-x_1,x_3-x_2,x_4-x_3,\dots,x_n-x_{n-1}\).
ccf(diff(log(ftse)), diff(log(dax)))
ccf(diff(log(ftse)), diff(log(cac)))
ccf(diff(log(ftse)), diff(log(smi)))
ccf(diff(log(smi)), diff(log(ftse)))Consideriamo ftse come fattore esogeno.
smi <- diff(log(smi))
ftse <- diff(log(ftse))
pacf(smi)
ccf(smi, ftse)
O = smi[2:l]
I1 = smi[1:(l - 1)]
I2 = ftse[1:(l - 1)]
summary(lm(O ~ I1 + I2))Il modello è terribile. Le due serie erano correlate dal solo trend comune.
SeatbeltsIncidenti stradali in Uk, dal 1/1969 al 12/1984:
DriversKilled: conducenti morti,drivers: conducenti morti o seriamente feriti,front: passeggeri anteriori morti o seriamente feriti,rear: passeggeri posteriori morti o seriamente feriti,kms: distanze percorse,PetrolPrice: prezzo carburante,VanKilled: conducenti di furgoni.dvr = Seatbelts[, 2]
kms = Seatbelts[, 5]
price = Seatbelts[, 6]
van = Seatbelts[, 7]
l = length(dvr)Consideriamo la serie dei conducenti. I lag 1,11,12,13,14 mostrano le più alte correlazioni.
in effetti,
Proviamo il modello autocorrelato a mano.
Op = dvr[15:l]
I1 = dvr[14:(l - 1)] # un mese
I2 = dvr[4:(l - 11)] # 11 mesi
I3 = dvr[3:(l - 12)] # 12 mesi
I4 = dvr[2:(l - 13)] # 13 mesi
I5 = dvr[1:(l - 14)] # 14 mesi
summary(lm(Op ~ I1 + I2 + I3 + I4 + I5))Sembra che si possano eliminare I4, I5, I2 senza perdere granché.
Consideriamo ora i fattori esogeni PetrolPrice, kms e van. Controlliamo le cross-correlazioni. C’è un massimo a 10 mesi tra drivers e kms.
Troppe correlazioni tra drivers e price.
Correlazioni complesse anche tra drivers e van
Calcoliamo la varianza spiegata del modello di base precedente, e poi la varianza spiegata dei tre modelli con il fattore esogeno kms, il fattore esogeno PetrolPrice e il fattore esogeno van. Il fattore kms fornisce un maggiore guadagno di varianza spiegata.
dvr.R = summary(lm(Op ~ I1 + I3))$r.squared
Ek = kms[5:(l - 10)] # kms a 10 mesi
Ep = price[1:(l - 14)] # price a 14 mesi
Ev = van[3:(l - 12)] # van a 12 mesi
summary(lm(Op ~ I1 + I3 + Ek))$r.squared - dvr.R
summary(lm(Op ~ I1 + I3 + Ep))$r.squared - dvr.R
summary(lm(Op ~ I1 + I3 + Ev))$r.squared - dvr.RGenerare due serie storiche, la prima più adatta per l’analisi con il metodo di Yule-Walker, la seconda più adatta per l’analisi con il metodo dei minimi quadrati. Implementare poi l’analisi delle due serie e ottenere per ognuna una previsione, corredata della misura di incertezza.
Siano date due serie storiche \(S_0\), \(S_e\) di \(170\) dati con periodicità mensile. Formulare e implementare un modello autoregressivo per la previsione della serie \(S_0\) che preveda l’uso dei dati del mese precedente e dello stesso mese dell’anno precedente di \(S_0\), e del dato del semestre precedente della variabile esogena \(S_e\).
Analizzare i residui dei modelli autoregressivo e Holt-Winters per il dataset UKDriverDeaths.
Svolgere una analisi per il confronto tra il metodo Yule-Walker e il metodo Holt-Winters per il dataset nottem, allo scopo di identificare il modello ritenuto più adeguato. Valutare l’incertezza sulla previsione per il modello identificato.
Svolgere una analisi per il confronto tra il metodo autoregressivo dei minimi quadrati, il metodo Holt-Winters, ed un opportuno modello autoregressivo implementato direttamente, per il dataset AirPassengers, allo scopo di identificare il modello ritenuto migliore. Valutare l’incertezza sulla previsione per il modello identificato.
Determinare, attraverso l’analisi dei residui e l’autovalidazione, del modello più adeguato per il dataset co2. Ricavare poi una predizione corredata di stime di incertezza.
Determinare, attraverso l’analisi dei residui e l’autovalidazione, del modello più adeguato per il dataset JohnsonJohnson. Ricavare poi una predizione corredata di stime di incertezza.
Determinare, attraverso l’analisi dei residui e l’autovalidazione, del modello più adeguato per il dataset Nile. Ricavare poi una predizione corredata di stime di incertezza.
Determinare, attraverso l’analisi dei residui e l’autovalidazione, del modello più adeguato per il dataset UKgas. Ricavare poi una predizione corredata di stime di incertezza.