%%%%% Esercizio. I file prezzi1.txt e prezzi2.txt descrivono gli indici dei prezzi al consumo (tipo NIC, mensili) di alimentari ed energia, relativamente agli anni 1996-2005. Sono presi dal sito Istat. Conviene familiarizzare con il problema cercando i dati sul sito e leggento un po' di definizioni collegate. 1) Caricare i dati dei due file, visualizzarli, applicare a ciascuno alcuni metodi di previsione descrivendo i risultati trovati. 2) Analizzare il legame tra i due indici, anche traslati di tempo, vedendo se sia possibile prevederne uno sulla base dell'altro a tempi precedenti o allo stesso tempo. 3) Mettere insieme le due analisi cercando di prevedere l'indice degli alimentari sualla base dei valori passati di alimentari e energia. 4) Facoltativamente, cercare sul sito Istat altri dati che aiutino nella previsione. 5) Esaminare l'efficacia dei modelli di previsione trovati, applicandoli agli anni successivi. %%%%% Frammenti di soluzione A<-read.table(file="prezzi1.txt") E<-read.table(file="prezzi2.txt") A <- A[,1] E <- E[,1] ts.plot(A) ts.plot(E) ts.plot(c(99,A,140)) lines(E, col="red") %%% Considerazioni visive: ..... %%% Metodi SE ecc... %%% Autocorrelazione: acf(A,40) acf(E,40) %%% in entrambi i casi non ci sono picchi, cioè periodicità, ma non è un white noise, c'è una memoria, più lunga per gli alimentari, quindi usiamo un ARMA legato ai dati di precedenti. Ecco il risultato (per un ARMA facile) con parametri lievemente ottimizzati per tentativi: L <- length(A) pA <- matrix(nrow=L) pA[1] <- 0 pA[2] <- 0 a <- 1.7 b <- 1-a c <- 0 for (k in 3:L) { pA[k] <- a*A[k-1]+b*A[k-2]+c } ts.plot(A) lines(pA, col="red") sqrt((1/(L-2))*sum((pA[3:L]-A[3:L])^2)) Calcoliamo anche l'analogo di R^2: > 1-(1/(L-2))*sum((pA[3:L]-A[3:L])^2)/sd(A)^2 [1] 0.9955642 E' molto alto, è un fit ottimo. Vediamo per E: pE <- matrix(nrow=L) pE[1] <- 0 pE[2] <- 0 a <- 1.2 b <- 1-a c <- 0 for (k in 3:L) { pE[k] <- a*E[k-1]+b*E[k-2]+c } ts.plot(E) lines(pE, col="red") sqrt((1/(L-2))*sum((pE[3:L]-E[3:L])^2)) %%% Si può usare la regressione per automatizzare alcune ricerche: x1 <- A[2:(L-1)] x2 <- A[1:(L-2)] y <- A[3:L] reg1 <- lm(y ~ x1+x2) summary(reg1) Il risultato è ancora migliore di quello ottenuto a mano (Multiple R-Squared: 0.9988). Proviamo a migliorarlo: x1 <- A[3:(L-1)] x2 <- A[2:(L-2)] x3 <- A[1:(L-3)] y <- A[4:L] reg2 <- lm(y ~ x1+x2+x3) > summary(reg2) Call: lm(formula = y ~ x1 + x2 + x3) Residuals: Min 1Q Median 3Q Max -0.638421 -0.159322 -0.005136 0.115791 1.006831 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.03672 0.37068 0.099 0.921 x1 1.60566 0.09402 17.078 < 2e-16 *** x2 -0.68034 0.16645 -4.087 8.19e-05 *** x3 0.07507 0.09396 0.799 0.426 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.2582 on 113 degrees of freedom Multiple R-Squared: 0.9988, Adjusted R-squared: 0.9988 F-statistic: 3.094e+04 on 3 and 113 DF, p-value: < 2.2e-16 %%% si vede che x3 è inessenziale, mentre x2 era fondamentale. %%% Legame tra le due serie: > cor(A,E) [1] 0.8552864 le due serie sono molto correlate, nonostante ad occhio questo non sia evidente. Vediamo se la correlazione aumenta ritardandone una. Dobbiamo però accorciare le serie. Chiamiamo N l'accorciamento. N<-50 cor <- matrix(nrow=(N+1)) for (G in 0:N) { cor[G+1]<-cor(A[(1+G):(120-N+G)],E[1:(120-N)]) } plot(cor) max(cor) > max(cor) [1] 0.9305253 > > which.max(cor) [1] 14 quindi con un ritardo di circa 14, che in pratica può essere 12 (un anno), c'è una correlazione molto alta. Costruiamo un modello di regressione che utilizzi sia i due mesi passati (che erano molto importanti) sia l'indice dell'energia: x1 <- A[12:(L-1)] x2 <- A[11:(L-2)] e12 <- A[1:(L-12)] y <- A[13:L] reg3 <- lm(y ~ x1+x2+e12) summary(reg3) %%% anche qui si vede che e1 è inessenziale. Invece, osserviamo che si ottengono ottimi risultati sostituendo x1 e x2 con e12: e12 <- A[1:(L-12)] y <- A[13:L] reg4 <- lm(y ~ e12) summary(reg4) %%% Multiple R-Squared: 0.9391. %%% C'è però un'osservazione da fare. Anche il modello banale Media Mobile con finestra uno, MM1, x1 <- A[2:(L-1)] y <- A[3:L] reg5 <- lm(y ~ x1) summary(reg5) funziona bene: Multiple R-Squared: 0.9982. Quindi la bontà del modello va confrontata con questa. In quest'ottica, il modello y <- e12 è pessimo. Invece il modello x1 <- A[12:(L-1)] e12 <- A[1:(L-12)] y <- A[13:L] reg6 <- lm(y ~ x1+e12) summary(reg6) ha un Multiple R-Squared: 0.9984, migliore del modello banale MM1, e si vede che e12 ha due stelle. %%%%% Cerchiamo di capire più intuitivamente questi modelli. Da un lato, il modello banale MM1 ha già un buon R^2. Si tratta della previsione: p[k]=A[k-1]. Davvero la previsione p.reg1[k] <- (1.56)*A[k-1]+(-0.56)*A[k-2]-0.08 è migliore? R^2 direbbe di sì. Ma intuitivamente, perché dev'essere così? La ragione è che (1.56)*A[k-1]+(-0.56)*A[k-2] è un'estrapolazione dei valori A[k-1] e A[k-2] nella direzione di A[k-1]. Ed il grafico mostra che il comportamento è molto spesso di tipo estrapolatorio. Inoltre, a livello numerico, possiamo osservare a confronto i due metodi con i comandi: Ver <- matrix(nrow=51,ncol=3) k0<-10 for (k in k0:(k0+50)) { Ver[k-k0+1,1]<-A[k] Ver[k-k0+1,2]<-A[k-1] Ver[k-k0+1,3]<-(1.56)*A[k-1]+(-0.56)*A[k-2] } Ver %%% Torniamo invece al modello p.reg6[k] <- (1.05)*A[k-1]+(-0.05)*E[k-12] Perché intuitivamente dovrebbe essere migliore di quello banale? Perché c'è una generale tendenza a crescere ( (1.05)*A[k-1] ); i dati mostrano che dove E[k-12] è grande, la crescita è inferiore; però questo sembra un po' un caso. Si potrebbe invece immaginare, dai grafici, che se E cresce 12 mesi prima, allora A dovrebbe crescere di più. Proviamo il modello (abbiamo cambiato il ritardo ottimizzandoloi per tentativi): x1 <- A[13:(L-1)] e13 <- A[2:(L-12)] e14 <- A[1:(L-13)] y <- A[14:L] reg7 <- lm(y ~ x1+e13+e14) summary(reg7) > summary(reg7) Call: lm(formula = y ~ x1 + e13 + e14) Residuals: Min 1Q Median 3Q Max -0.690427 -0.168879 0.005689 0.165127 0.940441 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.52818 0.44148 1.196 0.234294 x1 1.03864 0.01640 63.339 < 2e-16 *** e13 0.25598 0.09304 2.751 0.007015 ** e14 -0.29896 0.08799 -3.398 0.000967 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.2756 on 103 degrees of freedom Multiple R-Squared: 0.9986, Adjusted R-squared: 0.9985 F-statistic: 2.416e+04 on 3 and 103 DF, p-value: < 2.2e-16 Il risultato è migliore di quello con solo e12 ed il modello mostra, come prevedevamo, che y tende a crescere rispetto a x1, ed ancor di più se 12 mesi prima circa, l'incremento di E era positivo.