R, l’analisi delle serie storiche partendo da Copenaghen

Eccoci di nuovo a parlare di R, linguaggio ed un ambiente open source per il calcolo statistico e la rappresentazione grafica in ambito statistico.

Anche questa volta faccio frutto del “dovere” per scrivere questo articolo per “piacere”.

Nel titolo leggete Copenaghen, mi spiego meglio.

Questa volta parliamo di come effettuare analisi delle serie storiche con il software R.

Per analizzare una serie storica dobbiamo averne una. Io che ogni tanto sbarco sul sito del settimanale internazionale ho memoria di una serie storica il cui grafico è apparso in un post sul sito di internazionale.

Ho cercato l’articolo e l’ho ritrovato: http://www.internazionale.it/home/?p=12879 .

(Potremmo anche riportare l’intero articolo dato che Internazionale fa uso di Creative Commons Attribution-Non-Commercial-Share Alike 2.5.)

Bene, ho pensato, perchè mai cercare i dati delle vendite di una gelateria alla ricerca di trend e pattern stagionali quando possiamo cercare serie storiche sulla concentrazione della co2 visto che siamo nel dopo Copenaghen ed il problema del riscaldamento terrestre e più che mai di interesse?

Ho cercato quindi tramite i motori di ricerca un documento che contenesse i dati sui rilevamenti di co2 di diversi anni campionati mensilmente ed ho trovato questo documento: http://cdiac.ornl.gov/ftp/trends/co2/maunaloa.co2 .

Proviamo quindi a studiare con R i dati relativi alle concentrazioni mensili di Co2 dal 2000 al 2008.

Prima di poter lavorare sui dati dobbiamo importarli in R e per farlo copiamo i valori mensili in un’unica riga in un foglio di calcolo; io ho usato Calc.

Una volta trascritti i dati in un’unica riga, salviamo il foglio in formato csv e scegliamo, al momento di salvare, la Tabulazione come delimitatore di campo.

Adesso possiamo importare i dati in R e iniziarne lo studio.

Con il seguente comando importiamo i dati in R e assegniamo ad un oggetto data frame i valori della serie storica.

>A=read.table(file=”/home/andrea/Documents/Flandoli/homework2/SeriStoricaCo2.csv”)

Per visualizzare una serie storica abbiamo bisogno di avere un vettore da passare come argomento al comando ts.plot().

Quindi digitiamo:

>X=t(A[1,])

cosi da trasporre il vettore.

X ora conterrà i 108 valori corrispondenti ai valori mensili di co2 dal 2000 al 2008.

Ora usiamo il comando ts.plot.

>ts.plot(X)

A questo punto ci apparirà il grafico della serie storica.

ts.plot(X) - X contiene la serie storica delle concentrazioni di co2 dal 2000 al 2008

Già da una analisi visiva si desume un trend per cui i valori crescono quasi linearmente e una periodicità stagionale con un periodo vicino a 12 mesi.

Con R possiamo produrre il grafico della funzione di autocorrelazione della serie storica con il comando acf(X), ottenendo il grafico della figura seguente.

>acf(X)

acf(X) - grafico della funzione di autocorrelazione con R

Come è giusto che sia, in zero ovvero con una traslazione di 0, la correlazione è pari a 1.

Il grafico della funzione di autocorrelazione mostra che vi è una alta correlazione tra la serie storica e la sua traslata di 11 e 12 unità e questo ci è utile per impostare poi il modello di previsione lineare. I picchi sono oltre l’intervallo di confidenza quindi possiamo prenderli in considerazione.

L’impressione sulla periodicità della serie storica non era errata ma tutto sommato era evidente, ed in questo modo ne abbiamo avuto conferma tramite la funzione di autocorrelazione calcolata da R.

Se assegniamo i risultati del comando acf(X) ad una variabile Y con il comando

>Y<-acf(X)

oltre ad avere il grafico della funzione di autocorrelazione avremo in Y i valori della funzione di autocorrelazione che appunto risultano corrispondenti al grafico e valere 0.711 per un periodo di 11 mesi e 0.706 per un periodo di 12 mesi.

Grafico della funzione di autocorrelazione in R

