Atividade Prática 4

Introdução e Tratamento dos Dados

Usaremos o conjunto de pacotes "tidyverse" para realizar algumas operações, sobretudo utilizando recursos gráficos.

##### ----- Importando Pacotes ------ #####

library("tidyverse")

## Warning: package 'ggplot2' was built under R version 4.0.5

## Warning: package 'readr' was built under R version 4.0.4

library("tidyquant")
library("Quandl")
library("timetk")

## Warning: package 'timetk' was built under R version 4.0.5

library("devtools")
library("readxl")
library("AER")
library("fBasics")

Importação e Construção da Base de Dados:

##### ------ Importando Dados ------ #####

# Cotações
ativos <- c("VALE3.SA", "^BVSP") # Localizar e Substituir o nome do ativo por outros
inicio <- as.Date("2010-01-04")
fim <- as.Date("2020-12-30")

cotacao <- ativos %>% 
    tq_get(get = "stock.prices", from = inicio, to = fim) %>% 
    mutate(adjusted = na.locf(adjusted))

# Retornos
df_retornos <- cotacao %>% 
    group_by(symbol) %>% 
    select(date, adjusted) %>% 
    tq_transmute(select = adjusted,
                 mutate_fun = periodReturn,
                 period = "daily",
                 type = "arithmetic",
                 col_rename = "retorno") %>% 
    spread(key = symbol, value = retorno)

##### ------ Importando a taxa SELIC da API do BACEN (de 2010 até hoje) ------ #####

Quandl.api_key('TC1ow5j6G7s4SFHTzgDz') # Acrescentando a chave da API - Comando necessário pra acessar o Quandl
selic_diaria <-  as.data.frame(Quandl("BCB/11", type = "xts")/100) # Importando a serie do selic do Bacen
df_selic <- selic_diaria %>% 
    rownames_to_column() %>% 
    tibble() %>% 
    rename("date" = rowname, "selic" = V1) %>% 
    mutate(selic = na.locf(selic))

df_selic$date <- as.Date(df_selic$date)

# acréscimo na base de dados
df_retornos <- df_retornos %>% left_join(df_selic, by = "date")

Estipulando modelos de regressão:

##### ------ Montando regressão do CAPM ------ #####

df_capm <- df_retornos %>% 
    transmute("date" = date,
              "Rp-Rf" = VALE3.SA - selic,
              "Rm-Rf" = `^BVSP` - selic)

capm1 <- lm(df_capm$`Rp-Rf` ~ df_capm$`Rm-Rf`)
summary(capm1)

## 
## Call:
## lm(formula = df_capm$`Rp-Rf` ~ df_capm$`Rm-Rf`)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.221786 -0.010506 -0.000477  0.010061  0.114318 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     0.0003861  0.0003983   0.969    0.332    
## df_capm$`Rm-Rf` 1.0576230  0.0250455  42.228   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02078 on 2719 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.3961, Adjusted R-squared:  0.3959 
## F-statistic:  1783 on 1 and 2719 DF,  p-value: < 2.2e-16

##### ------ Montando regressão de Fama-French ------ #####

MKT_Factor <- read_excel("D:/Dropbox/UFMG/Programação/R/Métodos/Ativ Prát 3/Market_Factor.xls")

## New names:
## * `` -> ...6

HML_Factor <- read_excel("D:/Dropbox/UFMG/Programação/R/Métodos/Ativ Prát 3/HML_Factor.xls")

## New names:
## * `` -> ...6

SMB_Factor <- read_excel("D:/Dropbox/UFMG/Programação/R/Métodos/Ativ Prát 3/SMB_Factor.xls")

## New names:
## * `` -> ...6

MKT <- MKT_Factor %>% select(date, Rm_minus_Rf) %>% rename("MKT" = Rm_minus_Rf)
MKT$date <- as.Date(MKT$date)

HML <- HML_Factor %>% select(date, HML)
HML$date <- as.Date(HML$date)

