Usaremos o conjunto de pacotes "tidyverse" para realizar algumas operações, sobretudo utilizando recursos gráficos.
# ---- Importação das Bibliotecas ----
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("quantmod")
library("readxl")
library("GetQuandlData")
library("AER")
library("strucchange")
# ---- Importação dos Dados ----
# Cotações
ativos <- c("VALE3.SA", "^BVSP")
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")
# ---- 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 ----
# (adaptar os arquivos diretamente no Excel, manipulando a data para o formato YYYY-MM-DD)
# Posso usar a função as.Date(paste(year, month, day, sep = "-")) também
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
HML <- HML_Factor %>% select(date, HML)
HML$date <- as.Date(HML$date)
MKT <- MKT_Factor %>% select(date, Rm_minus_Rf) %>% rename("MKT" = Rm_minus_Rf)
MKT$date <- as.Date(MKT$date)
SMB <- SMB_Factor %>% select(date, SMB)
SMB$date <- as.Date(SMB$date)
Montando gráfico com os retornos da ação escolhida
# Gráfico CAPM
graf_1 <- df_capm %>%
ggplot(aes(x = `Rm-Rf`, y = `Rp-Rf`)) +
geom_point() +
labs(x = "Prêmio de Risco do Mercado",
y = "Prêmio de Risco do Ativo",
title = "CAPM",
subtitle = "Empresa: VALE S.A",
caption = "Dados coletados do Yahoo Finance e da API do Quandl") +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_smooth(method = "lm", color = "red") +
theme_bw()
graf_1
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 6 rows containing non-finite values (stat_smooth).
## Warning: Removed 6 rows containing missing values (geom_point).
graf_2 <- MKT %>%
ggplot(aes(x = date, y = MKT)) +
geom_point() +
labs(x = "",
y = "MKT",
title = "Fator MKT",
caption = "Dados coletados do site NEFIN-USP") +
geom_hline(yintercept = 0, linetype = "dashed") +
theme_bw()
graf_2
graf_3 <- HML %>%
ggplot(aes(x = date, y = HML)) +
geom_point() +
labs(x = "",
y = "HML",
title = "Fator HML",
caption = "Dados coletados do site NEFIN-USP") +
geom_hline(yintercept = 0, linetype = "dashed") +
theme_bw()
graf_3
graf_4<- SMB %>%
ggplot(aes(x = date, y = SMB)) +
geom_point() +
labs(x = "",
y = "SMB",
title = "Fator SMB",
caption = "Dados coletados do site NEFIN-USP") +
geom_hline(yintercept = 0, linetype = "dashed") +
theme_bw()
graf_4
SMB <- SMB_Factor %>% select(date, SMB)
SMB$date <- as.Date(SMB$date)
# Gráfico Retornos
graf_5 <- df_retornos %>%
ggplot(aes(x = date, y = VALE3.SA)) +
geom_line() +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_y_continuous(labels = scales::percent) +
labs(title = "Retornos Mensais",
subtitle = "Empresa: Petrobrás S.A",
y = "Retornos Mensais",
x = "") +
theme_tq() +
scale_fill_tq()
graf_5
# ---- Montando regressão de três fatores de Fama-French ----
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
O teste F verifica se o modelo restrito é valido ou não. Quando a estatística F for alta, assim como obtivemos no resultado do teste, terá um p-valor significativo. Isso demonstra que a média condicional (que é estimada pelo modelo) é superior a média incondicional, logo, rejeita-se a hipótese numa que os dois modelos tem a mesma eficácia.
# Teste F entre os dois modelos
anova(capm1,ffrench1)
## Warning in anova.lmlist(object, ...): models with response '"Rp-Rf"' removed
## because response differs from model 1
## Analysis of Variance Table
##
## Response: df_capm$`Rp-Rf`
## Df Sum Sq Mean Sq F value Pr(>F)
## df_capm$`Rm-Rf` 1 0.76968 0.76968 1783.2 < 2.2e-16 ***
## Residuals 2719 1.17359 0.00043
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Valores médios para VIF superiores a 5 s?o um indício de multicolinearidade, sendo assim não há multicolineariade nos testes pois todos os valores são menores que 5.
# VIF do modelo Fama-French
vif(ffrench1)
## MKT HML SMB
## 1.101118 1.221621 1.141517
Foi obtido estatísticas de teste altas, de forma que o p-valor foi significativo para os dois modelos. Isso indica que há uma forma funcional ou especificação alternativa que melhor valide o modelo. Fazendo o teste para os regressores também identificamos p-valores significativos, onde devemos verificar se as variáveis independentes são melhores ajustadas ao modelo se tratadas de forma polinomial.
# Teste RESET em cada modelo
resettest(capm1, power = c(2,3,4))
##
## RESET test
##
## data: capm1
## RESET = 6.5551, df1 = 3, df2 = 2716, p-value = 0.0002056
resettest(ffrench1, power = c(2,3,4))
##
## RESET test
##
## data: ffrench1
## RESET = 7.153, df1 = 3, df2 = 2711, p-value = 8.784e-05
# Avalia o teste como potências dos regressores, e não potências de todas as variáveis em conjunto
resettest(ffrench1, power = c(2,3,4), type = "regressor")
##
## RESET test
##
## data: ffrench1
## RESET = 6.7626, df1 = 9, df2 = 2705, p-value = 1.191e-09
# Teste de Chow (quebra estrutural)
sctest(df_capm$`Rp-Rf` ~ df_capm$`Rm-Rf`, data = df_capm, type = "Chow")
##
## Chow test
##
## data: df_capm$`Rp-Rf` ~ df_capm$`Rm-Rf`
## F = 3.2026, p-value = 0.04081
sctest(`Rp-Rf` ~ MKT + HML + SMB, data = fatores, type = "Chow")
##
## Chow test
##
## data: `Rp-Rf` ~ MKT + HML + SMB
## F = 6.64, p-value = 2.579e-05
# Gráfico do teste F para subgrupos da amostra
capm_f <- Fstats(df_capm$`Rp-Rf` ~ df_capm$`Rm-Rf`, data = df_capm, from = 0.1, to = 0.9)
plot(capm_f, main = "CAPM Estatísticas F") # Gráfico do F
plot(capm_f, pval = T, main = "CAPM p-valores") # Gráfico do p-valor
ffrench1_f <- Fstats(`Rp-Rf` ~ MKT + HML + SMB, data = fatores, from = 0.1, to = 0.9)
plot(ffrench1_f, main = "Fama-French Estatísticas F") # Gráfico do F
plot(ffrench1_f, pval = T, main = "Fama-French p-valores") # Gráfico do p-valor
O modelo de Fama-French é preferível, dado o teste F;
O teste RESET indica que há formas funcionais dos REGRESSORES que se ajustam melhor ao modelo, tanto no CAPM quanto no Fama-French. Isso indica que o fator “prêmio de risco do mercado” pode ser ajustado para a obtenção de melhores resultados;
O Teste de Chow indica que há quebras estruturais em subgrupos das amostras, ou seja, outra regressão pode ser aplicada em algum ponto do modelo para que este seja melhor ajustado.