Avendo 108 valori nella serie storica e volendo visualizzare la presenza di un trend lineare, possiamo costruire un vettore con gli indici temporali ed effettuando la regressione lineare visualizzare intercetta e coefficiente così da poter costruire una retta da sovrapporre al grafico della serie storica che evidenzi il trend.

>Z<-c(1:108)

>reg<-lm(X~Z) // con questo comando R applica il modello di regressione lineare tra la variabile Z che contiene i valori temporali e X che contiene la serie storica ed in più associamo il risultato prodotto da R ad un oggetto reg.

>cor(X,Z) // con il comando cor è possibile calcolare la correlazione tra due variabili
>summary(reg) // con questo comando visualizziamo una serie di informazioni risultati dal calcolo della regressione lineare effettuato da R
F<-(0:1080)/10
b=reg$coefficients[1] //associamo il valore dell’intercetta ad un oggetto b
a=reg$coefficients[2]
G=a*F+b //costruiamo la retta con coefficiente angolare a
ts.plot(X) //
lines(F,G) // il comando lines permette di sovrapporre un grafico ad un altro grafico precedentemente plottato difatti se abbiamo chiuso la finestra del grafico precedentemente plottato il comano restituirà un errore.

Con la sequenza di comandi appena scritti otteniamo il seguente risultato:

Serie di comandi per sovrapporre un retta ottenuta tramite la regressione lineare al grafico della serie storica

Ora possiamo costruire un modello che descriva i dati della serie storica.

Ho scelto di costruire e confrontare tre modelli lineari che descrivano i dati visto anche il grafico della funzione di autocorrelazione il quale mostra una periodicità alta tra un periodo di 11 e 12.

1)xk=a1x(k-1)+a12x(k-12)+b+ek

2)xk=a1x(k-1)+a11x(k-11)+b+ek

3)xk=a1x(k-1)+a11x(k-11)+a12x(k-12)+b+ek

modello2 dei dati

modello 3 dei dati

La varianza spiegata nel terzo modello, che prende in considerazione una periodicità 11 e 12 insieme è maggiore.

Dai p-value molto vicini al valore zero si deduce che le ipotesi che le variabili non siano utili come predittori sono quasi nulle.

Possiamo usare quindi i valori a1, a2, a3 e b presi dalla regressione lineare del terzo modello per effettuare le previsioni poiché il terzo modello spiega maggiore varianza dei dati x0 da 13 a n e “fidandomi” anche del grafico della funzione di autocorrelazione.

Previsione dei valori futuri

Associando ad un vettore P i valori della serie storica presenti nel vettore X con un ciclo for possiamo calcolare i valori dei mesi del 2009 utilizzando come modello di previsione lineare quello costruito con i coefficienti e l’intercetta del terzo modello di descrizione dei dati della serie storica che ho scelto.

>x01<-X[13:n]
>x11<-X[12:(n-1)]
>x21<-X[1:(n-12)]
>REG1<-lm(x01~x11+x21)
>summary(REG1)
>x02<-X[13:n]
>x12<-X[12:(n-1)]
>x22<-X[2:(n-11)]
>REG2<-lm(x02~x12+x22)
>summary(REG2)

>x03<-X[13:n]
>x13<-X[12:(n-1)]
>x23<-X[2:(n-11)]

>x33<-X[1:(n-12)]
>REG3<-lm(x03~x13+x23+x33)
>summary(REG3)

>P<-X
>length(P)
>a11<-REG3$coefficients[2]
>a21<-REG3$coefficients[3]
>a31<-REG3$coefficients[4]
>b<-REG3$coefficients[1]
>for (k in (n+1) : (n+12) ) {
P[k]=a11*P[k-1]+a21*P[k-12]+a31*P[k-11]+b
}
>length(P) // ora p dovrebbe risultarvi composto da 120 valori, i primi 108 della serie + dodici dalla previsione per

// il 2009 basata sul modello lieneare

ts.plot(P,col=”red”)
lines(X,col=”blue”)

Dal grafico seguente si distingue la previsione per l’anno 2009 di valori di co2 che è di colore rosso.

Previsione con modello lineare

Tra le  mlle potenzialità di R c’è la possibilità di calcolare le previsioni con il metodo di Holt-Winters.

