Comandi e funzioni di R per le Serie storiche
di Michele Faggion
Funzioni
Descrizione
plot(y,xlab=” “,ylab=” “)
Disegna un grafico.
y(serie) =window(se- Estraiamo dalla serie una dei dati, ovvero facciamo partire la
rie,start=X, end=Y)
serie da una determinata osservazione o la facciamo finire
ad unʼaltra determinata (end).
230
200
BJsales
260
abline(h=mean(serie),- Disegna nel grafico una linea che rappresenta la media della
lty=”dashed”)
serie di dati.
0
50
100
150
t
deltat(serie)
Distanza fra le osservazioni, è il reciproco della frequenza.
frequency(serie)
Frequenza
0
-4 -2
y
2
4
lag.plot(y,1)
Grafico di autodispersione
lag.plot(y,set.lags=c(1, Mi dice se cʼè un legame (di qualsiasi tipo, lineare, esponen17,18))
ziale, ecc.) fra le coppie di valori ad un determinato ritardo
(nellʼesempio, ritardo 1).
Se la dipendenza per ritardi superiore alla distanza fra le osservazioni non si osserva, allora significa che il processo
osservato ha una memoria “corta”: osservazioni vicine
sono fra loro correlate, ma più lontane non lo sono più.
-10
-5
0
5
lag 1
acf(y,20)
1.0
0.0
-1.0
Distanza nel tempo
Autocorrelogramma
0
50
100
Autocorrelazione campionaria
150
Autocorrelazione (funzione di) - Autocorrelogramma
I coefficienti di autocorrelazione sono la stima della correlazione che si osserva nei diagrammi di autodispersione.
Il secondo termine (20) indica il numero di coefficienti di correlazione che vogliamo stimare, ovvero il numero di ritardi.
Lʼasse delle ascisse riporta la distanza nel tempo. Ogni semento rappresenta la correlazione fra due osservazioni distanti x nel tempo.
a) se il segmento ad una certa distanza esce dalla banda,
allora cʼè correlazione, cʼè dipendenza fra le osservazioni a quella distanza.
b) quando le correlazioni (ovvero i segmenti) sono dentro le
bande, allora significa che sono = 0.
Si ricorda che per ritardi maggiori, gli errori sono maggiori.
Funzioni
Descrizione
Box.test(y,10)
Test Ljung-Box - Stima della significatività dei coefficienti di autocorrelazione
Il secondo termine sono il numero di coeff. di autocorrelazione che vogliamo stimare.
Lʼipotesi nulla è H:coeffi = 0
Quando il p-value è > 0.05 accettiamo lʼipotesi nulla per la
quale non cʼè autocorrelazione fra le osservazioni (serie white noise).
Se p-value < 0.05 allora cʼè autocorreazione (rifiuto lʼipotesi
nulla), si nota una struttura di dipendenza dei dati. Quindi,
si dice, esiste una correlazione seriale fra le osservazioni.
Box-Pierce test
data: y
X-squared = 224.0319
df = 10
p-value < 2.2e-16
seaplot(serie)
Disegna i grafici delle sotto serie mensili (tutte le osservazioni di gennaio, febbraio, ecc.). Lʼordine è da sx a dx, dal
basso allʼalto. Il grafico mostra lʼandamento stagionale della
serie, ovvero mostra, ad esempio, se le osservazioni di gennaio sono più basse o più alte do quelle di marzo, ecc.
monthplot(serie)
Grafico che mostra la serie in modo diverso per ciascun mede. Sono evidenziate le medie aritmetiche delle serie mensili attraverso le linee orizzontali.
esId(serie)
Stima del modello
La funzione attraverso la tecnica della massima verosomiglianza, stima tutti i modelli basati sul lisciamento esponenziale compatibili con la serie osservata e mostra i modelli osservati secondo il criterio BIC.
Nellʼoutput vengono indicati il tipo di deriva, stagionalità e innovazione, il numero di parametri (basata sulla valutazione
della parsimonia del modello, se è più parsimonioso, rischia
meno di cogliere caratteristiche spurie della serie osservata)
e la log-verosomiglianza (ovvero la somma dei quadrati dei
residui usata dai criteri per valutare la bontà dʼadattamento
del modello).
Funzioni
Descrizione
esFit(serie, drift=” “, Stima dei parametri del modello
seasonality= “ “, in- La funzione mi restituisce la stima dei parametri del modello
novation= “ “)
(non sono i parametri veri), che dobbiamo interpretare:
- alfa:
- phi:
- beta:
- l.start:
- d.start:
- sigma:
- beta:
tsdiag(modellostima- Verifica dellʼadattamento
to, h)
La funzione restituisce tre grafici:
- i residui standardizzati (stima);
- lʼautocorrelazione campionaria dei residui (lʼautocorrelogramma): se il modello è buono, allora le autocorrelazioni
devono rimanere allʼinterno delle bande e quindi non sono
correlate.
- il loro livello di significatività tramite il test Ljung-Box calcolato sul numero h dei coefficienti di correlazione.
residuals(modellostima- Calcola la stima dei residui del modello.
to)
acf(residuals(modello- Ci fornisce uno “zoom” dellʼautocorrelogramma dei residui (il
stimato))
secondo grafico della funzione tsdiag).
Box.test(residuals(mod La funzione mi fornisce unʼulteriore stima dei p-value sotto
ellostimato),type=”L”)
lʼipotesi nulla di una serie white noise (ovvero che lʼautocorrealazione fra osservazioni sia = 0). Se > 0.05, allora accetto
lʼipotesi nulla e posso procedere con la previsione.
yc=esFilter(m)
La funzione determina lʼevoluzione stimata del livello, della
deriva, della stagionalità e dellʼinnovazione in tutti gli istanti di
tempo (sulle componenti stimate)
plot(yc)
Posso disegnare il grafico dellʼevoluzione delle componenti
stimate della serie.
Funzioni
Descrizione
yp<-predict(modellostimato,passo in
avanti, method =
“gauss”)
Previsione
Con questo comando, R calcola le previsioni per un numero
x di passi in avanti (secondo argomento), condizionatamente
al modello stimato (primo argomento), ovvero condizionato al
passato, utilizzando il metodo gaussiano, ovvero assumendo
che le innovazioni abbino distribuzioni normali.
Il comando restituisce previsioni puntuali e intervallari.
plot(yp[,c("5%", "mean", I comandi servono a disegnare il grafico di previsione delle
"95%")], plot.type = "s", serie previste, mentre con points sovrascrivo i valori osservalty = "dotdash")
ti.
points(BJsales, pch =
"*", cex = 2)
Migliorare le previsioni già esistenti
w=y-yp
Calcolo i residui fra le osservazioni vere e quelle previste dallʼente esterno.
g= y - residuals(m)
eqm1 = mean (w^2)
Errore quadratico medio
Lʼinnovazione è addittiva
eqm2 = mean (y-g)^2)
Errore quadratico medio
Lʼinnovazione è moltiplicativa
h= w - residuals (m.w)
Stima dellʼerrore quadratico medio per il terzo modello (m.w)
eqm3 = mean (y - che elaboriamo per migliorare le previsioni.
(yp+h)^2)