Considerações Preliminares
Modelos ARIMA (AR - Auto Regressivos;
I - Integrados; MA - Médias
Móveis)
Evidentemente, para entender os modelos, precisamos entender o
que são: (a) Processos Autoregressivos; (b) Processos Integrados e; (c)
Processos de Médias Móveis.
Eles representam uma família de modelos para previsão de Séries
temporais.
Juntos com os modelos de Suavização Exponencial foram os
modelos mais difundidos para previsões. Essa família de modelos consegue
fazer previsão através do que chamamos de autocorrelação dos dados (O
passado traz uma memória para o presente).
Dessa forma, antes de avançar, vamos começar discutindo
Estacionariedade e técnicas de
Diferenciação de Séries Temporais.
Considerações Preliminares: Estacionariedade
Uma série temporal estacionária é aquela na qual suas
propriedades não dependem do momento que observamos as séries.
Assim, séries com tendência ou sazonalidade não devem ser
estacionárias, uma vez que ambos os fatores afetam o valor que
observamos a série ao longo do tempo.
Um exemplos de série estacionária é o que chamamos de ruído
branco. Ela é uma série estocástica, sem memória temporal e não importa
em que momento no tempo iremos observá-la, as propriedades são as
mesmas.
Ruído Branco: \[y_t=\epsilon_t,
\;\;\;\epsilon_t \sim N(0,1)\]
Considerações Preliminares: Estacionariedade
Em geral, uma série estacionária não terá padrões reconhecíveis
no longo prazo. O gráfico da série contra o tempo irá mostrar uma série
praticamente horizontal, com variância constante.
Formalmente, uma série temporal é Fracamente
Estacionária se são satisfeitas as seguintes
definições:
- \(E|y_t|^2<\infty\) (Variância
Finita)
- \(E(y_t)=\mu,\;\;\forall t\in
\mathbb{Z}\) (Média igual para qualquer período)
- \(E(y_t-\mu).(y_{t-j}-\mu)=\gamma_j\)
(Variância é igual e a Autocovariância depende apenas da distância
temporal entre as observações)
- O que seria então uma série não estacionária? Visualmente é uma
série com tendência (seja ela estocástica ou determinística).
- Veja os exemplo de Tendência Estocástica:
Passeio Aleatório ou Random Walk
\[y_t=y_{t-1}+\epsilon_t\;\;\;\epsilon_t \sim
N(0,1)\]
- Veja o exemplo de Tendência Determinística: \[y_t=(0.02).t+\epsilon_t\;\;\;\epsilon_t \sim
N(0,2)\]
Diferenciação
Vamos agora verificar visualmente a questão da estacionariedade
em uma série econômica.
Escolhemos o IBC-BR sem dessazonalização, por ser um indicador de
atividade econômica (proxy do PIB).
Uma maneira de tentar transformar uma série não estacionária em
estacionária é aplicando a 1ª Diferença: \[\Delta y_t=y_t-y_{t-1}\]
Séries em Diferenças geralmente podem ser utilizadas, desde que a
variância não continue crescendo constantemente ao longo do tempo.
Quando isso ocorre, podemos trabalhar com taxas de crescimento: \[g_y=\frac{\Delta
y_t}{y_{t-1}}=y_t/y_{t-1}-1\] ou na sua aproximação por uso de
logaritmo: \[g_y \approx
\Delta(\ln{y_t})=\ln{y_t}-\ln{y_{t-1}}\]
library(ipeadatar)
ibc_br <- ipeadata(code = "SGS12_IBCBR12")
ibc_ts <- ts(ibc_br$value, start = c(2003,1), frequency = 12)
plot(ibc_ts, main = "Índice em Nível do IBC-BR",
ylab = "Número Índice",
xlab = "Tempo",
col = "blue")
plot(diff(ibc_ts), main = "Índice em Primeira Diferença do IBC-BR",
ylab = "1ª Diff do Número Índice",
xlab = "Tempo",
col = "tomato")
plot(diff(ibc_ts)/stats::lag(ibc_ts), main = "Taxa de Crescimento do IBC-BR",
ylab = "Taxa de Crescimento do Número Índice",
xlab = "Tempo",
col = "magenta")
Diferenciação
Observações:
Séries não estacionárias devem apresentar um
lento decaimento na Função de Autocorrelação
(FAC)
Séries estacionárias podem apresentar ou não
um truncamento na FAC, mas não o deicaimento
lento.
acf(ibc_ts,
main = "FAC - IBC-BR em Número Índice",
ylab = "FAC",
xlab = "Defasagem (anos)")
acf(diff(ibc_ts),
main = "FAC - 1ª Diferença do IBC-BR em Número Índice",
ylab = "FAC",
xlab = "Defasagem (anos)")
acf(diff(ibc_ts)/stats::lag(ibc_ts),
main = "FAC - Taxa de Crescimento do IBC-BR em Número Índice",
ylab = "FAC",
xlab = "Defasagem (anos)")
Testes Estatísticos para Estacionariedade:
Uma forma mais objetiva de determinar se a diferenciação é
necessária para deixar a série estacionária é usar um teste de
raiz unitária.
Esses testam a hipótese estatística de que as séries possuam raiz
unitária. Se sim, aceitamos a hipótese da série ser não estacionária. A
sua rejeição é a aceitação da estacionariedade.
Um grande números de testes pode ser utilizado para esse fim.
Alguns exemplos:
- Teste ADF (Dickey-Fuller Aumentado)
- Teste PP (Phillips-Perron)
- Teste KPSS (Kwiatkowski-Phillips-Schmidt-Shin)
- Testando a raiz unitária no IBC-BR (em nível) através do ADF:
library(urca)
### Especificação 1 - Sem constante e sem tendência
summary(ur.df(ibc_ts, type = "none"))
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.4631 -2.9664 -0.6138 2.2527 18.2527
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 0.001338 0.002309 0.579 0.562893
## z.diff.lag -0.244903 0.062384 -3.926 0.000113 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.786 on 243 degrees of freedom
## Multiple R-squared: 0.06011, Adjusted R-squared: 0.05237
## F-statistic: 7.77 on 2 and 243 DF, p-value: 0.0005359
##
##
## Value of test-statistic is: 0.5793
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62
### Especificação 2 - Com constante e sem tendência
summary(ur.df(ibc_ts, type = "drift"))
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.3226 -3.0473 -0.8233 1.9921 18.6332
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.41650 2.96940 2.834 0.004979 **
## z.lag.1 -0.06184 0.02241 -2.760 0.006222 **
## z.diff.lag -0.21979 0.06214 -3.537 0.000485 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.718 on 242 degrees of freedom
## Multiple R-squared: 0.08859, Adjusted R-squared: 0.08106
## F-statistic: 11.76 on 2 and 242 DF, p-value: 1.335e-05
##
##
## Value of test-statistic is: -2.76 4.1896
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81
### Especificação 3 - Com constante e com tendência
summary(ur.df(ibc_ts, type = "trend"))
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.3689 -2.9342 -0.7385 2.0385 17.2586
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.265649 3.829980 3.725 0.000244 ***
## z.lag.1 -0.120330 0.033081 -3.637 0.000337 ***
## tt 0.015011 0.006296 2.384 0.017893 *
## z.diff.lag -0.189160 0.062869 -3.009 0.002901 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.673 on 241 degrees of freedom
## Multiple R-squared: 0.1096, Adjusted R-squared: 0.09851
## F-statistic: 9.887 on 3 and 241 DF, p-value: 3.577e-06
##
##
## Value of test-statistic is: -3.6375 4.7419 6.7247
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2 6.22 4.75 4.07
## phi3 8.43 6.49 5.47
CONCLUSÃO: Aceitamos a H0. Portanto, em nível, o
IBC-BR é não estacionário.
Testes Estatísticos para Estacionariedade:
- Agora testando na 1ª Diferença do IBC-BR.
library(urca)
### Especificação 1 - Sem constante e sem tendência
summary(ur.df(diff(ibc_ts), type = "none"))
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.0341 -2.6483 -0.3731 2.4174 16.7861
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -1.61202 0.09672 -16.667 < 2e-16 ***
## z.diff.lag 0.29528 0.06134 4.814 2.61e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.573 on 242 degrees of freedom
## Multiple R-squared: 0.6555, Adjusted R-squared: 0.6527
## F-statistic: 230.2 on 2 and 242 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -16.6674
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62
### Especificação 2 - Com constante e sem tendência
summary(ur.df(diff(ibc_ts), type = "drift"))
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.3459 -2.9849 -0.6784 2.1065 16.4600
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.31911 0.29329 1.088 0.278
## z.lag.1 -1.61902 0.09689 -16.709 < 2e-16 ***
## z.diff.lag 0.29885 0.06141 4.867 2.05e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.571 on 241 degrees of freedom
## Multiple R-squared: 0.6572, Adjusted R-squared: 0.6543
## F-statistic: 231 on 2 and 241 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -16.7092 139.5988
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81
### Especificação 3 - Com constante e com tendência
summary(ur.df(diff(ibc_ts), type = "trend"))
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.1584 -3.0684 -0.7691 2.2258 16.7212
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.597220 0.593167 1.007 0.315
## z.lag.1 -1.621467 0.097143 -16.692 < 2e-16 ***
## tt -0.002248 0.004165 -0.540 0.590
## z.diff.lag 0.300101 0.061542 4.876 1.97e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.578 on 240 degrees of freedom
## Multiple R-squared: 0.6576, Adjusted R-squared: 0.6533
## F-statistic: 153.6 on 3 and 240 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -16.6916 92.8892 139.3336
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2 6.22 4.75 4.07
## phi3 8.43 6.49 5.47
CONCLUSÃO: Rejeitamos a H0. Portanto, em 1ª
Diferença, o IBC-BR é estacionário.
Dessa forma, se a série em nível é não estacionária e a sua
primeira diferença é estacionária, temos uma série que possui
ordem de integração igual a 1.
Se fosse necessário trabalhar com a segunda diferença da série
para estacionarizar, teríamos uma série com ordem de
integração igual a 2.
ps.: Em economia, dificilmente temos séries que são I(2)!
Assim, resolvemos a discussão para a determinação do termo
\(I\) do
modelo ARIMA.
Modelo (ou parte) Autoregressivo(a)
Em uma regressão linear múltipla, nós prevemos a variável de
interesse com base em uma combinação linear das variáveis explicativas
(preditores).
Em um modelo autoregressivo, nós prevemos a
variável de interesse usando uma combinação linear dos valores
defasados da variável de interesse.
O termo autoregressivo indica uma regressão de uma variável
contra ela mesma.
Assim, podemos escrever um modelo autoregressivo de ordem \(p\) como: \[y_t=c+\phi_1.y_{t-1}+\phi_2.y_{t-2}+\ldots+\phi_p.y_{t-p}+\epsilon_t,
\] onde \(\epsilon_t\) é um
ruído branco.
Iremos nos referir como modelo \(AR(p)\), um modelo autoregressivo
de ordem \(p\).
Alguns comportamentos podem surgir conforme o coeficiente dos
modelos. No caso de um \(AR(1)\), temos
que:
- Quando \(\phi_1=0\), temos um ruído
branco.
- QUando \(\phi_1=1\) e \(c=0\) temos um passeio aleatório.
- Quando \(\phi_1=1\) e \(c\neq0\), temos um passeio aleatório com
drift (tendência).
- Quando \(\phi_1<0\), \(y_t\) tende a oscilar em torno da
média.
- Os modelos Autoregressivos podem ter restrições em um ou no conjunto
de parâmetros.
Exemplo de Simulação de um Modelo \(AR(2)\), com \(\phi_1=0.65,\phi_2=0.3\)
# Simulação de um modelo AR(2), phi_1 = 0.65, phi_2 = 0.3.
ts.sim <- arima.sim(list(order = c(2,0,0),
ar = c(0.65,0.3)),
n = 2000)
ts.plot(ts.sim, col = "purple",
main = "Modelo AR(2) Simulado",
ylab = "Valor",
xlab = "Tempo")
Modelos (ou parte) de Médias Móveis
- Agora, ao invés de utilizar valores defasados da variável de
interesse, o modelo de médias móveis usa os erros de previsão passados
como variáveis explicativas.
\[y_t=c+\epsilon_t+\theta_1.\epsilon_{t-1}+\theta_2.\epsilon_{t-2}+\ldots+\theta_q.\epsilon_{t-q},\]
onde \(\epsilon_t\) é um ruído
branco. Chamaremos então de modelo \(MA(q)\), como sendo um modelo de
médias móveis de ordem \(q\).
Claro que nós não observamos os valores de \(\epsilon_t\), portanto não podemos entender
como uma regressão no sentido usual.
O modelo é chamado de médias móveis porque pode ser visto como
uma média ponderada dos erros passados. Não confundir com modelos de
suavização por médias móveis. (Esse último é utilizado para estimar
ciclos e tendências).
Exemplo de Simulação de um Modelo \(MA(3)\), com \(\theta_1=0.70,\theta_2=0.2,\theta_3=-0.10\)
# Simulação de um modelo MA(3), theta_1 = 0.7, theta_2 = 0.2, theta_3 = -0.1.
ts.sim <- arima.sim(list(order = c(0,0,3),
ma = c(0.7,0.2,-0.1)),
n = 2000)
ts.plot(ts.sim, col = "darkcyan",
main = "Modelo MA(3) Simulado",
ylab = "Valor",
xlab = "Tempo")
Juntando as Partes - Modelo ARIMA
Precisamos agora juntar o que aprendemos com os Modelos AR, Ordem
de Integração e Modelos MA.
O ARIMA é a combinação dos três, de modo que dizemos que temos um
modelo \(ARIMA(p,d,q)\), como sendo um
modelo que possui \(q\) ordem de
defasagem no efeito AR; \(d\) ordem de
integração; e \(q\) ordem de defasagem
para a parte MA.
O modelo completo pode ser escrito como: \[y'_t=c+\phi_1.y'_{t-1}+\ldots+\phi_p.y'_{t-p}+\theta_1.\epsilon_{t-1}+\ldots+\theta_q.\epsilon_{t-q}+\epsilon_t,
\] onde \(y'_t\) é a série
temporal em diferença (pode ser que seja preciso mais de uma
diferença).
Casos Especiais do ARIMA:
| \(ARIMA(0,0,0)\) |
Ruído Branco |
| \(ARIMA(0,1,0)\) |
Passeio Aleatório c/s constante |
| \(ARIMA(p,0,0)\) |
Autoregressão |
| \(ARIMA(0,0,q)\) |
Média Móvel |
Sazonalidade nos Dados
O que é sazonalidade? É um evento cíclico com regularidade definida e
esperança matemática nula no término do ciclo.
Na prática, temos em dados mensais, meses em que há esperança de
elevação dos resultados acima da média e em outros meses há a esperaça
de queda dos resultados abaixo da média.
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
ggsubseriesplot(ibc_ts,
main='Sub-série Sazonal para o IBC-BR')
ggseasonplot(ibc_ts,
polar=TRUE,
main='Sub-série Sazonal por ano para o IBC-BR em escala Polar')
Modelos SARIMA
Até agora havíamos nos preocupado apenas com modelos sem o fator
Sazonalidade.
A sazonalidade, no âmbito da econometria tradicional, deve ser
tratada de modo a obter uma série dessazonalida.
Ocorre que nosso objetivo aqui é lidar com previsões. Com esse
intuito, não iremos tratar a série e sim melhorar a nossa previsão com
base na carga de sazonalidade prevista para um determinado mês.
Para isso, iremos agora lidar com modelos
SARIMA, onde o \(S\) indica justamente a
sazonalidade. Agora o modelo para a ser descrito como:
\[ARIMA(p,d,q)(P,D,Q)_m\]
onde \(m\) é o número de observações
por ano. Usamos agora a notação maiúscula para a parte sazonal do modelo
que pode também ser modela por efeito AR e MA.
- Para o diagnóstico dos efeitos, devemos fazer uso das:
- Função de Autocorrelação (FAC)
- Função de Autocorrelação Parcial (FACP)
Modelos SARIMA
Exemplos:
- Modelo \(ARIMA(0,0,0)(0,0,1)_{12}\)
- SMA(1)
o comportamento observado de ser:
- Um truncamento (pico) na 12ª Defasagem da FAC, sem nenhum outro
truncamento relevante.
- Decaimento exponencial nas defasagens sazonais da FACP (isto é, 12,
24, 36, etc).
## Carregando pacotes exigidos: stats4
##
## Attaching package: 'sarima'
## The following object is masked from 'package:stats':
##
## spectrum
# SAR(1), 12 seasons
SAR <- sim_sarima(n=2000,
model=list(sar=0.8,
nseasons=12,
sigma2 = 1))
# FAC e FACP do SAR
acf(SAR, lag.max = 60, main = "FAC do SAR(1)")
pacf(SAR, lag.max = 60, main = "FACP do SAR(1)")
Modelos SARIMA
Exemplos:
- Modelo \(ARIMA(0,0,0)(1,0,0)_{12}\)
- SAR(1)
o comportamento observado de ser:
- Decaimento exponencial nas defasagens sazonais da FAC (isto é, 12,
24, 36, etc).
- Um truncamento (pico) na 12ª Defasagem da FACP, sem nenhum outro
truncamento relevante.
# SMA(1), 12 Meses
SMA <- sim_sarima(n=2000,
model=list(sma=0.8,
nseasons=12,
sigma2 = 1))
# FAC e FACP do SMA
acf(SMA, lag.max = 60, main = "FAC do SMA(1)")
pacf(SMA, lag.max = 60, main = "FACP do SMA(1)")
Sazonalidade - SARIMA - IBC-BR
- Vamos agora verificar a Sazonalidade do IBC-BR via FAC e FACP (em
nível e em 1ª diferença).
# FAC e FACP do IBC-BR
acf(ibc_ts, lag.max = 60, main = "FAC do IBC-BR")
pacf(ibc_ts, lag.max = 60, main = "FACP do IBC-BR")
# FAC e FACP da 1ª Diff do IBC-BR
acf(diff(ibc_ts), lag.max = 60, main = "FAC do Diff IBC-BR")
pacf(diff(ibc_ts), lag.max = 60, main = "FACP do Diff IBC-BR")
A Escolha das Ordens \(p\) e \(q\).
Até agora definimos como escolher os termos \(d\) através dos testes de raízes unitárias,
bem como do comportamento da FAC.
Vimos os padrões que mostram a sazonalidade para os modelos SAR
e/ou SMA.
Agora como proceder com o truncamento de \(p\) e \(q\)?
Os dados podem ser de um \(ARIMA(p,d,0)\) se a FAC e a FACP dos dados
em diferença obtiverem o seguinte comportamento:
- A FAC tiver decaimento exponencial ou sinoidal;
- Houver um truncamento na defasagem \(p\) na FACP, mas nenhum além da \(p\).
Os dados podem ser de um \(ARIMA(0,d,q)\) se a FAC e a FACP dos dados
em diferença obtiverem o seguinte comportamento:
- A FACP tiver decaimento exponencial ou sinoidal;
- Houver um truncamento na defasagem \(q\) na FAC, mas nenhum além da \(q\).
Cuidado: Por vezes os gráficos das FAC e FACP podem
não ser tão conclusivos. Nesse caso, tenha parcimônia!
- No nosso modelo, temos a suspeita até agora que tenhamos um \(ARIMA(2,1,2)(1,1,1)_{12}\).
Estimação do Modelo
Estimação por Máxima Verossimilhança
Uma vez que o modelo tenha sido identificado, ou seja, as ordens
de \(p,d,q\), nós precisamos estimar os
parâmetros \(c,\phi_1,\ldots,\phi_p,\theta_1,\ldots,\theta_q\).
Quando estimamos o modelo ARIMA no R, ele usa a
Estimação de Máxima Verossimilhança (EMV).
Essa estimador encontra os valores dos parâmetros que maximizem a
probabilidade de obtermos os dados que temos obvervado.
Para os modelos ARIMA, a EMV é similar ao que teríamos por MQO,
ao minimizar: \[\sum_{t=1}^T
\epsilon_t^2\]
Na prática, o R reporta o valor da log
verossimilhança dos dados, ou seja, o logaritmo da
probabilidade dos dados observados serem oriundos do modelo
estimado.
Para dados valores de \(p,d,q\)
a função de estimação irá tentar maximizar o logaritmo da
verossimilhança quendo os parâmetros forem estimados.
Fit.sarima1 <- arima(ibc_ts, order = c(2,1,2), seasonal = c(1,1,0))
Fit.sarima2 <- arima(ibc_ts, order = c(2,1,2), seasonal = c(0,1,1))
Fit.sarima3 <- arima(ibc_ts, order = c(2,1,2), seasonal = c(1,1,1))
#Fit.sarima4 <- arima(ibc_ts, order = c(3,1,3), seasonal = c(1,1,0))
Fit.sarima5 <- arima(ibc_ts, order = c(3,1,3), seasonal = c(0,1,1))
#Fit.sarima6 <- arima(ibc_ts, order = c(3,1,3), seasonal = c(1,1,1))
summary(Fit.sarima1)
##
## Call:
## arima(x = ibc_ts, order = c(2, 1, 2), seasonal = c(1, 1, 0))
##
## Coefficients:
## ar1 ar2 ma1 ma2 sar1
## -0.9928 -0.6974 0.7621 0.4229 -0.4064
## s.e. 0.1198 0.1513 0.1548 0.1990 0.0610
##
## sigma^2 estimated as 9.364: log likelihood = -594.97, aic = 1201.93
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.002219466 2.978512 2.159694 -0.01508774 1.622025 0.5980308
## ACF1
## Training set -0.0169677
stargazer::stargazer(list(Fit.sarima1,
Fit.sarima2,
Fit.sarima3,
#Fit.sarima4,
Fit.sarima5),
type = 'html')
|
|
|
|
Dependent variable:
|
|
|
|
|
|
ibc_ts
|
|
|
(1)
|
(2)
|
(3)
|
(4)
|
|
|
|
ar1
|
-0.993***
|
-1.013***
|
-0.970***
|
-0.633***
|
|
|
(0.120)
|
(0.119)
|
(0.127)
|
(0.241)
|
|
|
|
|
|
|
|
ar2
|
-0.697***
|
-0.812***
|
-0.722***
|
-0.377
|
|
|
(0.151)
|
(0.157)
|
(0.198)
|
(0.280)
|
|
|
|
|
|
|
|
ar3
|
|
|
|
0.525**
|
|
|
|
|
|
(0.241)
|
|
|
|
|
|
|
|
ma1
|
0.762***
|
0.826***
|
0.757***
|
0.481**
|
|
|
(0.155)
|
(0.186)
|
(0.172)
|
(0.205)
|
|
|
|
|
|
|
|
ma2
|
0.423**
|
0.604***
|
0.488*
|
0.128
|
|
|
(0.199)
|
(0.226)
|
(0.263)
|
(0.241)
|
|
|
|
|
|
|
|
sar1
|
-0.406***
|
|
0.099
|
|
|
|
(0.061)
|
|
(0.089)
|
|
|
|
|
|
|
|
|
ma3
|
|
|
|
-0.724***
|
|
|
|
|
|
(0.202)
|
|
|
|
|
|
|
|
sma1
|
|
-0.850***
|
-0.884***
|
-0.853***
|
|
|
|
(0.052)
|
(0.061)
|
(0.053)
|
|
|
|
|
|
|
|
|
|
Observations
|
234
|
234
|
234
|
234
|
|
Log Likelihood
|
-594.967
|
-569.097
|
-568.483
|
-560.107
|
|
sigma2
|
9.364
|
7.105
|
7.041
|
6.425
|
|
Akaike Inf. Crit.
|
1,201.934
|
1,150.193
|
1,150.967
|
1,136.213
|
|
|
|
Note:
|
p<0.1; p<0.05;
p<0.01
|
Alternativa de Estimação
E se eu não quiser ter todo esse trabalho para escolher o melhor
modelo?
Tem maneira rápida? Tem a função auto.arima() do
pacote forecast.
Ela não é um remédio universal! Não há garantia de que seja o
melhor modelo.
Todavia, o algoritmo de Hyndman-Khandakar (2008) é implementado
para selecionar o modelo. Aqui é escolhido o melhor modelo até a segunda
ordem \(p,d,q\), com base no ajuste
pela minimização do AICc.
fit_auto <- auto.arima(ibc_ts)
summary(fit_auto)
## Series: ibc_ts
## ARIMA(2,1,2)(1,1,1)[12]
##
## Coefficients:
## ar1 ar2 ma1 ma2 sar1 sma1
## -0.9697 -0.7218 0.7570 0.4885 0.0988 -0.8841
## s.e. 0.1272 0.1982 0.1719 0.2625 0.0888 0.0614
##
## sigma^2 = 7.226: log likelihood = -568.48
## AIC=1150.97 AICc=1151.46 BIC=1175.15
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.06664837 2.582726 1.850837 -0.06819569 1.396626 0.3805575
## ACF1
## Training set -0.02440057
Previsão
Tendo estimado o modelo, podemos agora realizar as
previsões.
Quanto mais no futuro, mais imprecisa é a nossa previsão, tendo
em vista que ela irá realizar projeçoes sobre projeções.
A previsão nunca é pontual. É sempre com intervalos de confiança
que vão se abrindo ao longo do tempo.
Aqui começam grandes discussões sobre a arte da previsão.
Separação de períodos, validação cruzada, melhoria dos modelos,
etc.
previsao <- forecast(fit_auto, h=24)
previsao
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Aug 2023 151.6521 148.2044 155.0997 146.3794 156.9248
## Sep 2023 148.1628 143.7749 152.5507 141.4521 154.8735
## Oct 2023 148.5976 143.4866 153.7087 140.7809 156.4143
## Nov 2023 147.2329 141.1809 153.2848 137.9773 156.4885
## Dec 2023 147.1065 140.4767 153.7363 136.9671 157.2459
## Jan 2024 141.3483 134.1604 148.5363 130.3553 152.3414
## Feb 2024 143.6630 135.8374 151.4886 131.6948 155.6312
## Mar 2024 154.9047 146.6148 163.1947 142.2263 167.5832
## Apr 2024 148.3618 139.5942 157.1294 134.9530 161.7706
## May 2024 147.8612 138.5952 157.1271 133.6901 162.0322
## Jun 2024 147.5626 137.8906 157.2346 132.7706 162.3546
## Jul 2024 153.0925 142.9987 163.1862 137.6554 168.5296
## Aug 2024 153.4675 142.7226 164.2124 137.0346 169.9004
## Sep 2024 149.8330 138.5652 161.1008 132.6004 167.0656
## Oct 2024 151.1181 139.3299 162.9063 133.0896 169.1466
## Nov 2024 149.4285 137.1014 161.7555 130.5759 168.2811
## Dec 2024 149.0031 136.2113 161.7949 129.4397 168.5665
## Jan 2025 143.6949 130.4335 156.9562 123.4134 163.9764
## Feb 2025 145.4746 131.7454 159.2038 124.4776 166.4717
## Mar 2025 156.1349 141.9813 170.2884 134.4889 177.7808
## Apr 2025 150.0964 135.5145 164.6784 127.7952 172.3976
## May 2025 149.5388 134.5371 164.5406 126.5957 172.4820
## Jun 2025 149.3387 133.9433 164.7341 125.7935 172.8839
## Jul 2025 155.0198 139.2292 170.8105 130.8701 179.1696
Final da Apresentação
Obrigado!
