O estudo desse laboratório consiste na análise da série temporal apresentada no artigo apresentado em sala. Iremos estudar o conjunto de dados e utilizar as técnicas conhecidas para poder propor um modelo SARIMA, analisar os resíduos do modelo, e por fim, realizar previsões. O conjunto de dados a seguir consiste em medidas (em dBm) de microondas do início de 2006 até o fim de 2010.

##        MONTH YEAR2006 YEAR2007 YEAR2008 YEAR2009 YEAR2010
## 1    January -71.0625 -70.5000 -74.3750 -74.0625 -63.9375
## 2   February -67.8750 -71.2500 -69.3750 -73.8125 -62.8750
## 3      March -26.5625 -28.2500 -28.8125 -29.0625 -28.7500
## 4      April -27.3125 -27.1875 -29.6875 -27.3750 -27.0000
## 5        May -27.7500 -27.9375 -27.3750 -27.6875 -27.7500
## 6       June -28.0625 -27.4375 -30.2500 -28.0625 -27.8125
## 7       July -27.6875 -27.6875 -26.8125 -27.6250 -27.1875
## 8     August -28.2500 -30.7500 -27.1875 -31.1875 -27.1875
## 9  September -26.3750 -26.3750 -26.3750 -26.3750 -26.3125
## 10   October -26.6250 -26.6250 -26.3750 -26.6250 -32.8125
## 11  November -26.3750 -26.3750 -26.7500 -26.3750 -30.2500
## 12  December -71.3750 -73.5625 -73.6250 -63.1875 -71.8750


Agora, iremos plotar o gráfico desta série temporal para observar o seu comportamento e realizaremos o teste Augmented Dickey-Fuller para verificar se a série é estacionária.

datats = as.vector(unlist(microwavedata[,2:6]))
microwavets = ts(datats,start=2006,frequency=12)
plot(microwavets)

adf.test(microwavets)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  microwavets
## Dickey-Fuller = -3.599, Lag order = 3, p-value = 0.04089
## alternative hypothesis: stationary


Verificamos que a série não apresenta tendência(confirmamos a estacionaridade da série através do p-valor do teste ser < 0.05), porém há claramente um componente sazonal. O ciclo observado é a cada 12 meses. A ACF e PACF será calculada a seguir e irá nos ajudar nesta análise.

acf2(microwavets)

##         ACF  PACF
##  [1,]  0.51  0.51
##  [2,]  0.07 -0.25
##  [3,] -0.30 -0.32
##  [4,] -0.29  0.07
##  [5,] -0.29 -0.24
##  [6,] -0.26 -0.23
##  [7,] -0.30 -0.24
##  [8,] -0.30 -0.36
##  [9,] -0.32 -0.59
## [10,]  0.11  0.10
## [11,]  0.50  0.08
## [12,]  0.80  0.37
## [13,]  0.41 -0.20
## [14,]  0.04 -0.05
## [15,] -0.24 -0.07
## [16,] -0.22  0.01
## [17,] -0.24 -0.05
## [18,] -0.22 -0.06
## [19,] -0.25  0.19
## [20,] -0.25  0.08
## [21,] -0.26  0.01
## [22,]  0.10 -0.09
## [23,]  0.39 -0.16
## [24,]  0.60 -0.05
## [25,]  0.28 -0.09
## [26,]  0.03  0.00
## [27,] -0.16  0.01
## [28,] -0.16  0.02
## [29,] -0.16  0.04
## [30,] -0.15 -0.08
## [31,] -0.17 -0.01
## [32,] -0.17  0.01
## [33,] -0.18  0.06
## [34,]  0.07 -0.03
## [35,]  0.27  0.01
## [36,]  0.38 -0.09
## [37,]  0.16 -0.04
## [38,]  0.00 -0.02
## [39,] -0.08  0.01
## [40,] -0.08  0.00
## [41,] -0.09  0.02
## [42,] -0.08  0.00
## [43,] -0.09 -0.02
## [44,] -0.10 -0.02
## [45,] -0.10  0.06
## [46,]  0.06  0.05
## [47,]  0.15 -0.01
## [48,]  0.18 -0.03

Com o efeito sazonal verificado, o artigo sugere aplicar a diferença para1 sazonalidade. Como o ciclo observado era de 12 meses, aplicaremos a primeira diferença com lag = 12. Depois, através do ACF e PACF desta nova série, iremos observar as ordens do modelo SARIMA.

plot(diff(microwavets,12))

acf2(diff(microwavets,12),max.lag = 40)

