Econometria das Séries Temporais - Atividade 4
Questão 1
O poder do teste de raiz unitária
Primeiro, criamos uma função capaz de gerar séries de acordo com as especificações que quisermos. Também adicionamos a possibilidade de definir o set.seed para propósitos de reproducibilidade dos resultados.
gerar_serie <- function(beta, var_erro, n_obs){
erro <- rnorm(n_obs, mean = 0, sd = sqrt(var_erro))
y_t <- vector(length = n_obs)
y_t[1] <- erro[1] + 1
for (i in 2:length(y_t)) {
y_t[i] <- 1 + beta * y_t[i-1] + erro[i]
}
return(y_t)
}Definição dos Parâmetros
O exercício define um processo autorregressivo AR(1) \(y_t = \alpha + \beta y_{t-1} + \mu_t\). Também define a parametrização: \(\alpha = 1; \beta = \{0.95, 0.9, 0.85 \}; \mu \sim N(0,1)\: \& \: \mu \sim N(0,25); n = \{50, 100\}\).
grid <- expand.grid(beta = c(0.95, 0.9, 0.85), var_erro = c(1, 25), n_obs = c(50, 100))Construção das séries
lista_1000_series <- list()
lista_suporte <- vector("list", length = 1000)
for (j in 1:nrow(grid)) {
params <- grid[j, ]
nome <- paste(as.character(params), collapse = "-")
for (i in 1:1000) {
lista_suporte[[i]] <- gerar_serie(as.numeric(params[1]), as.numeric(params[2]), as.numeric(params[3]))
}
lista_1000_series[[nome]] <- as.data.frame(do.call("cbind", lista_suporte))
}## Length Class Mode
## 0.95-1-50 1000 data.frame list
## 0.9-1-50 1000 data.frame list
## 0.85-1-50 1000 data.frame list
## 0.95-25-50 1000 data.frame list
## 0.9-25-50 1000 data.frame list
## 0.85-25-50 1000 data.frame list
## 0.95-1-100 1000 data.frame list
## 0.9-1-100 1000 data.frame list
## 0.85-1-100 1000 data.frame list
## 0.95-25-100 1000 data.frame list
## 0.9-25-100 1000 data.frame list
## 0.85-25-100 1000 data.frame listTestes ADF e ADF-GLS
pvalor_adf <- lapply(lista_1000_series, function(df) do.call("c", lapply(df, function(x) attr(ur.df(x, type = c("drift"), lags = 0), 'testreg')$coefficients[2,3])))
pvalor_adf_gls <- lapply(lista_1000_series, function(df) do.call("c", lapply(df, function(x) attr(ur.ers(x, model = "constant", lag.max = 0), 'testreg')$coefficients[1,3]))) Valores Críticos de cada Teste
serie_50_obs <- gerar_serie(0.95, 25, 50)
serie_100_obs <- gerar_serie(0.95, 25, 100)Poder dos Testes
adf_prop5 <- vector(length = length(pvalor_adf))
adf_prop1 <- vector(length = length(pvalor_adf))
for (i in seq_along(pvalor_adf)) {
vetor_cval <- pvalor_adf[[i]]
if(i <= 6){ # n = 50
adf_prop5[i] <- sum(vetor_cval < val_critc[1,2]) / length(vetor_cval)
adf_prop1[i] <- sum(vetor_cval < val_critc[1,1]) / length(vetor_cval)
} else { # n = 100
adf_prop5[i] <- sum(vetor_cval < val_critc[2,2]) / length(vetor_cval)
adf_prop1[i] <- sum(vetor_cval < val_critc[2,1]) / length(vetor_cval)
}
}adf_gls_prop5 <- vector(length = length(pvalor_adf_gls))
adf_gls_prop1 <- vector(length = length(pvalor_adf_gls))
for (i in seq_along(pvalor_adf_gls)) {
vetor_cval <- pvalor_adf_gls[[i]]
if(i <= 6){ # n = 50
adf_gls_prop5[i] <- sum(vetor_cval < val_critc[3,2]) / length(vetor_cval)
adf_gls_prop1[i] <- sum(vetor_cval < val_critc[3,1]) / length(vetor_cval)
} else { # n = 100
adf_gls_prop5[i] <- sum(vetor_cval < val_critc[4,2]) / length(vetor_cval)
adf_gls_prop1[i] <- sum(vetor_cval < val_critc[4,1]) / length(vetor_cval)
}
}Questão 2
Utilize uma série de dados diários para a realização desse exercício. O interesse é modelar o retorno diário da série, portanto, caso a série seja uma taxa de câmbio, deverão logaritmizá-la e depois tirar a primeira diferença (ou trabalhar com a taxa de variação percentual diária). Deverão proceder da mesma maneira se a série for um índice financeiro ou o preço de um ativo financeiro.
Importamos a série de preços da Apple.
getSymbols("AAPL", from = '2005-01-01',
to = "2021-12-31",warnings = FALSE,
auto.assign = TRUE)
apple <- Ad(AAPL)apple$AAPL <- append(NA, diff(log(apple$AAPL)))## Data AAPL
## Min. :2005-01-04 Min. :-0.197470
## 1st Qu.:2009-04-03 1st Qu.:-0.008260
## Median :2013-07-04 Median : 0.001121
## Mean :2013-07-03 Mean : 0.001219
## 3rd Qu.:2017-10-01 3rd Qu.: 0.011891
## Max. :2021-12-30 Max. : 0.130194## 'data.frame': 4278 obs. of 2 variables:
## $ Data: Date, format: "2005-01-04" "2005-01-05" ...
## $ AAPL: num 0.010218 0.008719 0.000775 0.070282 -0.004196 ...Enunciado A
Mostre o correlograma da série de retornos e o correlograma do quadrado da mesma série. O correlograma da série sugere algum padrão? E o correlograma do quadrado?
Correlograma da série
### Correlograma da série ao quadrado
Enunciado B
Modele a série de retorno através de GARCH(p,q) e também através de pelo menos uma variação do modelo que contemple a possibilidade de ajuste assimétrico da série. Mostre que a distribuição empírica aproxima-se da distribuição teórica utilizada para a construção de verossimilhança.
Para a série de retornos, iremos considerar as seguintes especificações ARMA:
- ARMA(15, 3)
- ARMA(15, 1)
- ARMA(10, 1)
- ARMA(8, 1)
mod <- list()
mod[[1]] <- arima(apple$AAPL, order = c(15, 0, 3))
mod[[2]] <- arima(apple$AAPL, order = c(15, 0, 1))
mod[[3]] <- arima(apple$AAPL, order = c(10, 0, 1))
mod[[4]] <- arima(apple$AAPL, order = c(8, 0, 1))Com base nos resultados acima, iremos selecionar a especificação ARMA(15,3).
res_arma15_3 <- residuals(mod[[1]])
Agora analisamos os correlogramas dos resíduos ao quadrado do modelo ARMA(10,10).
res_arma15_3_quad <- res_arma15_3^2
Para essa série (quadrado dos resíduos do modelo acima) iremos analisar as seguintes especificações:
- ARMA(10, 10)
- ARMA(10, 6)
- ARMA(6, 3)
- ARMA(3, 3)
mod <- list()
mod[[1]] <- arima(res_arma15_3_quad, order = c(10, 0, 10))
mod[[2]] <- arima(res_arma15_3_quad, order = c(10, 0, 6))
mod[[3]] <- arima(res_arma15_3_quad, order = c(6, 0, 3))
mod[[4]] <- arima(res_arma15_3_quad, order = c(3, 0, 3))Com base nos resultados acima, iremos selecionar a especificação ARMA(6,3) para construir os modelos GARCH e EGARCH (GARCH Exponencial). Para construção da verossimilhança vamos utilizar uma distribuição normal.
Desse modo, a variância condicional será modelada como um GARCH(6,3). A média será modelada por um ARMA(10,10).
Estimação do GARCH
garch_mod <- ugarchspec(mean.model = list(armaOrder = c(10,10)),
variance.model = list(model = "sGARCH", garchOrder = c(6, 3)),
distribution.model = 'norm')
garch <- ugarchfit(apple$AAPL, spec = garch_mod)## Level Fixed Include Estimate LB UB
## mu 0.002 0 1 1 -0.122 0.122
## ar1 0.350 0 1 1 -4.000 4.000
## ar2 -0.220 0 1 1 -4.000 4.000
## ar3 -0.418 0 1 1 -4.000 4.000
## ar4 -0.462 0 1 1 -4.000 4.000
## ar5 0.321 0 1 1 -4.000 4.000
## ar6 0.329 0 1 1 -4.000 4.000
## ar7 0.378 0 1 1 -4.000 4.000
## ar8 -0.230 0 1 1 -4.000 4.000
## ar9 0.521 0 1 1 -4.000 4.000
## ar10 0.340 0 1 1 -4.000 4.000
## ma1 -0.357 0 1 1 -4.000 4.000
## ma2 0.219 0 1 1 -4.000 4.000
## ma3 0.414 0 1 1 -4.000 4.000
## ma4 0.494 0 1 1 -4.000 4.000
## ma5 -0.336 0 1 1 -4.000 4.000
## ma6 -0.367 0 1 1 -4.000 4.000
## ma7 -0.328 0 1 1 -4.000 4.000
## ma8 0.192 0 1 1 -4.000 4.000
## ma9 -0.516 0 1 1 -4.000 4.000
## ma10 -0.364 0 1 1 -4.000 4.000
## arfima 0.000 0 0 0 NA NA
## archm 0.000 0 0 0 NA NA
## mxreg 0.000 0 0 0 NA NA
## omega 0.000 0 1 1 0.000 0.431
## alpha1 0.098 0 1 1 0.000 1.000
## alpha2 0.040 0 1 1 0.000 1.000
## alpha3 0.000 0 1 1 0.000 1.000
## alpha4 0.000 0 1 1 0.000 1.000
## alpha5 0.021 0 1 1 0.000 1.000
## alpha6 0.001 0 1 1 0.000 1.000
## beta1 0.630 0 1 1 0.000 1.000
## beta2 0.007 0 1 1 0.000 1.000
## beta3 0.163 0 1 1 0.000 1.000
## gamma 0.000 0 0 0 NA NA
## eta1 0.000 0 0 0 NA NA
## eta2 0.000 0 0 0 NA NA
## delta 0.000 0 0 0 NA NA
## lambda 0.000 0 0 0 NA NA
## vxreg 0.000 0 0 0 NA NA
## skew 0.000 0 0 0 NA NA
## shape 0.000 0 0 0 NA NA
## ghlambda 0.000 0 0 0 NA NA
## xi 0.000 0 0 0 NA NA
Agora partimos para o EGARCH(6,3).
egarch_mod <- ugarchspec(mean.model = list(armaOrder = c(10,10)),
variance.model = list(model = "eGARCH", garchOrder = c(2, 1)),
distribution.model = 'norm')
egarch <- ugarchfit(apple$AAPL, spec = egarch_mod)## Level Fixed Include Estimate LB UB
## mu 0.004 0 1 1 -0.122 0.122
## ar1 0.994 0 1 1 -4.000 4.000
## ar2 -0.593 0 1 1 -4.000 4.000
## ar3 0.300 0 1 1 -4.000 4.000
## ar4 0.028 0 1 1 -4.000 4.000
## ar5 0.575 0 1 1 -4.000 4.000
## ar6 -0.438 0 1 1 -4.000 4.000
## ar7 0.522 0 1 1 -4.000 4.000
## ar8 -0.300 0 1 1 -4.000 4.000
## ar9 0.464 0 1 1 -4.000 4.000
## ar10 -0.552 0 1 1 -4.000 4.000
## ma1 -0.996 0 1 1 -4.000 4.000
## ma2 0.607 0 1 1 -4.000 4.000
## ma3 -0.295 0 1 1 -4.000 4.000
## ma4 -0.002 0 1 1 -4.000 4.000
## ma5 -0.602 0 1 1 -4.000 4.000
## ma6 0.439 0 1 1 -4.000 4.000
## ma7 -0.495 0 1 1 -4.000 4.000
## ma8 0.224 0 1 1 -4.000 4.000
## ma9 -0.415 0 1 1 -4.000 4.000
## ma10 0.536 0 1 1 -4.000 4.000
## arfima 0.000 0 0 0 NA NA
## archm 0.000 0 0 0 NA NA
## mxreg 0.000 0 0 0 NA NA
## omega -0.291 0 1 1 -10.000 10.000
## alpha1 -0.140 0 1 1 -10.000 10.000
## alpha2 0.045 0 1 1 -10.000 10.000
## beta1 0.963 0 1 1 -1.000 1.000
## gamma1 0.186 0 1 1 -10.000 10.000
## gamma2 -0.028 0 1 1 -10.000 10.000
## eta1 0.000 0 0 0 NA NA
## eta2 0.000 0 0 0 NA NA
## delta 0.000 0 0 0 NA NA
## lambda 0.000 0 0 0 NA NA
## vxreg 0.000 0 0 0 NA NA
## skew 0.000 0 0 0 NA NA
## shape 0.000 0 0 0 NA NA
## ghlambda 0.000 0 0 0 NA NA
## xi 0.000 0 0 0 NA NA
Questão 3
Projeção Local
Importação dos Dados
Voltando para os dados utilizados na questão 3, apresentamos quais são essas séries e de onde extraimos os dados:
- PIB (IBG - SIDRA);
- Inflação (IBGE - SIDRA);
- Taxa de Desemprego (Investing);
- Indicador de Incerteza Econômica (Yahii);
- SELIC (Receita Federal).
PIB
pib <- sidrar::get_sidra(api = "/t/6613/n1/all/v/all/p/all/c11255/90707/d/v9319%202") %>%
dplyr::select(`Trimestre (Código)`, Valor) %>%
set_names("Data", "Valores")
pib$Data <- ymd(19960301) %m+% months(seq(0, 308, by = 3))pib$Valores <- append(NA, diff(pib$Valores)/pib$Valore[-length(pib$Valore)])IPCA
ipca <- sidrar::get_sidra(api = "/t/1737/n1/all/v/63/p/all/d/v63%202") %>%
dplyr::select(`Mês (Código)`, Valor) %>%
set_names("Data", "Valores")
ipca$Data <- as.Date(ym(ipca$Data))
ipca$Valores <- ipca$Valores / 100
ipca <- ipca[-1, ]A série de inflação possui um “problema” a mais que a série de PIB: a frequência é mensal. Para transformar essa série em trimestral, primeiro calculamos a taxa acumulada; em seguida, transformamos essa série em trimestral; por fim, calculamos a variação percentual dessa série. Ou seja, no fim teremos a variação percentual da inflação em cada trimestre.
ipca$Valores <- cumprod(ipca$Valores + 1)
ipca <- xts(ipca[,-1], ipca$Data)
ipca <- ipca[endpoints(ipca, on = "quarters"), ]
ipca <- PerformanceAnalytics::Return.calculate(ipca)Desemprego
desemprego <- read_excel("Desemprego.xlsx")
data <- strsplit(desemprego$Data, ".", fixed = TRUE)
data <- lapply(data, function(x) paste0(x[3], x[2], x[1]))
data <- do.call("c", data)
desemprego$Data <- as.Date(lubridate::ymd(data))
valor <- substr(desemprego$Atual,1,nchar(desemprego$Atual)-1)
desemprego$Atual <- as.numeric(gsub(",", ".", valor))/100Para o desemprego, apenas transformamos os dados de mensal para trimestral.
desemprego <- xts(desemprego[,-1], desemprego$Data)
desemprego <- desemprego[endpoints(desemprego, on = "quarters"), ]Incerteza
theurl <- "http://www.yahii.com.br/iieBR.html"
file<- read_html(theurl)
tables<- html_nodes(file, "table")
lista_tabelas <- list()
for (i in 1:3) {
aux <- html_table(tables[i+5], fill = TRUE)[[1]]
colnames(aux) <- aux[1, ]
aux <- aux[-1, ]
aux$Data_mes <- 1:12
aux <- aux %>% dplyr::select(Data_mes, everything()) %>% dplyr::select(!`M/A`)
lista_tabelas[[i]] <- aux
}
incerteza <- merge(merge(lista_tabelas[[1]], lista_tabelas[[2]], by = "Data_mes"), lista_tabelas[[3]], by = "Data_mes")
incerteza <- incerteza[,c(-(ncol(incerteza)-1), -ncol(incerteza))]
incerteza <- incerteza %>% gather(year, value, -Data_mes)
incerteza$Data <- ym(200001) %m+% months(0:(nrow(incerteza)-1))
incerteza <- incerteza %>% dplyr::select(Data, value)
incerteza$value <- as.numeric(gsub(",", ".", incerteza$value))Para o índice de incerteza calculada pela FGV, primeiro transformamos a série em trimestral e, em seguida, calculamos a variação percentual.
incerteza <- xts(incerteza[,-1], incerteza$Data)
incerteza <- incerteza[endpoints(incerteza, on = "quarters"), ]
incerteza <- PerformanceAnalytics::Return.calculate(incerteza)Selic
theurl <- "https://www.gov.br/receitafederal/pt-br/assuntos/orientacao-tributaria/pagamentos-e-parcelamentos/taxa-de-juros-selic"
file<-read_html(theurl)
tables<-html_nodes(file, "table")
lista_tabelas <- list()
for (i in 1:3) {
aux <- html_table(tables[5-i], fill = TRUE)[[1]]
colnames(aux) <- aux[1, ]
aux <- aux[-1, ]
aux$Data_mes <- 1:12
aux <- aux %>% dplyr::select(Data_mes, everything()) %>% dplyr::select(!`Mês/Ano`)
lista_tabelas[[i]] <- aux
}
selic <- merge(merge(lista_tabelas[[1]], lista_tabelas[[2]], by = "Data_mes"), lista_tabelas[[3]], by = "Data_mes")
selic <- selic[,c(-(ncol(selic)-1), -ncol(selic))]
selic <- selic %>% gather(year, value, -Data_mes)
selic$Data <- ym(199501) %m+% months(0:(nrow(selic)-1))
selic <- selic %>% dplyr::select(Data, value)
selic$value <- substr(selic$value,1,nchar(selic$value)-1)
selic$value <- as.numeric(gsub(",", ".", selic$value))/100Para a Selic adotamos a mesma abordagem utilizada para a inflação.
selic$value <- cumprod(selic$value + 1)
selic <- xts(selic[,-1], selic$Data)
selic <- selic[endpoints(selic, on = "quarters"), ]
selic <- PerformanceAnalytics::Return.calculate(selic)Dados Final
Como podemos perceber, temos séries na mesma frequência (trimestral), mas com datas de inicio e fim diferentes. Para corrigir esses problemas, primeiro fazemos um “merge” em todas essas bases de dados. Como podemos perceber pelo “resumo” dos nossos dados, temos valores nulos frutos da diferença supracitada.
dados_lista <- list(pib, ipca, desemprego, incerteza, selic)
dados <- do.call("merge.xts", dados_lista)Para corrigir esse problema, restringimos nossa base de dados para trabalharmos apenas com aqueles trimestres que são comuns a todas as séries.
dados <- dados["2001-12-01/2019-12-01"]
dados <- data.frame(dados)Função de Resposta a Impulso
Iremos analisar o impacto de um impulso da taxa de crescimento do PIB sobre a taxa do IPCA.
period_frente <- 8
lista_irf <- list()
for (i in 1:(period_frente+1)) {
mod <- lm(data = dados, IPCA ~ lag(PIB, i-1) +
lag(PIB, i) + lag(IPCA, i) + lag(Desemprego, i) + lag(Incerteza, i) + lag(Selic, i) +
lag(PIB, i+1) + lag(IPCA, i+1) + lag(Desemprego, i+1) + lag(Incerteza, i+1) + lag(Selic, i+1))
matriz_rob_mod <- coeftest(mod, vcov = vcovHAC)
lista_irf[[i]] <- c(matriz_rob_mod[2], matriz_rob_mod[14])
}df_irf <- as.data.frame(t(do.call("cbind", lista_irf))) %>% set_names(c("Coef", "Std. Error"))
df_irf <- df_irf %>% mutate(Superior = Coef + `Std. Error`, Inferior = Coef - `Std. Error`, Tempo = 0:8)