Possiamo usare i comandi di R per predire l’anno 2009 e generare un grafico che sovrapponga la previsione con il metodo di Holt-Winters (in giallo) con quella calcolata con il modello lineare ottenendo il risultato visibile nella figura seguente.

>x.new<-ts(data=X,frequency=12) // il comando HoltWinters ha bisogno che all’oggetto contenente la serie storica sia  //associato il parametro che ne indica la periodicità quindi in realtà stiamo confrontando il modello lineare con quello di //Holt-Winters considerando una frequenza pari a 12
>HW<-HoltWinters(x.new)
>predHW<-predict(HW,12)
>C<-X
>for (k in 1 : 12 ) {
C[n+k]=predHW[k]
}
>ts.plot(P,col=”red”)
>lines(C,col=”yellow”)
>lines(X,col=”blue”)

http://farm3.static.flickr.com/2790/4277118593_88e58e4012_o.png

Si nota bene già graficamente che le due previsioni per quanto simili in grandi linee non sono le stesse.

Training Set e Set Test per il confronto tra il modello lineare e quello di Holt-Winters

Avendo a disposizione dati per molto anni possiamo prendere i valori della serie storica fino al 2007 e calcolare le previsioni per il 2008 con il terzo modello lineare di cui prima e con il metodo Holt-Winters per metterle a confronto con i dati veri della serie storica del 2008.

Quindi partendo da n=n1+n2 dati della serie ho usiamo n1=96 dati per training set e n2=12 come test set.

>X1<-X[1:96]
>n1<-length(X1)
>P1<-X1
>for (k in (n1+1) : (n1+12) ) {
P1[k]=a11*P1[k-1]+a21*P1[k-12]+a31*P1[k-11]+b
}

>x1.new<-ts(data=X1,frequency=12)
>HW1<-HoltWinters(x1.new)
>predHW1<-predict(HW1,12)
>C1<-X1
>for (k in 1 : 12 ) {
C1[n1+k]=predHW1[k]
}

length(C1)

P12<-P1[97:108]
X12<-X[97:108]
C12<-C1[97:108]

>ts.plot(P12,col=”red”)
>lines(C12,col=”yellow”)
>lines(X12,col=”blue”)

Plottando soltanto i valori del 2008, quelli veri e quelli previsti con i due modelli, otteniamo il seguente risultato mostrato in figura dove con il colore blu sono rappresentati i dati della vera serie storica con il rosso i dati previsti con il terzo modello lineare e in giallo i dati previsti con il modello Holt-Winters.

Confronto previsioni

Di prima impressione i dati previsti con il modello lineare sembrano essere più aderenti a quelli della serie storica vera.

In tutti i casi conviene calcolare la deviazione standard sigma usandola come metodo per verificare l’errore e per vedere quale dei due modelli di previsione si è comportato meglio in questo caso presentando una deviazione minore.

>sigma1<-sqrt(mean((P12-X12)^2))
>sigmaHW<-sqrt(mean((C12-X12)^2))
>sigma1
>sigmaHW

Sigma calcolato con R vale circa 0.49 per il modello lineare e 0.58 con il modello HW quindi l’intuizione grafica era corretta.

Intervallo di confidenza

Il valore di sigma può essere utilizzato per stabilire un intervallo di confidenza entro il quale accettare future previsioni.

Ipotizzando errori gaussiani e un largo margine di confidenza possiamo ipotizzare che il valore vero futuro ricada in una previsione pk + o – un valore 3*sigma scelto tra quelli precedentemente calcolato.

Scegliendo il sigma calcolato con il modello lineare si può decidere che i prossimi valori della serie storica molto probabilmente rientreranno in un valore p compreso tra pk+1.47 e pk-1.47.

Abbiamo visto soltanto qualche altra piccola funzionalità offerta da R.

Cercando dati a proposito della serie storica sulla co2 ho trovato un articolo molto interessante che riguarda lo strano crescere del trend delle rilevazioni noaa http://globalwarming.blog.meteogiornale.it/2009/06/17/anomalie-oceaniche-lo-strano-caso-dei-dati-noaa/ . (sempre in tema Copenaghen)

Vi ringrazio dell’attenzione e vi lascio esplorare R sperando, che al di là della propria maggiore o minore esperienza con il calcolo statistico, riusciate a sfruttare questo ambiente open source, a presto

Andrea Stani