%%%%%% TENTATIVI DI PREVISIONE IF <- read.table(file="import_food.txt") ts.plot(IF) % sono dati mensili sulle importazioni di beni alimentari tra due nazioni. IF.new<-ts(data = IF, frequency = 12) IF.W <- HoltWinters(IF.new) plot(IF.W) % così abbiamo visto i dati con le previsioni ottimizzate. % Ecco i valori numerici delle previsioni, insieme con s,m,F: > IF.W$fitted xhat level trend season Jan 2 7507.887 7619.629 34.09057 -145.832639 Feb 2 7055.686 7675.505 34.10908 -653.928472 Mar 2 7771.454 7725.523 34.12259 11.809028 Apr 2 8312.149 7770.724 34.13199 507.292361 ..... ..... ..... Feb 11 10620.920 11208.256 33.98016 -621.315632 % e se vogliamo solo i valori previsti: > IF.W$fitted[,1] Jan Feb Mar Apr May Jun Jul 2 7507.887 7055.686 7771.454 8312.149 7671.842 8426.428 7508.912 3 8100.718 7596.198 8275.261 8815.989 8174.776 8949.605 8040.228 4 8560.218 8093.104 8829.362 9283.897 8660.387 9374.977 8469.070 5 8779.441 8186.993 8882.425 9352.483 8898.212 9652.658 8569.023 6 9112.320 8720.427 9337.242 9779.040 9633.251 10383.737 9025.246 7 9587.906 9412.191 9859.228 10354.782 10229.319 10825.607 9716.736 8 9892.185 9692.595 9921.828 10609.818 10480.126 10942.997 10034.719 9 10259.901 9975.591 10303.762 11108.438 10943.242 11305.840 10412.786 10 10739.440 10364.751 10950.839 11672.715 11262.629 11713.234 10743.693 11 11016.537 10620.920 Aug Sep Oct Nov Dec 2 7524.584 8036.424 8494.482 8238.416 8301.632 3 8070.397 8587.318 9086.270 8799.327 8826.127 4 8363.566 8938.026 9228.578 8889.109 9042.896 5 8602.750 9507.006 9410.271 9166.837 9551.184 6 9306.979 10137.777 9704.198 9934.220 10084.313 7 9601.402 10486.541 10020.926 10462.398 10477.875 8 9913.544 10740.098 10402.541 10683.742 10659.408 9 10315.030 11094.121 10874.751 10939.454 11030.077 10 10632.131 11565.150 11186.435 11137.593 11461.146 11 %%% Calcoliamo a mano l'errore quadratico medio, ricordando che sono 122 dati a cui vanno tolti i primi 12, come si vede da: length(IF.W$fitted[,1]) % che fornisce il valore 110. IF <- IF[,1] L<- length(IF) MSE <- (1/(L-12))*sum((IF.W$fitted[1:(L-12),1]-IF[13:L])^2) IF.W$SSE/(L-12) % forniscono lo stesso valore. Quindi è proprio vero che R minimizza MSE (chiamato SSE). % Nota 1: calcola MSE su tutta la finestra possibile; se volessimo usarne una più breve, dovremmo fare a mano. % Nota 2: non minimizza MAPE: se volessimo usarlo, dovremmo fare a mano. sqrt(IF.W$SSE/(L-12)) % = 307.9, è l'analogo della devizione standard. Su dati dell'ordine dei 10.000 è un errore percentuale abbastanza piccolo. %%%%% Metodo tipo ARMA elementare: a <- 1.5 b <- 1-a c <- 0 p <- matrix(nrow=L) p[1:12]<- IF[1:12] for (k in 13:L) { p[k] <- a*IF[k-1]+b*IF[k-12]+c } ts.plot(IF) lines(p, col="red") sqrt((1/(L-12))*sum((p[13:L]-IF[13:L])^2)) % è peggiore di Winters. Si cerchino valori migliori dei parametri, per tentativi, cambiando a,b,c, e copiano ogni volta le seguenti istruzioni: a <- 0.2 b <- 1-a c <- 200 p <- matrix(nrow=L) p[1:12]<- IF[1:12] for (k in 13:L) { p[k] <- a*IF[k-1]+b*IF[k-12]+c } ts.plot(IF) lines(p, col="red") sqrt((1/(L-12))*sum((p[13:L]-IF[13:L])^2)) % qui abbiamo già inserito dei valori di a,b,c, che migliorano molto l'errore: = 359.6. % Non è ancora come Winters. Ne trovate uno migliore?