setwd("C:/Users/USER/OneDrive/UFMT")

require(forecast) #pacote para analise de series temporais
require(snpar) #pacote para fazer teste de tendencia
source("Series/funcoes auxiliares.R") #funcoes auxiliares
require(tidyverse)
require(highcharter) #pacote para gráfico

Descrição:

Série temporal referente à produção de soja em toneladas no estado de Mato Grosso, no anos de 1988 a 2020.

Fonte: IBGE

Conjunto de dados

# Dados producao de soja em toneladas no estado de Mato Grosso

dados <- xlsx::read.xlsx("C:/Users/USER/OneDrive/UFMT/SERIES/Dados_Soja.xlsx",1)
DT:: datatable(dados)

Dados convertidos em série temporal:

#Converter em serie temporal
y1 <- ts(dados$producao,start=1988,frequency=1) #Frequecia anual
y1
## Time Series:
## Start = 1988 
## End = 2020 
## Frequency = 1 
##  [1]  2694718  3795435  3064715  2738410  3642743  4118726  5319793  5491426
##  [9]  5032921  6060882  7228052  7473028  8774470  9533286 11684885 12965983
## [17] 14517912 17761444 15594221 15275087 17802976 17962819 18787783 20800544
## [25] 21841292 23416774 26495884 27850954 26277303 30479870 31608562 32242463
## [33] 35070044

Gráfico da Série

hchart(dados,type = "line", hcaes(x = Ano, y = producao),
       color = "green", name = "Toneladas") %>% 
       hc_title(text = "Produção de Soja em Mato Grosso") %>% 
       hc_yAxis(title = list(text = "Produção (em toneladas)"))

Gráfico da Funcao de autocorrelacao e autocorrelacao parcial:

##Funcao de autocorrelacao e autocorrelacao parcial
par(mfrow=c(1,2))
acf(y1, main = "FAC")
pacf(y1, main = "FACP")

par(mfrow=c(1,1))

Pelo gráfico da série temporal e os gráficos de da função de autocorrelação parcial, a série aparenta ter tendência.

Testes estatísticos

Teste homogeneidade da variância:

gma(y1,k=8)  

##     correlacao valor-p
## cor     -0.127  0.8387

Com p-valor > 0.05, não rejeita-se a hipótese de variância constante. Ou seja, não é necessária a transformação da série.

Teste de tendência:

Teste de Cox-Stuart

cs.test(y1) 
## 
##  Exact Cox-Stuart trend test
## 
## data:  y1
## S = 0, p-value = 3.052e-05
## alternative hypothesis: data have a monotonic trend

p-valor < 0,01.

Dado o p-valor, a série possui tendência.

Retirar Tendência:

Método de regressão

n <- length(y1)
t <- seq(1,n,1)

model <- lm(y1 ~ t) 
summary(model)
## 
## Call:
## lm(formula = y1 ~ t)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2720837 -1635136  -667641  1737043  4064982 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2299044     732235   -3.14   0.0037 ** 
## t            1014748      37579   27.00   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2056000 on 31 degrees of freedom
## Multiple R-squared:  0.9592, Adjusted R-squared:  0.9579 
## F-statistic: 729.2 on 1 and 31 DF,  p-value: < 2.2e-16
e <- residuals(model) %>% 
      as.ts() 

cs.test(e)
## 
##  Exact Cox-Stuart trend test
## 
## data:  e
## S = 8, p-value = 0.8036
## alternative hypothesis: data have a monotonic trend

p-valor > 0.05

Foi possível remover a tendência pelo método de regressão.

Método da Diferença

yd1 <- diff(y1,1)
cs.test(yd1)
## 
##  Exact Cox-Stuart trend test
## 
## data:  yd1
## S = 6, p-value = 0.4545
## alternative hypothesis: data have a monotonic trend

p-valor > 0.05

Foi possível remover a tendência com 1 diferença.

par(mfrow=c(1,2))
plot(e, main = "Método Regressão")
plot(yd1, main = "Método Diferença")

par(mfrow=c(1,1))

Teste de Sazonalidade

Teste Fisher

#Periodograma 
P <- periodograma(yd1)
plot(P,type="l", main = "Periodograma")

#Teste
Fisher.test(P)
##        g valor-p periodo
##  0.30978 0.06145     3.2

Dado p-valor > 0.05. A série não possui sazonalidade.

Conclusão: A Série possui tendência, mas não possui sazonalidade.

Funcao de autocorrelacao e autocorrelacao parcial

##Funcao de autocorrelacao e autocorrelacao parcial
par(mfrow=c(1,2))
acf(yd1, main = "FAC")
pacf(yd1, main = "FACP")

par(mfrow=c(1,1))

Modelo

Modelo de Suavizacao Exponencial de Holt

##Ajuste do modelo
seh=HoltWinters(y1,gamma=FALSE)
seh
## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = y1, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 0.7872417
##  beta : 0
##  gamma: FALSE
## 
## Coefficients:
##       [,1]
## a 34718246
## b  1100717
#Valores ajustados
#fitted(seh) 

##Plotar o valores ajustados
plot(seh, ylab = "Produção(t)", main = "Observado vs Ajustado", xlab = "Ano")

Predição:

Predição para o período de 5 anos:

#Predicao
p=predict(seh, 5, prediction.interval = TRUE)
p
## Time Series:
## Start = 2021 
## End = 2025 
## Frequency = 1 
##           fit      upr      lwr
## 2021 35818963 38448451 33189475
## 2022 36919680 40266213 33573147
## 2023 38020397 41955414 34085380
## 2024 39121114 43567397 34674831
## 2025 40221831 45126371 35317292
plot(seh, p, xlab = "Ano", ylab = "Produção(t)", main = "Modelo SEH - Previsão")

Modelo ARIMA

Modelo Arima(1,1,1)

m <- Arima(y1, order = c(1,1,1)) 
m
## Series: y1 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.9930  -0.9099
## s.e.  0.0161   0.0906
## 
## sigma^2 estimated as 2.004e+12:  log likelihood=-498.27
## AIC=1002.54   AICc=1003.4   BIC=1006.94
summary.arima(m)
##     estimativa erro padrão estatistica t valorp
## ar1  0.9930070  0.01611723      61.61153      0
## ma1 -0.9098836  0.09059367     -10.04357      0

Predição:

Predição para o período de 5 anos:

prev <- forecast(m, h=5, level = 95)
prev
##      Point Forecast    Lo 95    Hi 95
## 2021       36282890 33507812 39057968
## 2022       37487255 33395825 41578685
## 2023       38683197 33466676 43899719
## 2024       39870777 33608856 46132697
## 2025       41050051 33781620 48318482
plot(prev, main = "ARIMA(1,1,1) - Previsão")

Diagnóstico:

##Diagnostico do modelo
tsdiag(m)