1 Ajustando um modelo SARIMA

Queremos ajustar um modelo SARIMA na série z5.

#source("ajud.R", local = knitr::knit_global())
source("ajud.R")
load("zt.Rdata")
N <- length(z5)
zz <- z5
z <- z5[1:(N-10)]
w <- diff(diff(z, 12))

par(mfrow = c(1,2))
plot(ts(z))
plot(ts(w))

Sendo o modelo sarima:

fit <- arima(w, order=c(4,0,0), seasonal=list(order=c(2,0,0), period=12),
             include.mean=F)

2 Qualidade do ajuste

Verificando os gráficos de autocorrelações e autocorrelações parciais do resíduos, temos que:

plotap(fit)

Verificamos também a significância dos coeficientes estimados

signi(fit)
##             Coef     DesPad          T             P
## ar1   0.21868766 0.04499813   4.859928  1.174285e-06
## ar2  -0.09999916 0.04376555  -2.284883  2.231969e-02
## ar3  -0.32068348 0.04390224  -7.304490  2.783203e-13
## ar4  -0.11554557 0.04505042  -2.564806  1.032335e-02
## sar1  0.18278542 0.03314053   5.515464  3.478606e-08
## sar2 -0.68792605 0.03135846 -21.937496 1.140072e-106
all(signi(fit)[,4] < .05)
## [1] TRUE

Como podemos ver, não precisa zerar nenhum dos parâmetros.

Verificamos, então, o ajuste do modelo usando o teste Ljung-Box

knitr::kable(Ljung_Box(fit,6))
Lag \(\chi^2\) GdL P-valor
12 8.788248 6 0.1858417
24 28.856911 18 0.0501561
36 40.436825 30 0.0966772
48 48.105514 42 0.2393623

Como podemos ver, o modelo é bem ajustado (todos quatro testes não rejeitam a hipóse nula).

3 Modelo estimado:

Verificamos o modelo estimado

fit
## 
## Call:
## arima(x = w, order = c(4, 0, 0), seasonal = list(order = c(2, 0, 0), period = 12), 
##     include.mean = F)
## 
## Coefficients:
##          ar1      ar2      ar3      ar4    sar1     sar2
##       0.2187  -0.1000  -0.3207  -0.1155  0.1828  -0.6879
## s.e.  0.0450   0.0438   0.0439   0.0451  0.0331   0.0314
## 
## sigma^2 estimated as 9.165:  log likelihood = -1246.12,  aic = 2506.24

Logo,

\[\begin{align} &(1 - 0.2187B + 0.1000B^2 +0.3207B^3 +0.1155B^4)(1- 0.1828B^{12} + 0.6879B^{24})(1-B^{12})(1-B)Z_t = a_t, \\ \end{align}\] sendo \(\sigma^2_a = 9.165\).

Obtendo os valores no R:

p0 <- polynomial(c(1,rep(0,11),-1)) #diferenciação sazional 
pd <-  polynomial(c(1,-1))
coe<-fit$coef[1:4] # AUTOREGRESSIVO 
p1<-polynomial(c(1,-coe))
coe <- fit$coef[c(5,6)]
p2<-polynomial(c(1,rep(0,11),-coe[1], rep(0,11), -coe[2]))

coef <- coef(p0*p1*p2*pd)[-1]
theta0 <- 0

4 Previsões

O interesse agora é realizar previsões de um a dez passos a frente com o modelo estimado.

4.1 Previsão sem atualização

hmax <- 10  
prevm_ar_sa(z,zz,hmax,theta0,coef, fit$sigma2)
## $u
##           prev     orig      erro       dp      inf      sup
##  [1,] 2127.146 2125.270  1.876101 3.027313 2121.212 2133.079
##  [2,] 2150.316 2148.328  1.988669 4.772412 2140.963 2159.670
##  [3,] 2171.400 2166.893  4.506985 5.936892 2159.764 2183.036
##  [4,] 2185.883 2182.468  3.415886 6.426315 2173.288 2198.479
##  [5,] 2199.246 2193.025  6.220462 6.642083 2186.227 2212.264
##  [6,] 2201.468 2193.295  8.172994 6.829696 2188.082 2214.855
##  [7,] 2223.619 2216.900  6.718793 7.119384 2209.665 2237.573
##  [8,] 2216.069 2217.290 -1.220391 7.540686 2201.290 2230.849
##  [9,] 2216.261 2219.965 -3.703826 7.998071 2200.585 2231.938
## [10,] 2205.706 2208.801 -3.095612 8.389732 2189.262 2222.150
## 
## $EQMP
## [1] 21.48804

O erro quadrático médio da previsão é \(21.48804\)

4.2 Previsão com atualização

prevm_ar_ca(z, zz, 10, theta0, coef, fit$sigma2) 
## $u
##           prev     orig        erro       dp      inf      sup
##  [1,] 2127.146 2125.270  1.87610055 3.027313 2121.212 2133.079
##  [2,] 2148.030 2148.328 -0.29771195 3.027313 2142.097 2153.964
##  [3,] 2169.574 2166.893  2.68130763 3.027313 2163.641 2175.508
##  [4,] 2181.439 2182.468 -1.02893223 3.027313 2175.505 2187.372
##  [5,] 2196.573 2193.025  3.54788899 3.027313 2190.640 2202.507
##  [6,] 2195.346 2193.295  2.05068547 3.027313 2189.413 2201.280
##  [7,] 2215.241 2216.900 -1.65965667 3.027313 2209.307 2221.174
##  [8,] 2210.637 2217.290 -6.65260762 3.027313 2204.704 2216.571
##  [9,] 2220.023 2219.965  0.05754928 3.027313 2214.089 2225.956
## [10,] 2208.918 2208.801  0.11666638 3.027313 2202.985 2214.852
## 
## $EQMP
## [1] 7.56779

O erro quadrático médio da previsão é \(7.56779\)