SMB <- SMB_Factor %>% select(date, SMB)
SMB$date <- as.Date(SMB$date)

fatores <- df_capm %>% 
    left_join(HML, by = "date") %>% 
    left_join(SMB, by = "date") %>% 
    left_join(MKT, by = "date") %>% 
    select(-`Rm-Rf`)

ffrench1 <- lm(`Rp-Rf` ~ MKT + HML + SMB, data = fatores)
summary(ffrench1)

## 
## Call:
## lm(formula = `Rp-Rf` ~ MKT + HML + SMB, data = fatores)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.230730 -0.010448 -0.000378  0.009950  0.133562 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0002042  0.0003985   0.512    0.608    
## MKT          0.9891863  0.0294807  33.554  < 2e-16 ***
## HML          0.7748240  0.0551971  14.037  < 2e-16 ***
## SMB         -0.4052043  0.0507856  -7.979 2.16e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02077 on 2714 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.3976, Adjusted R-squared:  0.3969 
## F-statistic: 597.1 on 3 and 2714 DF,  p-value: < 2.2e-16

QUESTÃO 1:

# A) Teste de normalidade dos resíduos de Jarque-Bera

jarqueberaTest(residuals(capm1))

## 
## Title:
##  Jarque - Bera Normalality Test
## 
## Test Results:
##   STATISTIC:
##     X-squared: 6999.3136
##   P VALUE:
##     Asymptotic p Value: < 2.2e-16 
## 
## Description:
##  Mon Jul 12 22:59:18 2021 by user: Bruno Marcelino

jarqueberaTest(residuals(ffrench1))

## 
## Title:
##  Jarque - Bera Normalality Test
## 
## Test Results:
##   STATISTIC:
##     X-squared: 9779.084
##   P VALUE:
##     Asymptotic p Value: < 2.2e-16 
## 
## Description:
##  Mon Jul 12 22:59:18 2021 by user: Bruno Marcelino

# tanto para o modelo de Fama-French quanto para o CAPM, o teste rejeita a hipótese nula de que os resíduos seguem uma distribuição normal.

# B) Medidas de Influência das observações

outliers_capm <- influence.measures(capm1)
outliers_ffrench <- influence.measures(ffrench1)

# Gráficos de identificação dos outliers

par(mfrow=c(2,3))
plot(ffrench1, which=1:6)

# para o modelo de Fama-French, o principal gráfico (influência em relação a resíduos padronizados) indica que não há outliers tão influentes sobre o modelo, o que pode caracterizar sua robustez (pois envolve mais variáveis)

par(mfrow=c(2,3))
plot(capm1, which=1:6)

# para o modelo CAPM, houveram diversos resíduos que ultrapassaram o limite padrão estabelecido para a influência, interferindo no ajuste do modelo.

# C) Testes de heterocedasticidade de Breusch-Pagan

bptest(capm1)

## 
##  studentized Breusch-Pagan test
## 
## data:  capm1
## BP = 0.79666, df = 1, p-value = 0.3721

# o modelo CAPM aceita a hipótese nula do teste de Breusch-Pagan, o que o caracteriza como homocedástico

bptest(ffrench1)

## 
##  studentized Breusch-Pagan test
## 
## data:  ffrench1
## BP = 12.19, df = 3, p-value = 0.00676

# o modelo de Fama-French rejeita a hipótese nula, portanto há indícios de heterocedasticidade

# D) Testes de Correlação Serial

# Teste de Durbin-Watson

dwtest(capm1)

## 
##  Durbin-Watson test
## 
## data:  capm1
## DW = 1.9032, p-value = 0.005819
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(ffrench1)

## 
##  Durbin-Watson test
## 
## data:  ffrench1
## DW = 1.9333, p-value = 0.04096
## alternative hypothesis: true autocorrelation is greater than 0

# O teste de Durbin-Watson em ambos os modelos indica que há correlação serial, por rejeitarem a hipótese nula.