##         ACF  PACF
##  [1,]  0.48  0.48
##  [2,]  0.22 -0.02
##  [3,]  0.00 -0.12
##  [4,] -0.18 -0.17
##  [5,] -0.03  0.20
##  [6,]  0.05  0.05
##  [7,]  0.05 -0.06
##  [8,]  0.06  0.00
##  [9,] -0.15 -0.20
## [10,] -0.32 -0.22
## [11,] -0.27  0.01
## [12,] -0.35 -0.21
## [13,]  0.03  0.31
## [14,]  0.10 -0.10
## [15,]  0.02 -0.11
## [16,]  0.16  0.23
## [17,]  0.04  0.03
## [18,] -0.07 -0.19
## [19,]  0.01  0.02
## [20,] -0.07 -0.10
## [21,] -0.01 -0.10
## [22,]  0.03 -0.08
## [23,] -0.03  0.09
## [24,]  0.00 -0.11
## [25,] -0.11  0.00
## [26,] -0.12  0.03
## [27,] -0.04 -0.04
## [28,] -0.10  0.02
## [29,] -0.03 -0.05
## [30,]  0.06 -0.09
## [31,] -0.01  0.05
## [32,]  0.01 -0.14
## [33,]  0.00 -0.02
## [34,] -0.07 -0.03
## [35,] -0.01 -0.01
## [36,] -0.04 -0.11
## [37,]  0.01  0.07
## [38,]  0.03 -0.01
## [39,]  0.02 -0.01
## [40,]  0.03 -0.06

Observando o comportamento da ACF e PACF da diff(microwavets,12), uma sugestão para o modelo SARIMA é: p=1, d=1, q=1, P=0, D=1, Q=1. O artigo também sugere a análise de outros modelos, assim iremos construir alguns modelos e depois compará-los.

fit = sarima(microwavets,1,1,1,0,1,1,12, details = F)
fita1 = sarima(microwavets,1,1,1,0,1,2,12, details = F)
fita2 = sarima(microwavets,1,1,1,1,1,1,12, details = F)
fita3 = sarima(microwavets,1,1,1,1,1,0,12, details = F)
fita4 = sarima(microwavets,1,1,1,1,1,1,12, details = F)
  


fit$AICc
## [1] 4.009714
fita1$AICc
## [1] 4.038327
fita2$AICc
## [1] 4.034571
fita3$AICc
## [1] 4.01665
fita4$AICc
## [1] 4.034571


Após escolher o modelo com o menor AICc, iremos analisar os resíduos do modelo para verificar se o modelo é apropriado.

fit = sarima(microwavets,1,1,1,0,1,1,12, details = T,no.constant=TRUE)
## initial  value 1.223586 
## iter   2 value 1.057562
## iter   3 value 0.927674
## iter   4 value 0.927196
## iter   5 value 0.923603
## iter   6 value 0.923004
## iter   7 value 0.922641
## iter   8 value 0.922610
## iter   9 value 0.922609
## iter  10 value 0.922608
## iter  10 value 0.922608
## iter  10 value 0.922608
## final  value 0.922608 
## converged
## initial  value 1.034816 
## iter   2 value 1.029600
## iter   3 value 1.026995
## iter   4 value 1.020333
## iter   5 value 1.010982
## iter   6 value 1.002459
## iter   7 value 0.973885
## iter   8 value 0.965681
## iter   9 value 0.965594
## iter  10 value 0.965494
## iter  11 value 0.965470
## iter  12 value 0.965382
## iter  13 value 0.965339
## iter  14 value 0.965332
## iter  15 value 0.965328
## iter  16 value 0.965326
## iter  17 value 0.965322
## iter  18 value 0.965316
## iter  19 value 0.965310
## iter  20 value 0.965306
## iter  21 value 0.965306
## iter  22 value 0.965306
## iter  22 value 0.965306
## iter  22 value 0.965306
## final  value 0.965306 
## converged

fit
## $fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), include.mean = !no.constant, transform.pars = trans, fixed = fixed, 
##     optim.control = list(trace = trc, REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1     ma1     sma1
##       0.5705  -1.000  -0.9995
## s.e.  0.1239   0.115   1.0435
## 
## sigma^2 estimated as 4.235:  log likelihood = -112.06,  aic = 232.12
## 
## $degrees_of_freedom
## [1] 44
## 
## $ttable
##      Estimate     SE t.value p.value
## ar1    0.5705 0.1239  4.6058  0.0000
## ma1   -1.0000 0.1150 -8.6927  0.0000
## sma1  -0.9995 1.0435 -0.9579  0.3434
## 
## $AIC
## [1] 4.002051
## 
## $AICc
## [1] 4.009714
## 
## $BIC
## [1] 4.129648


Após verificar o modelo, iremos fazer previsões para 2 anos depois de 2010:

sarima.for(microwavets,24,1,1,1,0,1,1,12)

## $pred
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2011 -71.54565 -69.45433 -28.50965 -27.82363 -27.74787 -28.33692 -27.39165
## 2012 -70.87214 -69.06488 -28.28224 -27.68867 -27.66565 -28.28479 -27.35667
##            Aug       Sep       Oct       Nov       Dec
## 2011 -28.89299 -26.33735 -27.78541 -27.19902 -70.70355
## 2012 -28.86781 -26.31776 -27.76900 -27.18443 -70.69000
## 
## $se
##           Jan      Feb      Mar      Apr      May      Jun      Jul
## 2011 2.344666 2.738568 2.883439 2.944971 2.973908 2.988821 2.997414
## 2012 3.129092 3.158743 3.174178 3.182652 3.187626 3.190942 3.193739
##           Aug      Sep      Oct      Nov      Dec
## 2011 3.003363 3.008813 3.015386 3.024522 3.036998
## 2012 3.196920 3.201451 3.208574 3.219905 3.237026