1º Trabalho Prático
Henrique Drumond Pérez
08/07/2021
- Os dados
- Questão 1 - ACF & PACF
- Questão 2 - Teste de Dickey-Fuller Aumentado
- Questão 3 - Teste de Phillips-Perron PP
- Questão 4 - Relação entre DIV e LC
- Questão 5 - Primeira Diferença das Séries
- Questão 6 - Primeira Diferença DIV x LC
- Questão 7, 8, 9, 10, 11 - Modelo SARIMA(p,d,q)x(P,D,Q)
Os dados
Primeiro, importamos os dados que serão utilizados na resolução dos exercícios propostos e apresentamos o sumário destes(excluindo as colunas dos trimestres e anos).
dados <- read_excel("dados trabalho.xls")
summary(dados[,c(-1,-2)])## RPD PIB DCP LC
## Min. :1073 Min. : 1568 Min. :1017 Min. : 17.50
## 1st Qu.:1899 1st Qu.: 2701 1st Qu.:1698 1st Qu.: 32.42
## Median :3544 Median : 4773 Median :3156 Median : 129.90
## Mean :3946 Mean : 5345 Mean :3591 Mean : 241.00
## 3rd Qu.:5520 3rd Qu.: 7390 3rd Qu.:4966 3rd Qu.: 322.02
## Max. :8695 Max. :11676 Max. :8349 Max. :1441.40
## DIV
## Min. : 6.00
## 1st Qu.: 14.90
## Median : 44.75
## Mean :138.12
## 3rd Qu.:190.20
## Max. :829.40Sabe-se que:
RPD - Renda Real Pessoal Disponível
PIB - Produto Interno Bruto
DCP - Despesas Reais de Consumo Pessoal
LC - Lucros Corporativos
DIV - Dividendos Corporativos
Destaca-se que todas as variáveis são expressas em bilhões de dólares os dados possuem .
Agora, iremos converter os dados em uma série temporal, sabemos que começam no primeiro trimestre do ano de 1947, assim:
dados_ts <- ts(dados[,c(-1,-2)] , start = c(1947,1), frequency = 4)
head(dados_ts)## RPD PIB DCP LC DIV
## 1947 Q1 1096.0 1570.52 1017.2 21.2 6.0
## 1947 Q2 1072.8 1568.65 1034.0 20.0 6.3
## 1947 Q3 1102.8 1567.97 1037.5 19.8 6.5
## 1947 Q4 1089.7 1590.94 1037.7 21.7 6.4
## 1948 Q1 1107.3 1616.07 1042.6 22.9 7.0
## 1948 Q2 1145.3 1644.64 1054.3 24.1 6.7Apresentando os gráficos das séries:
Questão 1 - ACF & PACF
a)Correlogramas
Auto Correlation Function - ACF
Partial Auto Correlation Function
Optei por plotar separadamente os gráficos para melhor visualização.
b)Conclusões
Pela análise dos gráficos, percebe-se que a autocorrelação cai exponencialmente, isto é \(\phi>0\) , no limite tendendo à 0. Além disso, aparentam ser séries não estacionárias e que possuem algum tipo de tendência, ou seja, podem ser Processos Estocásticos em Tendência.
Questão 2 - Teste de Dickey-Fuller Aumentado
Explicações
Define-se a equação de Dickey-Fuler Aumentado como:
\[\Delta r_{t}=\mu+\beta t+(\phi-1)r_{t-1}+\sum_{k=1}^{p-1}\phi_k * \Delta r_{t-k}+\epsilon_t\]
No teste, temos as seguintes hipóteses:
- \(H_0\): possui raiz untária
- \(H_1\): \(r_t\) é estacionária
Utilizamos então a função ur.df para realização da questão, disponível no pacote do R library(urca)
Os valores das estatísticas de teste se referem, respectivamente, ao tau3, phi2 e phi3.
Definimos tau3 como o teste que tem como a hipótese nula de que há raiz unitária, há tendência determinística e há drift.
Já phi2 possui \(H_0\) de que existe raiz unitária e não existe tendência determinística e não existe drift.
Por fim, phi3 tem como \(H_0\) que existe raiz unitária, existe drift e não existe tendência determinística.
Iremos rejeitar \(H_0\):
- tau3 - critical value > test statistic
- phi2 - critical value < test statistic
- phi3 - critical value < test statistic
Adotaremos \(\alpha = 5\%\)
RPD - Renda Real Pessoal Disponível
summary(ur.df(dados_ts[,1], 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
## -142.009 -15.913 0.365 16.190 132.795
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.315849 5.214042 1.787 0.07526 .
## z.lag.1 0.001654 0.006155 0.269 0.78840
## tt 0.181394 0.191056 0.949 0.34337
## z.diff.lag -0.208767 0.063877 -3.268 0.00124 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 39.11 on 238 degrees of freedom
## Multiple R-squared: 0.1399, Adjusted R-squared: 0.129
## F-statistic: 12.9 on 3 and 238 DF, p-value: 7.737e-08
##
##
## Value of test-statistic is: 0.2687 51.3392 18.4812
##
## 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.47Assim, não rejeita-se tau3 porém rejeita-se phi2 e phi3. Conclui-se que \(\phi\) -1 = 0, \(\mu \neq 0\) e \(\beta \neq 0\) e, portanto, a série é caracterizada como uma random walk, com drift e também tendência determinística.
PIB - Produto Interno Bruto
summary(ur.df(dados_ts[,2], 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
## -142.648 -21.628 0.205 23.237 154.839
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.639e+00 5.519e+00 1.384 0.168
## z.lag.1 9.062e-06 4.414e-03 0.002 0.998
## tt 1.907e-01 1.812e-01 1.053 0.293
## z.diff.lag 2.569e-01 6.328e-02 4.060 6.65e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 40.39 on 238 degrees of freedom
## Multiple R-squared: 0.2158, Adjusted R-squared: 0.2059
## F-statistic: 21.83 on 3 and 238 DF, p-value: 1.599e-12
##
##
## Value of test-statistic is: 0.0021 25.4152 10.9829
##
## 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.47Semelhante à primeira conclusão, \(\phi\)-1= 0, \(\mu \neq 0\) e \(\beta \neq 0\), logo é caracterizada como uma random walk, com drift e também tendência determinística.
DCP - Despesas Reais de Consumo Pessoal
summary(ur.df(dados_ts[,3], 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
## -101.080 -12.250 1.033 12.242 66.042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.120986 2.914047 0.728 0.4674
## z.lag.1 0.004145 0.003016 1.375 0.1705
## tt 0.073149 0.087903 0.832 0.4062
## z.diff.lag 0.142789 0.064669 2.208 0.0282 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.39 on 238 degrees of freedom
## Multiple R-squared: 0.3472, Adjusted R-squared: 0.339
## F-statistic: 42.2 on 3 and 238 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: 1.3747 41.3054 29.3665
##
## 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\(\phi\)-1= 0, \(\mu \neq 0\) e \(\beta \neq 0\), logo é caracterizada como uma random walk, com drift e também tendência determinística.
LC - Lucros Corporativos
summary(ur.df(dados_ts[,4], 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
## -92.075 -2.720 0.303 3.281 185.592
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.767422 3.170825 -0.557 0.5778
## z.lag.1 0.018190 0.008034 2.264 0.0245 *
## tt 0.018770 0.033313 0.563 0.5737
## z.diff.lag 0.168061 0.066296 2.535 0.0119 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.72 on 238 degrees of freedom
## Multiple R-squared: 0.1491, Adjusted R-squared: 0.1384
## F-statistic: 13.9 on 3 and 238 DF, p-value: 2.204e-08
##
##
## Value of test-statistic is: 2.2642 9.0793 9.4205
##
## 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\(\phi\)-1 = 0, \(\mu \neq 0\) e \(\beta \neq 0\), logo é caracterizada como uma random walk, com drift e também tendência determinística.
DIV - Dividendos Corporativos
summary(ur.df(dados_ts[,4], 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
## -92.075 -2.720 0.303 3.281 185.592
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.767422 3.170825 -0.557 0.5778
## z.lag.1 0.018190 0.008034 2.264 0.0245 *
## tt 0.018770 0.033313 0.563 0.5737
## z.diff.lag 0.168061 0.066296 2.535 0.0119 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.72 on 238 degrees of freedom
## Multiple R-squared: 0.1491, Adjusted R-squared: 0.1384
## F-statistic: 13.9 on 3 and 238 DF, p-value: 2.204e-08
##
##
## Value of test-statistic is: 2.2642 9.0793 9.4205
##
## 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\(\phi\) -1 = 0, \(\mu \neq 0\) e \(\beta \neq 0\), logo é caracterizada como uma random walk, com drift e também tendência determinística.
Questão 3 - Teste de Phillips-Perron PP
Problema do ADF é de que considera que \(\epsilon_t\) é um Ruído Branco, ou seja, de que possuí variância constante e não apresenta o problema da autocorrelação. Isto é corrigido pelo teste de Phillips-Perron: \[\Delta r_t = \mu + \beta_t + \rho r_{t-1} + \epsilon_t \] Embora seja um teste mais simples, é mais robusto do que o teste ADF.
RPD - Renda Real Pessoal Disponível
pp.test(dados_ts[,1], type = "Z(t_alpha)")## Warning in pp.test(dados_ts[, 1], type = "Z(t_alpha)"): p-value greater than
## printed p-value##
## Phillips-Perron Unit Root Test
##
## data: dados_ts[, 1]
## Dickey-Fuller Z(t_alpha) = 0.018707, Truncation lag parameter = 4,
## p-value = 0.99
## alternative hypothesis: stationaryur.pp(dados_ts[,1], type="Z-tau", model="trend", lags="short")##
## ##################################################
## # Phillips-Perron Unit Root / Cointegration Test #
## ##################################################
##
## The value of the test statistic is: 0.018## 1pct 5pct 10pct
## critical values -3.998978 -3.429523 -3.137979Como o valor da estatística de teste é superior ao valor crítico com 1% e com 5% de significância, não podemos rejeitar a \(H_0\) de não estacionaridade.
PIB - Produto Interno Bruto
pp.test(dados_ts[,2], type = "Z(t_alpha)")## Warning in pp.test(dados_ts[, 2], type = "Z(t_alpha)"): p-value greater than
## printed p-value##
## Phillips-Perron Unit Root Test
##
## data: dados_ts[, 2]
## Dickey-Fuller Z(t_alpha) = 0.011409, Truncation lag parameter = 4,
## p-value = 0.99
## alternative hypothesis: stationaryur.pp(dados_ts[,2], type="Z-tau", model="trend", lags="short")##
## ##################################################
## # Phillips-Perron Unit Root / Cointegration Test #
## ##################################################
##
## The value of the test statistic is: 0.0126## 1pct 5pct 10pct
## critical values -3.998978 -3.429523 -3.137979Como o valor da estatística de teste é superior ao valor crítico com 1% e com 5% de significância, não podemos rejeitar a \(H_0\) de não estacionaridade.
DCP - Despesas Reais de Consumo Pessoal
pp.test(dados_ts[,3], type = "Z(t_alpha)")## Warning in pp.test(dados_ts[, 3], type = "Z(t_alpha)"): p-value greater than
## printed p-value##
## Phillips-Perron Unit Root Test
##
## data: dados_ts[, 3]
## Dickey-Fuller Z(t_alpha) = 1.2263, Truncation lag parameter = 4,
## p-value = 0.99
## alternative hypothesis: stationaryur.pp(dados_ts[,3], type="Z-tau", model="trend", lags="short")##
## ##################################################
## # Phillips-Perron Unit Root / Cointegration Test #
## ##################################################
##
## The value of the test statistic is: 1.2275## 1pct 5pct 10pct
## critical values -3.998978 -3.429523 -3.137979Como o valor da estatística de teste é superior ao valor crítico com 1% e com 5% de significância, não podemos rejeitar a \(H_0\) de não estacionaridade.
LC - Lucros Corporativos
pp.test(dados_ts[,4], type = "Z(t_alpha)")## Warning in pp.test(dados_ts[, 4], type = "Z(t_alpha)"): p-value greater than
## printed p-value##
## Phillips-Perron Unit Root Test
##
## data: dados_ts[, 4]
## Dickey-Fuller Z(t_alpha) = 2.9057, Truncation lag parameter = 4,
## p-value = 0.99
## alternative hypothesis: stationaryur.pp(dados_ts[,4], type="Z-tau", model="trend", lags="short")##
## ##################################################
## # Phillips-Perron Unit Root / Cointegration Test #
## ##################################################
##
## The value of the test statistic is: 2.9097## 1pct 5pct 10pct
## critical values -3.998978 -3.429523 -3.137979Como o valor da estatística de teste é superior ao valor crítico com 1% e com 5% de significância, não podemos rejeitar a \(H_0\) de não estacionaridade.
DIV - Dividendos Corporativos
pp.test(dados_ts[,5], type = "Z(t_alpha)")## Warning in pp.test(dados_ts[, 5], type = "Z(t_alpha)"): p-value greater than
## printed p-value##
## Phillips-Perron Unit Root Test
##
## data: dados_ts[, 5]
## Dickey-Fuller Z(t_alpha) = 5.6917, Truncation lag parameter = 4,
## p-value = 0.99
## alternative hypothesis: stationaryur.pp(dados_ts[,5], type="Z-tau", model="trend", lags="short")##
## ##################################################
## # Phillips-Perron Unit Root / Cointegration Test #
## ##################################################
##
## The value of the test statistic is: 5.6776## 1pct 5pct 10pct
## critical values -3.998978 -3.429523 -3.137979Como o valor da estatística de teste é superior ao valor crítico com 1% e com 5% de significância, não podemos rejeitar a \(H_0\) de não estacionaridade.
Conclusão
Tanto os testes de ADF quanto PP apontam para a não estacionaridade das séries analisadas. No caso de que apresentem resultados contraditórios, deve-se tentar corrigir a série de forma a apresentarem resultados semelhantes, mas, caso isso não seja possível, opta-se pelo teste de PP pois é mais robusto.
Questão 4 - Relação entre DIV e LC
Não esperaria a ocorrência de correlação espúria pois para a realização de análises econométricas via regressões lineares é necessária a forte intuição econômica por trás das variáveis analisadas. Como é amplamente descrito na vasta literatura de finanças, há uma relação direta entre dividendos corporativos e o lucro das companhias.
Questão 5 - Primeira Diferença das Séries
Tirando a primeira diferença:
dados_ts_dif <- apply(dados[,c(-1,-2)], 2, diff)
dados_ts_dif <- ts(dados_ts_dif , start = c(1947,2), frequency = 4)a) Gráficos
Aparentemente, são estacionárias, oscilando ao redor de um determinado ponto.
b) Correlogramas
Auto Correlation Function - ACF
Partial Auto Correlation Function
Pelos correlogramas das variações, vemos mais observações significativas ao longo das amostras, dando indícios de que um modelo ARMA seria razoável paras séries.
Questão 6 - Primeira Diferença DIV x LC
a) A regressão é espúria?
Não, uma vez que há, muitas vezes, uma relação direta entre acréscimos e decréscimos entre dividendos e lucros das ações. Assim, é de se esperar que caso haja um aumento significativo dos lucros também haverá um aumento significativo dos dividendos (considerando que a empresa não investirá em outros projetos).
b) Deve-se incluir o intercepto?
Ao removermos o intercepto, obrigaríamos a reta de regressão a passar pela origem, o que pode piorar o ajuste desta. Além disso, o intercepto seria o ponto em que os lucros são zero, necessariamente obrigando a distribuição de dividendos à ser zero, o que muitas vezes pode não ser verdade. Caso alguma empresa possua reservas de lucros significativas e perceba que o “não lucro” seja devido à, por exemplo, um efeito de sazonalidade, ela poderá distribuir seus lucros retidos.
c) Estimativas & Resultados
Analisando a primeira diferença
A primeira diferença:
head(dados_ts_dif)## RPD PIB DCP LC DIV
## 1947 Q2 -23.2 -1.87 16.8 -1.2 0.3
## 1947 Q3 30.0 -0.68 3.5 -0.2 0.2
## 1947 Q4 -13.1 22.97 0.2 1.9 -0.1
## 1948 Q1 17.6 25.13 4.9 1.2 0.6
## 1948 Q2 38.0 28.57 11.7 1.2 -0.3
## 1948 Q3 23.1 9.42 1.8 -0.3 0.4Vemos que queremos rodar a regressão com a quarta coluna explicando a quinta coluna:
summary(lm(dados_ts_dif[,5] ~ dados_ts_dif[,4]))##
## Call:
## lm(formula = dados_ts_dif[, 5] ~ dados_ts_dif[, 4])
##
## Residuals:
## Min 1Q Median 3Q Max
## -90.144 -3.326 -2.718 -0.286 123.416
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.444886 0.789931 4.361 1.92e-05 ***
## dados_ts_dif[, 4] -0.009761 0.032895 -0.297 0.767
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.95 on 241 degrees of freedom
## Multiple R-squared: 0.0003652, Adjusted R-squared: -0.003783
## F-statistic: 0.08805 on 1 and 241 DF, p-value: 0.7669Quando analisamos a diferença, a variação nos lucros corporativos foi estatísticamente igual à 0.
Analisando os dados sem diferença
head(dados_ts)## RPD PIB DCP LC DIV
## 1947 Q1 1096.0 1570.52 1017.2 21.2 6.0
## 1947 Q2 1072.8 1568.65 1034.0 20.0 6.3
## 1947 Q3 1102.8 1567.97 1037.5 19.8 6.5
## 1947 Q4 1089.7 1590.94 1037.7 21.7 6.4
## 1948 Q1 1107.3 1616.07 1042.6 22.9 7.0
## 1948 Q2 1145.3 1644.64 1054.3 24.1 6.7Vemos que queremos rodar a regressão com a quarta coluna explicando a quinta coluna:
summary(lm(dados_ts[,5] ~ dados_ts[,4]))##
## Call:
## lm(formula = dados_ts[, 5] ~ dados_ts[, 4])
##
## Residuals:
## Min 1Q Median 3Q Max
## -110.103 -9.092 -3.645 6.106 105.791
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.722594 2.537779 -0.285 0.776
## dados_ts[, 4] 0.576118 0.006409 89.892 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31.45 on 242 degrees of freedom
## Multiple R-squared: 0.9709, Adjusted R-squared: 0.9708
## F-statistic: 8081 on 1 and 242 DF, p-value: < 2.2e-16Percebe-se que há uma relação entre os lucros corporativos e o aumento no pagamento de dividendos, com R Quadrado ajustado bastante alto. Assim, variações pontuais em um único ano não explicam bem os dividendos daquele ano em específico, porém, ao longo do tempo, há uma relação bastante clara.
d) Testar estacionaridade dos resíduos da regressão
reg <- lm(dados_ts[,5] ~ dados_ts[,4])
summary(ur.df(reg$residuals, 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
## -187.415 -3.407 0.187 1.846 96.070
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.16368 2.31435 -0.935 0.3508
## z.lag.1 -0.15739 0.03840 -4.099 5.7e-05 ***
## tt 0.01810 0.01655 1.094 0.2752
## z.diff.lag -0.11539 0.06446 -1.790 0.0747 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.35 on 238 degrees of freedom
## Multiple R-squared: 0.1011, Adjusted R-squared: 0.0898
## F-statistic: 8.926 on 3 and 238 DF, p-value: 1.256e-05
##
##
## Value of test-statistic is: -4.0987 5.6006 8.3997
##
## 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.47pp.test(reg$residuals, type = "Z(t_alpha)")## Warning in pp.test(reg$residuals, type = "Z(t_alpha)"): p-value smaller than
## printed p-value##
## Phillips-Perron Unit Root Test
##
## data: reg$residuals
## Dickey-Fuller Z(t_alpha) = -4.4522, Truncation lag parameter = 4,
## p-value = 0.01
## alternative hypothesis: stationaryComo comprovado, deve-se rejeitar \(H_0\), então devemos considerar os resíduos como estacionários.
Questão 7, 8, 9, 10, 11 - Modelo SARIMA(p,d,q)x(P,D,Q)
RPD - Renda Real Pessoal Disponível
PIB - Produto Interno Bruto
DCP - Despesas Reais de Consumo Pessoal
LC - Lucros Corporativos
DIV - Dividendos Corporativos
RPD - Renda Real Pessoal Disponível
a) Passos & Conclusão
1º Analisar gráfico da série 2º Verificar correlogramas da série 3º Verificar se é estacionária 4º Verificar quais modelos possuem menor IC
Divide-se a amostra em duas, a primeira destinada à estimação dos parâmetros e a segunda destinada ao teste dos parâmetros para verificar se estão adequados para a realização de previsões e inferências.
Tira-se da amostra de estimação os 4 últimos trimestres:
dados_ts_estim <- ts((dados_ts[1:((nrow(dados_ts)/2)-2),]), start = c(1947,1), frequency = 4)
tail(dados_ts_estim)## RPD PIB DCP LC DIV
## 1975 Q3 3309.1 4340.87 2901.2 102.3 32.9
## 1975 Q4 3342.0 4397.81 2931.4 107.9 33.4
## 1976 Q1 3390.9 4496.76 2989.7 113.2 36.2
## 1976 Q2 3417.5 4530.33 3016.3 114.3 38.1
## 1976 Q3 3448.0 4552.03 3047.9 115.1 39.9
## 1976 Q4 3473.0 4584.62 3088.0 115.1 41.9dados_ts_test <- ts((dados_ts[-(1:((nrow(dados_ts)/2)))-4,]), start = c(1976,3), frequency = 4)
tail(dados_ts_test) ## RPD PIB DCP LC DIV
## 2005 Q3 8384.5 11336.7 8063.8 1381.0 711.1
## 2005 Q4 8510.7 11395.5 8141.2 1336.8 736.4
## 2006 Q1 8623.9 11412.6 8215.7 1363.3 759.4
## 2006 Q2 8607.1 11520.1 8244.3 1441.4 784.2
## 2006 Q3 8692.1 11658.9 8302.2 1410.2 807.7
## 2006 Q4 8695.2 11675.7 8349.1 1425.5 829.4dados_ts_fora_amostra <- ts(dados_ts[-(1:(nrow(dados_ts)-4)),], start = c(2007, 1), frequency = 4)
dados_ts_fora_amostra## RPD PIB DCP LC DIV
## 2007 Q1 8623.9 11412.6 8215.7 1363.3 759.4
## 2007 Q2 8607.1 11520.1 8244.3 1441.4 784.2
## 2007 Q3 8692.1 11658.9 8302.2 1410.2 807.7
## 2007 Q4 8695.2 11675.7 8349.1 1425.5 829.4b) Selecionando os modelos pelo IC & MAPE
Função no R para selecionar os modelos, com base no AIC
SARIMA_RPD <- auto.arima(dados_ts_estim[,1], ic = "aicc",seasonal = TRUE, trace = TRUE, stepwise = FALSE, approximation = FALSE)##
## ARIMA(0,2,0) : 1177.739
## ARIMA(0,2,0)(0,0,1)[4] : 1178.493
## ARIMA(0,2,0)(0,0,2)[4] : 1179.231
## ARIMA(0,2,0)(1,0,0)[4] : 1178.713
## ARIMA(0,2,0)(1,0,1)[4] : 1179.274
## ARIMA(0,2,0)(1,0,2)[4] : 1181.071
## ARIMA(0,2,0)(2,0,0)[4] : 1179.444
## ARIMA(0,2,0)(2,0,1)[4] : 1182.308
## ARIMA(0,2,0)(2,0,2)[4] : 1183.13
## ARIMA(0,2,1) : 1104.381
## ARIMA(0,2,1)(0,0,1)[4] : 1103.141
## ARIMA(0,2,1)(0,0,2)[4] : 1102.07
## ARIMA(0,2,1)(1,0,0)[4] : 1104.138
## ARIMA(0,2,1)(1,0,1)[4] : 1103.079
## ARIMA(0,2,1)(1,0,2)[4] : 1102.922
## ARIMA(0,2,1)(2,0,0)[4] : 1102.705
## ARIMA(0,2,1)(2,0,1)[4] : 1104.044
## ARIMA(0,2,1)(2,0,2)[4] : Inf
## ARIMA(0,2,2) : 1106.462
## ARIMA(0,2,2)(0,0,1)[4] : 1104.84
## ARIMA(0,2,2)(0,0,2)[4] : 1102.783
## ARIMA(0,2,2)(1,0,0)[4] : 1106.115
## ARIMA(0,2,2)(1,0,1)[4] : 1104.072
## ARIMA(0,2,2)(1,0,2)[4] : Inf
## ARIMA(0,2,2)(2,0,0)[4] : 1104.035
## ARIMA(0,2,2)(2,0,1)[4] : 1104.816
## ARIMA(0,2,3) : 1108.543
## ARIMA(0,2,3)(0,0,1)[4] : 1106.842
## ARIMA(0,2,3)(0,0,2)[4] : 1103.411
## ARIMA(0,2,3)(1,0,0)[4] : 1108.24
## ARIMA(0,2,3)(1,0,1)[4] : 1105.456
## ARIMA(0,2,3)(2,0,0)[4] : 1105.272
## ARIMA(1,2,0) : 1144.451
## ARIMA(1,2,0)(0,0,1)[4] : 1140.135
## ARIMA(1,2,0)(0,0,2)[4] : 1129.206
## ARIMA(1,2,0)(1,0,0)[4] : 1144.032
## ARIMA(1,2,0)(1,0,1)[4] : 1134.081
## ARIMA(1,2,0)(1,0,2)[4] : 1131.313
## ARIMA(1,2,0)(2,0,0)[4] : 1133.956
## ARIMA(1,2,0)(2,0,1)[4] : 1127.225
## ARIMA(1,2,0)(2,0,2)[4] : 1128.01
## ARIMA(1,2,1) : 1106.461
## ARIMA(1,2,1)(0,0,1)[4] : 1104.806
## ARIMA(1,2,1)(0,0,2)[4] : 1102.337
## ARIMA(1,2,1)(1,0,0)[4] : 1106.109
## ARIMA(1,2,1)(1,0,1)[4] : 1103.821
## ARIMA(1,2,1)(1,0,2)[4] : 1104.548
## ARIMA(1,2,1)(2,0,0)[4] : 1103.875
## ARIMA(1,2,1)(2,0,1)[4] : 1104.291
## ARIMA(1,2,2) : 1108.604
## ARIMA(1,2,2)(0,0,1)[4] : 1103.159
## ARIMA(1,2,2)(0,0,2)[4] : 1099.005
## ARIMA(1,2,2)(1,0,0)[4] : Inf
## ARIMA(1,2,2)(1,0,1)[4] : 1100.783
## ARIMA(1,2,2)(2,0,0)[4] : 1103.03
## ARIMA(1,2,3) : 1110.607
## ARIMA(1,2,3)(0,0,1)[4] : 1105.379
## ARIMA(1,2,3)(1,0,0)[4] : 1110.173
## ARIMA(2,2,0) : 1125.312
## ARIMA(2,2,0)(0,0,1)[4] : 1124.298
## ARIMA(2,2,0)(0,0,2)[4] : 1120.301
## ARIMA(2,2,0)(1,0,0)[4] : 1125.178
## ARIMA(2,2,0)(1,0,1)[4] : 1122.449
## ARIMA(2,2,0)(1,0,2)[4] : 1122.504
## ARIMA(2,2,0)(2,0,0)[4] : 1123.221
## ARIMA(2,2,0)(2,0,1)[4] : Inf
## ARIMA(2,2,1) : 1108.579
## ARIMA(2,2,1)(0,0,1)[4] : 1106.931
## ARIMA(2,2,1)(0,0,2)[4] : 1103.592
## ARIMA(2,2,1)(1,0,0)[4] : 1108.278
## ARIMA(2,2,1)(1,0,1)[4] : 1105.616
## ARIMA(2,2,1)(2,0,0)[4] : 1105.54
## ARIMA(2,2,2) : 1110.54
## ARIMA(2,2,2)(0,0,1)[4] : 1105.379
## ARIMA(2,2,2)(1,0,0)[4] : 1110.477
## ARIMA(2,2,3) : Inf
## ARIMA(3,2,0) : 1126.417
## ARIMA(3,2,0)(0,0,1)[4] : 1097.927
## ARIMA(3,2,0)(0,0,2)[4] : 1099.674
## ARIMA(3,2,0)(1,0,0)[4] : 1118.655
## ARIMA(3,2,0)(1,0,1)[4] : 1099.866
## ARIMA(3,2,0)(2,0,0)[4] : 1110.886
## ARIMA(3,2,1) : 1109.005
## ARIMA(3,2,1)(0,0,1)[4] : 1108.198
## ARIMA(3,2,1)(1,0,0)[4] : 1109.036
## ARIMA(3,2,2) : 1110.74
##
##
##
## Best model: ARIMA(3,2,0)(0,0,1)[4]Três modelos com menor IC & MAPE:
- SARIMA(3,2,0)(0,0,1) - IC 1097.927
arima_rpd <- Arima(dados_ts_estim[,1], order = c(3,2,0), seasonal = c(0,0,1))
forecast_rpd <- forecast(arima_rpd,122)
MAPE(forecast_rpd$mean, dados_ts_test[,1])## [1] 0.161823- SARIMA(1,2,2)(0,0,2) - IC 1099.005
arima_rpd <- Arima(dados_ts_estim[,1], order = c(1,2,2), seasonal = c(0,0,2))
forecast_rpd <- forecast(arima_rpd,122)
MAPE(forecast_rpd$mean, dados_ts_test[,1])## [1] 0.166092- SARIMA(3,2,0)(0,0,2) - IC 1099.674
arima_rpd <- Arima(dados_ts_estim[,1], order = c(3,2,0), seasonal = c(0,0,2))
forecast_rpd <- forecast(arima_rpd,122)
MAPE(forecast_rpd$mean, dados_ts_test[,1])## [1] 0.1640344O modelo que possui o menor MAPE é o SARIMA(3,2,0)(0,0,1).
c) Apresentando fora da amostra
arima_rpd <- Arima(dados_ts_test[,1], order = c(3,2,0), seasonal = c(0,0,1))
forecast(arima_rpd,4)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2007 Q1 8759.837 8442.136 9077.538 8273.955 9245.719
## 2007 Q2 8807.390 8351.663 9263.118 8110.415 9504.366
## 2007 Q3 8867.547 8302.067 9433.026 8002.720 9732.373
## 2007 Q4 8922.535 8266.297 9578.772 7918.906 9926.163dados_ts_fora_amostra[,1]## Qtr1 Qtr2 Qtr3 Qtr4
## 2007 8623.9 8607.1 8692.1 8695.2Percebe-se que há um erro significativo de previsão do modelo.
Para os demais dados não será considerada a questão a) pois seria redundante. Também não serão mostrados os modelos SARIMA devido a extensão destes na apresentação.
PIB
b) Selecionando os modelos pelo IC & MAPE
Função no R para selecionar os modelos, com base no AIC, todas contém drift
Três modelos com menor IC & MAPE:
- SARIMA(0,1,2)(0,0,0) - IC 1158.587
arima_pib <- Arima(dados_ts_estim[,2], order = c(0,1,2), seasonal = c(0,0,0), include.drift = TRUE )
forecast_pib <- forecast(arima_pib,122)
MAPE(forecast_pib$mean, dados_ts_test[,2])## [1] 0.2248325- SARIMA(3,1,1)(0,0,0) - IC 1158.922
arima_pib <- Arima(dados_ts_estim[,2], order = c(3,1,1), seasonal = c(0,0,0), include.drift = TRUE )
forecast_pib <- forecast(arima_pib,122)
MAPE(forecast_pib$mean, dados_ts_test[,2])## [1] 0.2288261- SARIMA(2,1,1)(0,0,0) - IC 1159.776
arima_pib <- Arima(dados_ts_estim[,2], order = c(2,1,1), seasonal = c(0,0,0), include.drift = TRUE )
forecast_pib <- forecast(arima_pib,122)
MAPE(forecast_pib$mean, dados_ts_test[,2])## [1] 0.2245943O modelo que possui o menor MAPE é o SARIMA(2,1,1)(0,0,0)
c) Apresentando fora da amostra
arima_pib <- Arima(dados_ts_test[,2], order = c(2,1,1), seasonal = c(0,0,0), include.drift = TRUE )
forecast(arima_pib,4)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2007 Q1 11742.92 11342.08 12143.76 11129.88 12355.96
## 2007 Q2 11814.14 11259.77 12368.51 10966.30 12661.98
## 2007 Q3 11885.28 11237.44 12533.11 10894.50 12876.06
## 2007 Q4 11956.04 11245.74 12666.34 10869.72 13042.35dados_ts_fora_amostra[,2]## Qtr1 Qtr2 Qtr3 Qtr4
## 2007 11412.6 11520.1 11658.9 11675.7Como no exemplo passado, o modelo superestimou os valores.
DCP
b) Selecionando os modelos pelo IC & MAPE
Função no R para selecionar os modelos, com base no AIC
Três modelos com menor IC & MAPE:
- SARIMA(2,2,1)(2,0,0) - IC 1007.858
arima_dcp <- Arima(dados_ts_estim[,3], order = c(2,2,1), seasonal = c(2,0,0))
forecast_dcp <- forecast(arima_dcp,122)
MAPE(forecast_dcp$mean, dados_ts_test[,3])## [1] 0.1793822- SARIMA(0,2,3)(2,0,0) - IC 1007.904
arima_dcp <- Arima(dados_ts_estim[,3], order = c(0,2,3), seasonal = c(2,0,0))
forecast_dcp <- forecast(arima_dcp,122)
MAPE(forecast_dcp$mean, dados_ts_test[,3])## [1] 0.1800532- SARIMA(2,2,1)(0,0,2) - IC 1010.039
arima_dcp <- Arima(dados_ts_estim[,3], order = c(2,2,1), seasonal = c(0,0,2))
forecast_dcp <- forecast(arima_dcp,122)
MAPE(forecast_dcp$mean, dados_ts_test[,3])## [1] 0.1710023O modelo que possui o menor MAPE é o SARIMA(2,2,1)(0,0,2)
c) Apresentando fora da amostra
arima_dcp <- Arima(dados_ts_test[,3], order = c(2,2,1), seasonal = c(0,0,2))
forecast(arima_dcp,4)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2007 Q1 8420.334 8149.884 8690.784 8006.716 8833.952
## 2007 Q2 8475.212 8093.166 8857.258 7890.924 9059.501
## 2007 Q3 8536.943 8069.874 9004.012 7822.622 9251.264
## 2007 Q4 8602.055 8061.707 9142.402 7775.664 9428.445dados_ts_fora_amostra[,3]## Qtr1 Qtr2 Qtr3 Qtr4
## 2007 8215.7 8244.3 8302.2 8349.1Modelo segue a superestimar os valores.
LC
b) Selecionando os modelos pelo IC & MAPE
Função no R para selecionar os modelos, com base no AIC
Três modelos com menor IC & MAPE:
- SARIMA(3,2,2)(0,0,0) - IC 601.4969
arima_lc <- Arima(dados_ts_estim[,4], order = c(3,2,2), seasonal = c(0,0,0))
forecast_lc <- forecast(arima_lc,122)
MAPE(forecast_lc$mean, dados_ts_test[,4])## [1] 0.4081371- SARIMA(3,2,1)(0,0,0) - IC 601.7582
arima_lc <- Arima(dados_ts_estim[,4], order = c(3,2,1), seasonal = c(0,0,0))
forecast_lc <- forecast(arima_lc,122)
MAPE(forecast_lc$mean, dados_ts_test[,4])## [1] 0.4012832- SARIMA(2,2,3)(0,0,0) - IC 601.8721
arima_lc <- Arima(dados_ts_estim[,4], order = c(2,2,3), seasonal = c(0,0,0))
forecast_lc <- forecast(arima_lc,122)
MAPE(forecast_lc$mean, dados_ts_test[,4])## [1] 0.4438252O modelo que possui o menor MAPE é o SARIMA(3,2,1)(0,0,0)
c) Apresentando fora da amostra
arima_dcp <- Arima(dados_ts_test[,4], order = c(3,2,1), seasonal = c(0,0,0))
forecast(arima_dcp,4)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2007 Q1 1461.787 1420.060 1503.514 1397.971 1525.603
## 2007 Q2 1488.802 1422.697 1554.907 1387.703 1589.900
## 2007 Q3 1502.735 1421.525 1583.945 1378.535 1626.934
## 2007 Q4 1519.668 1422.517 1616.819 1371.089 1668.247dados_ts_fora_amostra[,4]## Qtr1 Qtr2 Qtr3 Qtr4
## 2007 1363.3 1441.4 1410.2 1425.5Embora tenha superestimado novamente, Lo80 ficou bastante próximo dos valores que realmente ocorreram.
DIV
b) Selecionando os modelos pelo IC & MAPE
Função no R para selecionar os modelos, com base no AIC
Três modelos com menor IC & MAPE:
- SARIMA(0,2,3)(0,0,0) - IC 176.0856
arima_div <- Arima(dados_ts_estim[,5], order = c(0,2,3), seasonal = c(0,0,0))
forecast_div <- forecast(arima_div,122)
MAPE(forecast_div$mean, dados_ts_test[,5])## [1] 0.5907032- SARIMA(2,2,1)(0,0,0) - IC 176.6826
arima_div <- Arima(dados_ts_estim[,5], order = c(2,2,1), seasonal = c(0,0,0))
forecast_div <- forecast(arima_div,122)
MAPE(forecast_div$mean, dados_ts_test[,5])## [1] 0.6008486- SARIMA(1,2,3)(0,0,0) - IC 177.7502
arima_div <- Arima(dados_ts_estim[,5], order = c(1,2,3), seasonal = c(0,0,0))
forecast_div <- forecast(arima_div,122)
MAPE(forecast_div$mean, dados_ts_test[,5])## [1] 0.6063996O modelo que possui o menor MAPE é o SARIMA(0,2,3)(0,0,0)
c) Apresentando fora da amostra
arima_div <- Arima(dados_ts_test[,5], order = c(0,2,3), seasonal = c(0,0,0))
forecast(arima_div,5)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2007 Q1 851.3707 830.8929 871.8486 820.0526 882.6889
## 2007 Q2 873.4890 847.0305 899.9475 833.0243 913.9538
## 2007 Q3 895.6075 863.1300 928.0849 845.9374 945.2775
## 2007 Q4 917.7259 879.0934 956.3584 858.6426 976.8092
## 2008 Q1 939.8443 894.8798 984.8088 871.0770 1008.6116dados_ts_fora_amostra[,5]## Qtr1 Qtr2 Qtr3 Qtr4
## 2007 759.4 784.2 807.7 829.4Como os demais, apresentam desvios relativamente grandes.