# Teste de Breusch-Godfrey

bgtest(capm1)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  capm1
## LM test = 6.3642, df = 1, p-value = 0.01164

# O teste de Breusch-Godfrey indica que há correlação serial no modelo CAPM, rejeitando a hipótese nula. 

bgtest(ffrench1)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  ffrench1
## LM test = 2.9722, df = 1, p-value = 0.08471

# Já o teste de Breusch-Godfrey indica, por pouco, que não haveria correlação serial no modelo Fama-French, aceitando a hipótese nula.

QUESTÃO 2:

# Gerando CAPM e Fama-French sem outliers

outliers_capm <- which(apply(influence.measures(capm1)$is.inf, 1, any))
outliers_ffrench <- which(apply(influence.measures(ffrench1)$is.inf, 1, any))

capm2 <- lm(`Rp-Rf` ~ `Rm-Rf`, data = df_capm[-outliers_capm, ])
summary(capm2)

## 
## Call:
## lm(formula = `Rp-Rf` ~ `Rm-Rf`, data = df_capm[-outliers_capm, 
##     ])
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.222470 -0.009627 -0.000059  0.009771  0.076152 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0001153  0.0003430  -0.336    0.737    
## `Rm-Rf`      1.0063167  0.0261028  38.552   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01714 on 2496 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.3732, Adjusted R-squared:  0.373 
## F-statistic:  1486 on 1 and 2496 DF,  p-value: < 2.2e-16

ffrench2 <- lm(`Rp-Rf` ~ MKT + HML + SMB, data = fatores[-outliers_ffrench, ])
summary(ffrench2)

## 
## Call:
## lm(formula = `Rp-Rf` ~ MKT + HML + SMB, data = fatores[-outliers_ffrench, 
##     ])
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.231328 -0.009729  0.000040  0.009623  0.089260 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0003399  0.0003437  -0.989    0.323    
## MKT          0.8700760  0.0325098  26.763  < 2e-16 ***
## HML          0.8026690  0.0547316  14.666  < 2e-16 ***
## SMB         -0.3359015  0.0493906  -6.801  1.3e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01713 on 2480 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.3725, Adjusted R-squared:  0.3717 
## F-statistic: 490.6 on 3 and 2480 DF,  p-value: < 2.2e-16

# A) Teste de normalidade dos resíduos de Jarque-Bera

jarqueberaTest(residuals(capm2))

## 
## Title:
##  Jarque - Bera Normalality Test
## 
## Test Results:
##   STATISTIC:
##     X-squared: 14186.0263
##   P VALUE:
##     Asymptotic p Value: < 2.2e-16 
## 
## Description:
##  Mon Jul 12 22:59:21 2021 by user: Bruno Marcelino

jarqueberaTest(residuals(ffrench2))

## 
## Title:
##  Jarque - Bera Normalality Test
## 
## Test Results:
##   STATISTIC:
##     X-squared: 20108.7444
##   P VALUE:
##     Asymptotic p Value: < 2.2e-16 
## 
## Description:
##  Mon Jul 12 22:59:21 2021 by user: Bruno Marcelino

# Os dois testes indicam que os resíduos de ambos os modelos sem outliers seguem uma distribuição normal, consertando o problema anterior.

# B) Medidas de Influência das observações

# Gráficos de identificação dos outliers

par(mfrow=c(2,3))
plot(capm2, which=1:6)

par(mfrow=c(2,3))
plot(ffrench2, which=1:6)

# ainda que tenhamos retirado os outliers obtidos com a função influence.measures, os dois modelos apresentaram no mínimo três dados que apresentaram distância de Cook mais alta que os demais. O restante dos gráficos relata o mesmo. 

# C) Testes de heterocedasticidade de Breusch-Pagan

bptest(capm2)

## 
##  studentized Breusch-Pagan test
## 
## data:  capm2
## BP = 0.26599, df = 1, p-value = 0.606

bptest(ffrench2)

