Proviamo il metodo di smorzamento esponenziale su una serie con forte trend creata artificiosamente.
st = ts((1:100) + 10 * rnorm(100), frequency = 10)
se = HoltWinters(st, beta = F, gamma = F)
summary(se)
se
plot(se)
# provare con valore di alpha più opportuno
# se=HoltWinters(st,alpha=?,beta=F,gamma=F)
predict(se, 1)
plot(se, predict(se, 1))Passiamo a esaminare la serie esempio sulle anomalie di temperatura, che contiene dati con trend dominante.
tp.data <- read.csv("13anomalie.csv")
tp = ts(tp.data, frequency = 1, start = 1880)
tp.se = HoltWinters(tp, beta = F, gamma = F)
tp.se$alpha
plot(tp.se)Proviamo a variare il valore del parametro \(\alpha\in\{1,0.9,0.7,0.6,0.5,0.4,0.3,0.2,0.1\}\)
plot(tp.se)
points(HoltWinters(tp, alpha = 1, beta = F, gamma = F)$fitted, col = "blue",
type = "l")Non sembra essere utile variare la condizione iniziale.
plot(tp.se)
points(HoltWinters(tp, alpha = tp.se$alpha, beta = F, gamma = F, l.start = -0.4)$fitted,
col = "blue", type = "l")Predizione con smorzamento esponenziale.
Includiamo il trend nella analisi.
Proviamo a variare direttamente i valori dei parametri. Prima \(\alpha\in\{0.9,0.8,0.7,0.6,0.5\}\),
plot(tp.set)
points(HoltWinters(tp, alpha = 0.9, beta = tp.set$beta, gamma = F)$fitted,
col = "blue", type = "l")poi \(\beta\in\{0.1,0.2,0.3,0.4,0.5\}\),
plot(tp.set)
points(HoltWinters(tp, alpha = tp.set$alpha, beta = 0.1, gamma = F)$fitted,
col = "blue", type = "l")Riaggiustiamo la condizione iniziale. Ricordiamo che l.start è la intercetta iniziale, e b.start è la pendenza iniziale.
plot(HoltWinters(tp, gamma = F, l.start = -0.1))
plot(HoltWinters(tp, gamma = F, l.start = 0, b.start = 0.01))
plot(HoltWinters(tp, gamma = F, l.start = 0.05, b.start = -0.02))Proviamo una ``mini’’ regressione lineare.
x = 1:20
coefficients(lm(tp[1:20] ~ x))
plot(HoltWinters(tp, gamma = F, l.start = -0.08, b.start = -0.005))Abbiamo ancora la possibilità di calibrare \(\alpha\) e \(\beta\).
tp.set = HoltWinters(tp, gamma = F, l.start = -0.08, b.start = -0.005)
c(tp.set$alpha, tp.set$beta)
plot(HoltWinters(tp, gamma = F, l.start = -0.08, b.start = -0.005, alpha = 0.2,
beta = 0.1))Con tali valori, osserviamo tutte le componenti (in verde la serie iniziale, in blu la previsione, in nero le intercette e in rosso i coefficienti angolari).
tp.set = HoltWinters(tp, alpha = 0.2, beta = 0.1, gamma = F, l.start = -0.08,
b.start = -0.005)
ts.plot(tp, tp.set$fitted, col = c("green3", "blue", "black", "red"))Quale metodo prevede meglio? Confrontiamo le predizioni dei metodi di smorzamento esponenziale senza/con trend.
tp.end = window(tp, 2015)
tp.begin = window(tp, 1880, 2014)
tp.se = HoltWinters(tp.begin, beta = F, gamma = F)
tp.se.pt = predict(tp.se, 4)
tp.set = HoltWinters(tp.begin, alpha = 0.2, beta = 0.1, gamma = F, l.start = -0.08,
b.start = -0.005)
tp.set.pt = predict(tp.set, 4)
ts.plot(tp.end, tp.se.pt, tp.set.pt, col = c("black", "red", "blue"), type = "b")Realizziamo un confronto più strutturato mediante autovalutazione.
m = 20
l = length(tp)
res.se = rep(0, m)
res.set = rep(0, m)
j = 1
for (i in (l - m):(l - 1)) {
tp_sub = ts(tp.data[1:i, 1], frequency = 1, start = 1880)
tp_sub.se = HoltWinters(tp_sub, beta = F, gamma = F)
tp_sub.set = HoltWinters(tp_sub, alpha = 0.2, beta = 0.1, gamma = F,
l.start = -0.08, b.start = -0.005)
res.se[j] = tp.data[i + 1, 1] - predict(tp_sub.se, n.ahead = 1, se.fit = F)
res.set[j] = tp.data[i + 1, 1] - predict(tp_sub.set, n.ahead = 1, se.fit = F)
j = j + 1
}
plot(res.se, type = "b", col = "blue", pch = 20)
points(res.set, type = "b", col = "green3", pch = 20)
sqrt(mean(res.se^2))
sd(res.se)
sqrt(mean(res.se^2))
sd(res.set)Possiamo applicare i due metodi di smorzamento per la previsione del trend della serie UKDriverDeaths.
uk = UKDriverDeaths
uk_t = stl(uk, 7)$time.series[, 2]
uk_t.se = HoltWinters(uk_t, beta = F, gamma = F)
uk_t.se$alpha
uk_t.set = HoltWinters(uk_t, gamma = F)
c(uk_t.set$alpha, uk_t.set$beta)
ts.plot(uk_t.se$fitted[, 1], uk_t.set$fitted[, 1], col = c("blue", "green3"))
lines(uk_t, col = "red")Confrontiamo le previsioni dal trend della decomposizione, rispetto ad una previsione diretta.
Passiamo ad esaminare serie con carattere stagionale. Non ci aspettiamo in generale risultati soddisfacenti.
SE non è soddisfacente, segue la stagionalità (e la scelta della condizione iniziale non aiuta).
st = ts(sin(1:100) + 0.2 * rnorm(100), frequency = 10)
se = HoltWinters(st, beta = F, gamma = F)
plot(se)Anche SET segue la stagionalità, ma un po’ meno, almeno sui valori ottimali standard.
Naturalmente, poiché non c’è trend, una scelta opportuna delle condizioni iniziali cattura (inutilmente, non è lo strumento giusto) il trend al netto della stagionalità
Esaminiamo la serie delle temperature dell’aria di Nottingham (serie con carattere stagionale). Sia SE che SET sono fedeli alla serie senza catturare il trend.
nt.se = HoltWinters(nottem, beta = F, gamma = F)
plot(nt.se)
nt.set = HoltWinters(nottem, gamma = F)
plot(nt.set)Per catturare il trend, proviamo ad aggiustare la condizione iniziale.
Esaminiamo la serie sulle temperature includendo la stagionalità.
Il comando HoltWinters restituisce tra l’altro la sotto-struttura fitted, dove in particolare level sono le intercette e xhat sono i valori stimati. Può risultare utile graficamente per calibrare il modello.
Esploriamo i risultati confrontandoli con quelli della decomposizione
nt.stl = stl(nottem, "periodic")
# confronto grafico (delle intercette) con il trend di stl
ts.plot(nt.stl$time.series[, 2], nt.hw$fitted[, 2], col = c("black", "red"))
# confronto grafico con la stagionalità di stl
ts.plot(nt.stl$time.series[, 1], nt.hw$fitted[, 4], col = c("black", "red"))Proviamo a ottenere una previsione.
layout(t(1:2))
plot(nt.hw, predict(nt.hw, 12), main = "Previsione a 12 mesi")
plot(nt.hw, predict(nt.hw, 24), main = "Previsione a 24 mesi")
layout(1)Estraiamo i residui.
Calcoliamo l’incertezza per via non parametrica.
plot(nt.hw, predict(nt.hw, 12))
lines(predict(nt.hw, 12) + quantile(nt.hw.r, 0.05), col = "green3")
lines(predict(nt.hw, 12) + quantile(nt.hw.r, 0.95), col = "green3")Esaminiamo brevemente i residui del metodo HW.
plot(nt.hw.r, type = "p", pch = 20)
plot(nt.hw$fitted[, 1], nt.hw.r, pch = 20)
# Attenzione! La varianza spiegata va calcolata su una finestra comune.
start(nottem)
end(nottem)
start(nt.hw.r)
end(nt.hw.r)
var(nt.hw.r)/var(window(nottem, c(1921, 1)))
acf(nt.hw.r)
hist(nt.hw.r, 20, freq = F)
lines(density(nt.hw.r))
lines(sort(nt.hw.r), dnorm(sort(nt.hw.r), mean(nt.hw.r), sd(nt.hw.r)), col = "red")
qqnorm(nt.hw.r, pch = 20)
qqline(nt.hw.r)
shapiro.test(nt.hw.r)Calcoliamo l’incertezza per via parametrica.
plot(nt.hw, predict(nt.hw, 12))
lines(predict(nt.hw, 12) + qnorm(0.05, mean(nt.hw.r), sd(nt.hw.r)), col = "blue")
lines(predict(nt.hw, 12) + qnorm(0.95, mean(nt.hw.r), sd(nt.hw.r)), col = "blue")Visualizziamo le incertezze al 90% sovrapposte.
plot(nt.hw, predict(nt.hw, 12))
lines(predict(nt.hw, 12) + qnorm(0.05, mean(nt.hw.r), sd(nt.hw.r)), col = "blue")
lines(predict(nt.hw, 12) + qnorm(0.95, mean(nt.hw.r), sd(nt.hw.r)), col = "blue")
lines(predict(nt.hw, 12) + quantile(nt.hw.r, 0.05), col = "green3")
lines(predict(nt.hw, 12) + quantile(nt.hw.r, 0.95), col = "green3")Una nota a margine: il comando predict per Holt-Winters prevede l’opzione prediction.interval. L’interpretazione del risultato è sensata solo se si sa cosa significhi e come sia stato ottenuto. Non ha senso usarlo altrimenti. Ad esempio,
ts.plot(predict(nt.hw, 12, prediction.interval = T, level = 0.9))
lines(predict(nt.hw, 12) + qnorm(0.05, mean(nt.hw.r), sd(nt.hw.r)), col = "blue")
lines(predict(nt.hw, 12) + qnorm(0.95, mean(nt.hw.r), sd(nt.hw.r)), col = "blue")
lines(predict(nt.hw, 12) + quantile(nt.hw.r, 0.05), col = "green3")
lines(predict(nt.hw, 12) + quantile(nt.hw.r, 0.95), col = "green3")Testiamo la capacità previsiva del metodo di smorzamento esponenziale con trend e stagionalità. In mancanza di dati futuri noti (set di validazione) con i quali confrontare le previsioni, (eventualmente in caso di confronto tra modelli) possiamo utilizzare una procedura di autovalidazione. Tale procedura deve essere rispettosa della progressione temporale.
train = window(nottem, end = c(1938, 12))
test = window(nottem, 1939)
nt.hw <- HoltWinters(train)
nt.hw.p = predict(nt.hw, 12)
sqrt(mean((nt.hw.p - test)^2))Esaminiamo le differenze graficamente.
Consideriamo la serie sull’andamento dei passeggeri di voli internazionali tra il 1949 e il 1960. Valutiamo la differenza tra l’analisi e la predizione del metodo di Holt-Winter con modello stagionale additivo e moltiplicativo
ap = AirPassengers
ap.hwa = HoltWinters(ap, seasonal = "additive")
ap.hwm = HoltWinters(ap, seasonal = "multiplicative")
ts.plot(ap, ap.hwa$fitted[, 1], ap.hwm$fitted[, 1], col = c("black", "blue",
"red"))Valutiamo i due modelli prima in termini di residui
layout(t(1:2))
# estrazione dei residui
ap.hwa.r = resid(ap.hwa)
ap.hwm.r = resid(ap.hwm)
# proporzione di varianza non spiegata
var(ap.hwa.r)/var(window(ap, 1950))
var(ap.hwm.r)/var(window(ap, 1950))
# rappresentazione grafica rispetto al tempo
plot(ap.hwa.r, type = "p", pch = 20)
plot(ap.hwm.r, type = "p", pch = 20)
# rappresentazione grafica rispetto ai valori stimati
plot(as.numeric(ap.hwa$fitted[, 1]), as.numeric(ap.hwa.r), type = "p", pch = 20)
plot(as.numeric(ap.hwm$fitted[, 1]), as.numeric(ap.hwm.r), type = "p", pch = 20)
# autocorrelazione
acf(ap.hwa.r)
acf(ap.hwm.r)
# densità empiriche
hist(ap.hwa.r, 20, freq = F)
lines(density(ap.hwa.r), col = "blue")
lines(sort(ap.hwa.r), dnorm(sort(ap.hwa.r), mean(ap.hwa.r), sd(ap.hwa.r)),
col = "red")
hist(ap.hwm.r, 20, freq = F)
lines(density(ap.hwm.r), col = "blue")
lines(sort(ap.hwm.r), dnorm(sort(ap.hwm.r), mean(ap.hwm.r), sd(ap.hwm.r)),
col = "red")
# grafico quantile-quantile
qqnorm(ap.hwa.r, pch = 20)
qqline(ap.hwa.r)
qqnorm(ap.hwm.r, pch = 20)
qqline(ap.hwm.r)
# test
shapiro.test(ap.hwa.r)
shapiro.test(ap.hwm.r)
layout(1)Valutiamo poi i due modelli in termini di capacità di predizione.
ap_train = window(ap, end = c(1959, 12))
ap_test = window(ap, 1960)
ap.hwa = HoltWinters(ap_train, seasonal = "additive")
ap.hwm = HoltWinters(ap_train, seasonal = "multiplicative")
ap.hwa.p = predict(ap.hwa, 12)
ap.hwm.p = predict(ap.hwm, 12)
sqrt(mean((ap.hwa.p - ap_test)^2))
sqrt(mean((ap.hwm.p - ap_test)^2))Osserviamo il risultato graficamente.
Proviamo una autovalidazione più robusta.
l = length(ap)
res.hwa = rep(0, 24)
res.hwm = rep(0, 24)
j = 1
for (i in (l - 24):(l - 1)) {
ap_cv = ts(ap[1:i], frequency = 12, start = c(1949, 1))
ap.hwa = HoltWinters(ap_cv, seasonal = "additive")
ap.hwm = HoltWinters(ap_cv, seasonal = "multiplicative")
ap.hwa.p = predict(ap.hwa, 1)
ap.hwm.p = predict(ap.hwm, 1)
res.hwa[j] = ap.hwa.p - ap[i + 1]
res.hwm[j] = ap.hwm.p - ap[i + 1]
j = j + 1
}
sqrt(mean(res.hwa^2))
sqrt(mean(res.hwm^2))
plot(res.hwa, type = "b", pch = 20, col = "blue")
lines(res.hwm, type = "b", pch = 20, col = "green3")Generare una serie storica \(S\) dominata principalmente da trend, e perturbata da fluttuazioni casuali. Implementare una stima della inizializzazione della intercetta per lo smorzamento esponenziale con trend.
Generare una serie storica caratterizzata principalmente da trend e rumore, e tale che la previsione ottenuta dal metodo di smorzamento esponenziale con trend non sia coerente con l’aspettativa creatasi con la intuizione dello sperimentatore.
Generare una serie storica caratterizzata da una netta stagionalità e un rumore non trascurabile. Analizzare la serieper mezzo dello smorzamento esponenziale con trend. Determinare la scelta più opportuna dei parametri per migliorare l’analisi.
Data la serie storica
ottenere una previsione del primo periodo futuro per mezzo di un metodo opportunamente scelto. Tracciare il grafico della serie arricchito della previsione e di bande di confidenza al 95%.
Ripetere l’analisi del trend estratto da stl. Calibrare il valore di \(\alpha\) per il metodo SE e i valori di \(\alpha,\beta\) e delle condizioni iniziali per il metodo SET, trovando eventuali valori più convincenti.
Confrontare l’ampiezza dell’incertezza valutata per via parametrica attraverso quantili empirici (confidenza al 95%) con l’ampiezza del prediction.interval fornito dal comando predict.
Confrontare i residui ottenuti dal metodo di smorzamento esponenziale con trend e dal metodo di smorzamento esponenziale con trend e stagionalità ottenuti sul dataset delle anomalie di temperatura.
Considerare il dataset nottem e confrontare il modello di smorzamento esponenziale con trend e il modello di smorzamento esponenziale con trend e stagionalità
Confrontare le previsioni fornite dallo smorzamento esponenziale e dallo smorzamento esponenziale con trend sul trend ottenuto attraverso la decomposizione per il dataset nottem.
Considerare il dataset UKDriverDeaths. Selezionare il metodo più opportuno per analizzare i dati e calibrare il modello scelto. In corrispondenza dei valori dei parametri scelti analizzare i residui del modello e testarne le capacità predittive, includendo una misura di incertezza delle previsioni. Confrontare con i dati effettivi.