Forecast Arima
Fagner Sutel
03/07/2019
Carga de Dados
rm(list = ls(all.names = TRUE))
setwd("/Users/fagne/Desktop/")
fatais = read.csv("PLAN REG VÍTIMAS FATAIS.csv", header = T, sep = ";")
#fatais = fatais[10:19,]
fatais## ANO FATAIS
## 1 2000 163
## 2 2001 140
## 3 2002 156
## 4 2003 170
## 5 2004 173
## 6 2005 161
## 7 2006 157
## 8 2007 155
## 9 2008 148
## 10 2009 170
## 11 2010 143
## 12 2011 146
## 13 2012 105
## 14 2013 127
## 15 2014 141
## 16 2015 100
## 17 2016 92
## 18 2017 90
## 19 2018 75dim(fatais)## [1] 19 2x = fatais$ANO
y = fatais$FATAIS
ajuste = lm(y ~ x)
summary(ajuste)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.053 -13.474 5.842 9.303 32.526
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9072.2368 1495.1047 6.068 1.25e-05 ***
## x -4.4474 0.7442 -5.976 1.50e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.77 on 17 degrees of freedom
## Multiple R-squared: 0.6775, Adjusted R-squared: 0.6585
## F-statistic: 35.71 on 1 and 17 DF, p-value: 1.503e-05options(scipen = 999)
summary(ajuste)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.053 -13.474 5.842 9.303 32.526
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9072.2368 1495.1047 6.068 0.0000125 ***
## x -4.4474 0.7442 -5.976 0.0000150 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.77 on 17 degrees of freedom
## Multiple R-squared: 0.6775, Adjusted R-squared: 0.6585
## F-statistic: 35.71 on 1 and 17 DF, p-value: 0.00001503predict(lm(y ~ x))## 1 2 3 4 5 6 7
## 177.50000 173.05263 168.60526 164.15789 159.71053 155.26316 150.81579
## 8 9 10 11 12 13 14
## 146.36842 141.92105 137.47368 133.02632 128.57895 124.13158 119.68421
## 15 16 17 18 19
## 115.23684 110.78947 106.34211 101.89474 97.44737novo = c(2019, 2020,2021, 2022, 2023, 2024, 2025, 2026, 2027, 2028, 2029, 2030)
novo <- data.frame(x = novo)
novo## x
## 1 2019
## 2 2020
## 3 2021
## 4 2022
## 5 2023
## 6 2024
## 7 2025
## 8 2026
## 9 2027
## 10 2028
## 11 2029
## 12 2030res = predict(lm(y ~ x), novo, se.fit = TRUE)
res$fit## 1 2 3 4 5 6 7 8
## 93.00000 88.55263 84.10526 79.65789 75.21053 70.76316 66.31579 61.86842
## 9 10 11 12
## 57.42105 52.97368 48.52632 44.07895res$fit = as.integer(as.character(res$fit))
res$fit## [1] 92 88 84 79 75 70 66 61 57 52 48 44#INTERVALO DE PREVISAO
int_prev <- predict(lm(y ~ x), novo, interval = "prediction")
int_prev## fit lwr upr
## 1 93.00000 51.458365 134.54163
## 2 88.55263 46.392477 130.71279
## 3 84.10526 41.277930 126.93260
## 4 79.65789 36.116962 123.19883
## 5 75.21053 30.911815 119.50924
## 6 70.76316 25.664716 115.86160
## 7 66.31579 20.377857 112.25372
## 8 61.86842 15.053376 108.68347
## 9 57.42105 9.693348 105.14876
## 10 52.97368 4.299772 101.64760
## 11 48.52632 -1.125434 98.17807
## 12 44.07895 -6.580439 94.73833#INTERVALO DE CONFIANÇA
int_conf <- predict(lm(y ~ x), novo, interval = "confidence")
int_conf## fit lwr upr
## 1 93.00000 75.09781 110.90219
## 2 88.55263 69.25860 107.84667
## 3 84.10526 63.39388 104.81665
## 4 79.65789 57.50854 101.80725
## 5 75.21053 51.60635 98.81470
## 6 70.76316 45.69025 95.83606
## 7 66.31579 39.76255 92.86903
## 8 61.86842 33.82507 89.91177
## 9 57.42105 27.87930 86.96280
## 10 52.97368 21.92645 84.02092
## 11 48.52632 15.96749 81.08515
## 12 44.07895 10.00323 78.15467library(ggplot2)## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlangp <- ggplot(fatais, aes(ANO, FATAIS)) +
geom_point()
# adiciona a regressao
p + geom_smooth(method = lm)
# ajuste de regressao local
p + geom_smooth(method = "loess")
predito= cbind(novo, res$fit)
predito## x res$fit
## 1 2019 92
## 2 2020 88
## 3 2021 84
## 4 2022 79
## 5 2023 75
## 6 2024 70
## 7 2025 66
## 8 2026 61
## 9 2027 57
## 10 2028 52
## 11 2029 48
## 12 2030 44names(predito) = c("ANO", "FATAIS")
predito## ANO FATAIS
## 1 2019 92
## 2 2020 88
## 3 2021 84
## 4 2022 79
## 5 2023 75
## 6 2024 70
## 7 2025 66
## 8 2026 61
## 9 2027 57
## 10 2028 52
## 11 2029 48
## 12 2030 44totais = rbind(fatais, predito)
totais## ANO FATAIS
## 1 2000 163
## 2 2001 140
## 3 2002 156
## 4 2003 170
## 5 2004 173
## 6 2005 161
## 7 2006 157
## 8 2007 155
## 9 2008 148
## 10 2009 170
## 11 2010 143
## 12 2011 146
## 13 2012 105
## 14 2013 127
## 15 2014 141
## 16 2015 100
## 17 2016 92
## 18 2017 90
## 19 2018 75
## 20 2019 92
## 21 2020 88
## 22 2021 84
## 23 2022 79
## 24 2023 75
## 25 2024 70
## 26 2025 66
## 27 2026 61
## 28 2027 57
## 29 2028 52
## 30 2029 48
## 31 2030 44p <- ggplot(totais, aes(ANO, FATAIS)) +
geom_point()
# adiciona a regressao
p + geom_smooth(method = lm)
# ajuste de regressao local
p + geom_smooth(method = "loess")
dados.obs = read.csv("PLAN REG VÍTIMAS FATAIS.csv", header = T, sep = ";")
dados.obs## ANO FATAIS
## 1 2000 163
## 2 2001 140
## 3 2002 156
## 4 2003 170
## 5 2004 173
## 6 2005 161
## 7 2006 157
## 8 2007 155
## 9 2008 148
## 10 2009 170
## 11 2010 143
## 12 2011 146
## 13 2012 105
## 14 2013 127
## 15 2014 141
## 16 2015 100
## 17 2016 92
## 18 2017 90
## 19 2018 75Plotando dados
Fatais de 2000 a 2018
dados.ts = ts(dados.obs$FATAIS, frequency = 19, start = c(1, 19))
plot(dados.ts)
Predict
Predict para 10 anos:
library(forecast)## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Registered S3 methods overwritten by 'forecast':
## method from
## fitted.fracdiff fracdiff
## residuals.fracdiff fracdiffpredict(dados.obs$FATAIS)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 66.476546 56.415502 76.537591 51.089506 81.863586
## 21 55.041584 46.207627 63.875542 41.531213 68.551956
## 22 43.606622 35.303624 51.909621 30.908282 56.304963
## 23 32.171661 23.489169 40.854152 18.892936 45.450385
## 24 20.736699 10.830456 30.642941 5.586408 35.886989
## 25 9.301737 -2.424385 21.027859 -8.631820 27.235293
## 26 -2.133225 -16.048385 11.781934 -23.414625 19.148174
## 27 -13.568187 -29.900636 2.764262 -38.546512 11.410138
## 28 -25.003149 -43.908888 -6.097410 -53.916982 3.910684
## 29 -36.438111 -58.044809 -14.831413 -69.482704 -3.393517Arrima
model = auto.arima(dados.ts)
model## Series: dados.ts
## ARIMA(0,1,0)
##
## sigma^2 estimated as 375.6: log likelihood=-78.9
## AIC=159.79 AICc=160.04 BIC=160.68Projeção h=30
Impressão de gráfico com série temporal e projeção com intervalos de confiança de 80 e 95% Incluir no gráfico os valores do modelo para dados observados em azul
forec1 = forecast(model, level = c(80,95), h=30)
{plot(forec1)
lines(fitted(model),col="blue")}
Verificando os dados
class(forec1)## [1] "forecast"forec = as.data.frame(forec1)
forec1## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2.947368 75 50.1644356 99.83556 37.0172819 112.9827
## 3.000000 75 39.8772080 110.12279 21.2843250 128.7157
## 3.052632 75 31.9835406 118.01646 9.2120025 140.7880
## 3.105263 75 25.3288712 124.67113 -0.9654361 150.9654
## 3.157895 75 19.4659897 130.53401 -9.9319396 159.9319
## 3.210526 75 14.1655397 135.83446 -18.0382783 168.0383
## 3.263158 75 9.2912729 140.70873 -25.4928261 175.4928
## 3.315789 75 4.7544159 145.24558 -32.4313501 182.4314
## 3.368421 75 0.4933067 149.50669 -38.9481542 188.9482
## 3.421053 75 -3.5369505 153.53695 -45.1119008 195.1119
## 3.473684 75 -7.3702486 157.37025 -50.9744243 200.9744
## 3.526316 75 -11.0329188 161.03292 -56.5759950 206.5760
## 3.578947 75 -14.5459010 164.54590 -61.9486376 211.9486
## 3.631579 75 -17.9261731 167.92617 -67.1183176 217.1183
## 3.684211 75 -21.1877274 171.18773 -72.1064345 222.1064
## 3.736842 75 -24.3422577 174.34226 -76.9308723 226.9309
## 3.789474 75 -27.3996554 177.39966 -81.6067585 231.6068
## 3.842105 75 -30.3683761 180.36838 -86.1470251 236.1470
## 3.894737 75 -33.2557155 183.25572 -90.5628297 240.5628
## 3.947368 75 -36.0680206 186.06802 -94.8638791 244.8639
## 4.000000 75 -38.8108539 188.81085 -99.0586806 249.0587
## 4.052632 75 -41.4891228 191.48912 -103.1547394 253.1547
## 4.105263 75 -44.1071828 194.10718 -107.1587166 257.1587
## 4.157895 75 -46.6689206 196.66892 -111.0765566 261.0766
## 4.210526 75 -49.1778221 199.17782 -114.9135903 264.9136
## 4.263158 75 -51.6370276 201.63703 -118.6746206 268.6746
## 4.315789 75 -54.0493782 204.04938 -122.3639925 272.3640
## 4.368421 75 -56.4174543 206.41745 -125.9856522 275.9857
## 4.421053 75 -58.7436075 208.74361 -129.5431966 279.5432
## 4.473684 75 -61.0299886 211.02999 -133.0399148 283.0399dim(forec1)## NULLSem h
forec = forecast(model, level = c(80,95))
{plot(forec)
lines(fitted(model),col="blue")}
Verificando os dados residuais
#class(forec)
#forec = as.data.frame(forec)
#forec
#dim(forec)fit = auto.arima(dados.obs$FATAIS, seasonal = F)
tsdisplay(residuals(fit), lag.max = 45, main = 'Modelo de resíduos (1,1,1)')
fit2 = arima(dados.obs$FATAIS, order=c(1,1,7))
tsdisplay(residuals(fit2), lag.max = 45, main = 'Modelo de resíduos (1,1,7)')
Autocorrelação e correlação parcial
acf(dados.obs$FATAIS, lag.max = 18, type = c("correlation", "covariance",
"partial"), plot = TRUE, na.action = na.contiguous, demean = TRUE)
pacf(dados.obs$FATAIS, lag.max = 18, plot = TRUE, na.action = na.contiguous,
demean = TRUE)## Warning in plot.window(...): "demean" não é um parâmetro gráfico## Warning in plot.xy(xy, type, ...): "demean" não é um parâmetro gráfico## Warning in axis(side = side, at = at, labels = labels, ...): "demean" não é
## um parâmetro gráfico
## Warning in axis(side = side, at = at, labels = labels, ...): "demean" não é
## um parâmetro gráfico## Warning in box(...): "demean" não é um parâmetro gráfico## Warning in title(...): "demean" não é um parâmetro gráfico
Sem Arima
fcast = forecast(dados.obs$FATAIS, h=10)
plot(fcast)
Modelando com comparativo de buffer de dados reais em um intervalo de 6 anos (1/3 dos dados)
Reservando dados reais para comparar com o modelo:
dados = window(ts(dados.obs$FATAIS), start= 14)
ajuste_controle = arima(ts(dados.obs$FATAIS[-c(14:19)]), order = c(1,1,7))
forecast_ajuste = forecast(ajuste_controle, h=17, level = c(80))
{plot(forecast_ajuste)
lines(ts(dados.obs$FATAIS))}
forecast_ajuste## Point Forecast Lo 80 Hi 80
## 14 123.61156 107.391778 139.8313
## 15 96.58497 77.091045 116.0789
## 16 115.84523 88.621250 143.0692
## 17 86.60136 55.816028 117.3867
## 18 99.78787 62.211409 137.3643
## 19 71.09866 29.642780 112.5545
## 20 86.98577 34.448994 139.5225
## 21 73.80535 15.341074 132.2696
## 22 84.74022 18.497871 150.9826
## 23 75.66832 4.289750 147.0469
## 24 83.19464 5.624585 160.7647
## 25 76.95058 -5.260341 159.1615
## 26 82.13084 -5.325855 169.5875
## 27 77.83314 -13.897665 169.5639
## 28 81.39864 -14.948465 177.7457
## 29 78.44059 -21.886954 178.7681
## 30 80.89468 -23.599266 185.3886Reservando dados reais para comparar com o modelo em intervalos de confiança de 80% e 95%:
dados = window(ts(dados.obs$FATAIS), start= 14)
ajuste_controle = arima(ts(dados.obs$FATAIS[-c(14:19)]), order = c(0,1,7))
forecast_ajuste = forecast(ajuste_controle, h=17, level = c(80, 95))
{plot(forecast_ajuste)
lines(ts(dados.obs$FATAIS))}
forecast_ajuste## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 14 125.68108 108.882518 142.4796 99.989897 151.3723
## 15 103.32428 84.011831 122.6367 73.788439 132.8601
## 16 118.16484 92.817056 143.5126 79.398749 156.9309
## 17 95.66851 67.495297 123.8417 52.581297 138.7557
## 18 102.63235 68.506320 136.7584 50.441093 154.8236
## 19 84.20918 46.384726 122.0336 26.361669 142.0567
## 20 92.58936 45.331478 139.8472 20.314669 164.8640
## 21 92.58936 39.970496 145.2082 12.115755 173.0630
## 22 92.58936 35.107344 150.0714 4.678205 180.5005
## 23 92.58936 30.624696 154.5540 -2.177413 187.3561
## 24 92.58936 26.445147 158.7336 -8.569483 193.7482
## 25 92.58936 22.514440 162.6643 -14.580979 199.7597
## 26 92.58936 18.792804 166.3859 -20.272731 205.4514
## 27 92.58936 15.250048 169.9287 -25.690908 210.8696
## 28 92.58936 11.862620 173.3161 -30.871532 216.0502
## 29 92.58936 8.611721 176.5670 -35.843354 221.0221
## 30 92.58936 5.482063 179.6966 -40.629753 225.8085Auto Arima
Aplica-se o ajuste:
fit = auto.arima(dados.obs$FATAIS)Plota-se os dados
plot(forecast(fit, h=30))
Resume-se os dados
summary(fit)## Series: dados.obs$FATAIS
## ARIMA(0,1,0)
##
## sigma^2 estimated as 375.6: log likelihood=-78.9
## AIC=159.79 AICc=160.04 BIC=160.68
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -4.623 18.86242 14.53489 -5.248565 11.83527 0.9479279
## ACF1
## Training set -0.3425583Previsão
forecast(fit, h=13)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 75 50.1644356 99.83556 37.0172819 112.9827
## 21 75 39.8772080 110.12279 21.2843250 128.7157
## 22 75 31.9835406 118.01646 9.2120025 140.7880
## 23 75 25.3288712 124.67113 -0.9654361 150.9654
## 24 75 19.4659897 130.53401 -9.9319396 159.9319
## 25 75 14.1655397 135.83446 -18.0382783 168.0383
## 26 75 9.2912729 140.70873 -25.4928261 175.4928
## 27 75 4.7544159 145.24558 -32.4313501 182.4314
## 28 75 0.4933067 149.50669 -38.9481542 188.9482
## 29 75 -3.5369505 153.53695 -45.1119008 195.1119
## 30 75 -7.3702486 157.37025 -50.9744243 200.9744
## 31 75 -11.0329188 161.03292 -56.5759950 206.5760
## 32 75 -14.5459010 164.54590 -61.9486376 211.9486Correlação
print(fit_RSq <- cor(fit$fitted, fit$x)^2)## [1] 0.626536