## 
##  studentized Breusch-Pagan test
## 
## data:  ffrench2
## BP = 3.9857, df = 3, p-value = 0.263

# Os dois testes indicam a ausência do problema da heterocedasticidade nos dois modelos, o que é um grande benefício para o ajuste da regressão de maneira próxima da realidade.

# D) Testes de Correlação Serial

# Teste de Durbin-Watson

dwtest(capm2)

## 
##  Durbin-Watson test
## 
## data:  capm2
## DW = 2.066, p-value = 0.9505
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(ffrench2)

## 
##  Durbin-Watson test
## 
## data:  ffrench2
## DW = 1.9965, p-value = 0.4645
## alternative hypothesis: true autocorrelation is greater than 0

# Entretanto, o modelo der Durbin-Watson indica que há altíssima correlação serial nos modelos, sobretudo no CAPM (95% de aceitação da hipótese nula). Ainda que seja um problema, era um problema a se esperar dado que se tratam de dados de séries temporais, e o retorno passado tem grande influência no retorno futuro. 

# Teste de Breusch-Godfrey

bgtest(capm2)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  capm2
## LM test = 2.7254, df = 1, p-value = 0.09877

bgtest(ffrench2) 

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  ffrench2
## LM test = 0.0046572, df = 1, p-value = 0.9456

# Já no teste de Breusch-Godfrey, os modelos apresentaram ausência de correlação serial após a retirada dos outliers, tendo em vista que não rejeitaram a hipótese nula. No caso do CAPM foi por pouco, mas para o Fama-French houveram mudanças significativas no ajuste do que fizeram o teste apresentar p-valor de 95%)

QUESTÃO 3:

### Correção da heterocedasticidade -> erros-padrão de White ### 
 
# sem correção
summary(capm2)

## 
## Call:
## lm(formula = `Rp-Rf` ~ `Rm-Rf`, data = df_capm[-outliers_capm, 
##     ])
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.222470 -0.009627 -0.000059  0.009771  0.076152 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0001153  0.0003430  -0.336    0.737    
## `Rm-Rf`      1.0063167  0.0261028  38.552   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01714 on 2496 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.3732, Adjusted R-squared:  0.373 
## F-statistic:  1486 on 1 and 2496 DF,  p-value: < 2.2e-16

summary(ffrench2)

## 
## Call:
## lm(formula = `Rp-Rf` ~ MKT + HML + SMB, data = fatores[-outliers_ffrench, 
##     ])
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.231328 -0.009729  0.000040  0.009623  0.089260 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0003399  0.0003437  -0.989    0.323    
## MKT          0.8700760  0.0325098  26.763  < 2e-16 ***
## HML          0.8026690  0.0547316  14.666  < 2e-16 ***
## SMB         -0.3359015  0.0493906  -6.801  1.3e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01713 on 2480 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.3725, Adjusted R-squared:  0.3717 
## F-statistic: 490.6 on 3 and 2480 DF,  p-value: < 2.2e-16

# com correção 
coeftest(capm2, vcov. = vcovHAC)

## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value Pr(>|t|)    
## (Intercept) -0.00011528  0.00032998 -0.3493   0.7269    
## `Rm-Rf`      1.00631675  0.03287387 30.6115   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(ffrench2, vcov. = vcovHAC)

## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) -0.00033991  0.00034071 -0.9976    0.3186    
## MKT          0.87007602  0.04229203 20.5730 < 2.2e-16 ***
## HML          0.80266903  0.06151166 13.0491 < 2.2e-16 ***
## SMB         -0.33590147  0.06132341 -5.4775 4.747e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Nesse teste de erros-padrão de White, é gerado um modelo que atribui menores pesos a informações discrepantes. Sendo assim, são corrigidos os valores dos coeficientes e os seus erros-padrão, tornando-o mais robusto (ou seja, menos sensível à presença de outliers).

### Correção da normalidade dos resíduos -> dummy que equivale a 1 se a observação for identificada como um outlier ###

df_capm <- df_capm %>% 
    mutate("position" = seq(1:length(date))) %>% 
    mutate("D1" = ifelse(position %in% outliers_capm, 1, 0))

fatores <- fatores %>% 
    mutate("position" = seq(1:length(date))) %>% 
    mutate("D2" = ifelse(position %in% outliers_capm, 1, 0))

# Análise da regressão com dummy que só vale para os outliers

capm3 <- lm(`Rp-Rf` ~ `Rm-Rf` + D1, data = df_capm)
summary(capm3) 

## 
## Call:
## lm(formula = `Rp-Rf` ~ `Rm-Rf` + D1, data = df_capm)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.221259 -0.010287 -0.000292  0.010373  0.108657 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0001164  0.0004144  -0.281    0.779    
## `Rm-Rf`      1.0586713  0.0249690  42.399  < 2e-16 ***
## D1           0.0061311  0.0014476   4.235 2.36e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02071 on 2718 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:    0.4,  Adjusted R-squared:  0.3996 
## F-statistic: 906.1 on 2 and 2718 DF,  p-value: < 2.2e-16

ffrench3 <- lm(`Rp-Rf` ~ MKT + HML + SMB + D2, data = fatores)
summary(ffrench3)

## 
## Call:
## lm(formula = `Rp-Rf` ~ MKT + HML + SMB + D2, data = fatores)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.230059 -0.010303 -0.000156  0.010232  0.127411 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0003120  0.0004145  -0.753    0.452    
## MKT          0.9932494  0.0293984  33.786  < 2e-16 ***
## HML          0.7704764  0.0550243  14.002  < 2e-16 ***
## SMB         -0.3923388  0.0507044  -7.738 1.42e-14 ***
## D2           0.0063151  0.0014497   4.356 1.37e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0207 on 2713 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.4018, Adjusted R-squared:  0.4009 
## F-statistic: 455.5 on 4 and 2713 DF,  p-value: < 2.2e-16

# Concluimos que os outliers não estão contidos nos padrões gerais da amostra, pois a variável dummy é significativa. Isso significa que sua presença não contribui para as especificações do modelo, estando muito fora dos padrões gerais da amostra (podendo ser incluída como outra variável e assim melhorando o ajuste do modelo por ser significativa)

### Correção da correlação serial

# Segundo os testes realizados em modelos que retiraram os outliers, pudemos concluir que a correlação serial já foi corrigida. Entretanto, a função "coeftest" gera um modelo corrigido para esse problema com o parâmetro (vcov. = vcovHAC)

QUESTÃO 3:

O modelo Fama-French se mostrou melhor que o CAPM no sentido de incluir a mesmas variável (prêmio de risco do mercado) em junção com outras que se demonstraram também significantes ao longo do processo, como o SMB e o HML. Além disso, vale ressaltar que vale a pena estipular sua regressão sem intercepto pois ele já se apresenta estatisticamente igual a zero. Seguem abaixo os problemas enfrentados por nós no módulo e suas respectivas conclusões:

Heterocedasticidade

A retirada dos outliers do modelo Fama-French tornou seus resultados de testes de heteroedasticidade altamente conclusivos para homoedasticidade, o que é ótimo e um pressuposto para a regressão de variáveis pelo método linear e de mínimos quadrados.

Pressuposto de normalidade dos resíduos

Os testes indicaram que os modelos sem outliers indicam confirmação do pressuposto de normalidade dos resíduos no novo modelo, o que é um forte benefício para seu ajuste à vida real.

Correlação Serial

Os testes não forneceram estatísticas realmente conclusivas de que o modelo sem outliers apresentou ou não correlação serial, o que ainda deve ser estudado mais a fundo por nós, os responsáveis pelos métodos econométricos. Como não foi possível armazenar os dados do modelo que corrige esse problema por meio do método de erros-padrão de White, a correção ou não da correlação (que deveria ter sido feita) ainda é uma incógnita de forma empírica.