Dados

Ajusto um modelo da classe VAR aos log-retornos diários do Ibovespa e da Cemig. Usando métodos automatizados para a identificação da ordem do modelo encontro a possibilidade de dois modelos: VAR (1) e VAR(6). Estimo os dois modelos e verifico seus resíduos usando o teste Portmanteau, o teste Breusch-Pagan e o teste ARCH. Os resultados não são animadores para ambos os modelos, pois parece haver autocorrelação tanto nos resíduos como no quadrado dos resíduos, indicando heterocedasticidade condicional. Ao final são apresentadas as previsões (fora da amostra) de ambos os modelos.

# Para ler os dados em excel
suppressWarnings(library(readxl))
suppressWarnings(library(vars))

## Loading required package: MASS

## Loading required package: strucchange

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: sandwich

## Loading required package: urca

## Loading required package: lmtest

suppressWarnings(library(knitr))
#options(knitr.table.format = "latex")
# Arquivos
ibov <- read.table("D-IBV.txt", quote="\"", comment.char="")
cemig <- read.table("D-CEMIG.txt", quote="\"", comment.char="")

# Transforma em objeto time.series
ibov <- ts(ibov)
cemig <- ts(cemig)
# Gráfico das séries
plot.ts(ibov, xlab = "", ylab = "", main = "Ibovespa")

plot.ts(cemig, xlab = "", ylab = "", main = "Cemig")

# Computa o log-retorno das séries
ribov <- log(ibov) - log(lag(ibov, -1))
rcemig <- log(cemig) - log(lag(cemig, -1))

Abaixo segue o gráfico das séries dos log-retornos.

Identificação

Para visualizar a relação entre as séries, ploto o correlograma cruzado delas. Em seguida apresento as matrizes de correlação das seis primeiras defasagens e uma representação pictórica, que resume a informação exibida nestas primeiras defasagens.

# Junta as duas séries
x <- ts.intersect(ribov, rcemig)

# ACF das séries (correlograma cruzado)
acf(x, xlab ="", ylab = "")

# Grava os valores estimados das correlações
a <- acf(x, plot = FALSE)
correl <- matrix(a$acf, ncol = 4)

ro <- list()
for (i in 2:7){# Note que o loop começa em 2, pois o valor 1 corresponde à matriz de correlação de ordem 0
  
  ro_atual <- list(matrix(correl[i,], ncol = 2))
  ro <- c(ro, ro_atual)
}
ro

## [[1]]
##            [,1]       [,2]
## [1,] 0.05415033 0.05428396
## [2,] 0.12771439 0.13827518
## 
## [[2]]
##              [,1]        [,2]
## [1,] -0.004689109 -0.03224021
## [2,]  0.005565351 -0.02738206
## 
## [[3]]
##             [,1]        [,2]
## [1,] -0.05590189 -0.03487328
## [2,] -0.05922101 -0.05677269
## 
## [[4]]
##             [,1]        [,2]
## [1,] -0.05725903 -0.07202516
## [2,] -0.05125774 -0.08500187
## 
## [[5]]
##             [,1]        [,2]
## [1,] -0.07452499 -0.07617628
## [2,] -0.07194796 -0.10180511
## 
## [[6]]
##             [,1]        [,2]
## [1,] -0.09341139 -0.06325495
## [2,] -0.10918062 -0.10815005

Abaixo segue uma representação pictórica das matrizes acima em que um sinal de positivo ou negativo indica uma autocorrelação significante com o resepctivo sinal, e um ponto indica uma autocorrelação não-significante. Usou-se como intervalo de confiança \(2/\sqrt{T}\). Como se vê há correlação entre as séries e das séries consigo mesmas.

\[\begin{align*} \text{Lag 1 } & \text{ Lag 2 } & \text{ Lag 3 } & \text{ Lag 4 } & \text{ Lag 5 } & \text{ Lag 6}\\ \begin{bmatrix} + & +\\ + & + \end{bmatrix} & \begin{bmatrix} \cdot & \cdot \\ \cdot & \cdot \end{bmatrix} & \begin{bmatrix} - & \cdot \\ - & - \end{bmatrix} & \begin{bmatrix} - & -\\ - & \cdot \end{bmatrix} & \begin{bmatrix} - & -\\ - & - \end{bmatrix} & \begin{bmatrix} - & -\\ - & - \end{bmatrix} \end{align*}\]

Usando a função VARselect seleciono as ordens \(1\) e \(6\) baseado em diferentes critérios de informação. Este resultado permanece o mesmo para diferentes especificações do modelo: com constante, com tendência, com constante e tendência, e sem

# Seleciona a ordem do modelo VAR
VARselect(x, lag.max = 8, type = "const")$selection

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      6      1      1      6

VARselect(x, lag.max = 8, type = "trend")$selection

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      6      1      1      6

VARselect(x, lag.max = 8, type = "both")$selection

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      6      1      1      6

VARselect(x, lag.max = 8, type = "none")$selection

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      6      1      1      6

Estimação

Denote \(\text{Ibov}_{t}\) e \(\text{Cemig}_{t}\) o log-retorno, em \(t\), das séries Ibovespa e Cemig, respectivamente; sejam, também, \(Z_{t} = (\text{Ibov}_{t} \, \text{Cemig}_{t})'\) um vetor \(2\text{x}1\), \(\Phi_{0} = (\Phi_{10} \, \Phi_{20})'\) um vetor \(2\text{x}1\), \(\Phi_{1}\) uma matriz \(2\text{x}2\), \(\delta = (\delta_{1} \, \delta_{2})'\) um vetor \(2\text{x}1\) e \(a_{t} = (a_{1t} \, \, a_{2t})'\) um vetor \(2\text{x}1\), onde \(a_{t} = \text{RB}(\vec{0}, \Sigma)\) é ruído branco multivariado. O modelo geral, VAR(1), com constante e tendência tem a forma abaixo.

\[ \begin{bmatrix} \text{Ibov}_{t} \\ \text{Cemig}_{t} \end{bmatrix} = \begin{bmatrix} \Phi_{10}\\ \Phi_{20} \end{bmatrix} + \begin{bmatrix} \Phi_{11} & \Phi_{12}\\ \Phi_{21} & \Phi_{22} \end{bmatrix} \begin{bmatrix} \text{Ibov}_{t-1} \\ \text{Cemig}_{t-1} \end{bmatrix} + \begin{bmatrix} \delta_{1}\\ \delta_{2} \end{bmatrix}t + \begin{bmatrix} a_{1t} \\ a_{2t} \end{bmatrix} \] Ou, de maneira mais condensada,

\[ Z_{t} = \Phi_{0} + \Phi_{1}Z_{t-1} + \delta t + a_{t} \]

Sigo o critério SC e estimo um modelo VAR (1) com constante e tendência. Como os termos de constante e tendência são não-significantes, estimo um modelo somente com constante e outro somente com tendência. Novamente, as estimativas pontuais desses parâmetros não são significantes; de fato, com exceção de \(\Phi_{22}\), nenhum dos parâmetros estimados em qualquer um desses modelos foi significante a níveis convencionais. Por questão de parcimônia prefiro o modelo VAR (1) estimado abaixo.

\[\begin{align*} \text{Ibov}_{t} & = \underset{(0.04360)}{0.03108}\text{ Ibov}_{t-1} + \underset{(0.03445)}{0.02346}\text{ Cemig}_{t-1}\\ \text{Cemig}_{t} & = \underset{(0.05467)}{0.05968}\text{ Ibov}_{t-1} + \underset{(0.04320)}{0.10049}\text{ Cemig}_{t-1} \end{align*}\]

onde o número entre parênteses é erro-padrão estimado do coeficiente.

summary(fit <- VAR(x, p = 1, type = "none"))

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: ribov, rcemig 
## Deterministic variables: none 
## Sample size: 1497 
## Log Likelihood: 6845.376 
## Roots of the characteristic polynomial:
## 0.1168 0.01476
## Call:
## VAR(y = x, p = 1, type = "none")
## 
## 
## Estimation results for equation ribov: 
## ====================================== 
## ribov = ribov.l1 + rcemig.l1 
## 
##           Estimate Std. Error t value Pr(>|t|)
## ribov.l1   0.03108    0.04360   0.713    0.476
## rcemig.l1  0.02346    0.03445   0.681    0.496
## 
## 
## Residual standard error: 0.0286 on 1495 degrees of freedom
## Multiple R-Squared: 0.00334, Adjusted R-squared: 0.002007 
## F-statistic: 2.505 on 2 and 1495 DF,  p-value: 0.08199 
## 
## 
## Estimation results for equation rcemig: 
## ======================================= 
## rcemig = ribov.l1 + rcemig.l1 
## 
##           Estimate Std. Error t value Pr(>|t|)  
## ribov.l1   0.05968    0.05467   1.092   0.2751  
## rcemig.l1  0.10049    0.04320   2.326   0.0201 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.03586 on 1495 degrees of freedom
## Multiple R-Squared: 0.02003, Adjusted R-squared: 0.01872 
## F-statistic: 15.28 on 2 and 1495 DF,  p-value: 2.701e-07 
## 
## 
## 
## Covariance matrix of residuals:
##            ribov   rcemig
## ribov  0.0008171 0.000827
## rcemig 0.0008270 0.001285
## 
## Correlation matrix of residuals:
##        ribov rcemig
## ribov  1.000  0.807
## rcemig 0.807  1.000

Diagnóstico dos resíduos

Primeiro faço o gráfico dos resíduos, das suas correlações e do qudrado das suas correlações. Nesta inspeção visual já é possível suspeitar de que há autocorrelação no quadrado dos resíduos do modelo, e também de alguma autocorrelação nos resíduos que não foi capturada pelo modelo. Este resultado inicial não é animador.

## Análise de resíduos
# Modelo 1: VAR(1)
# ACF dos resíduos
acf(resid(fit), xlab = "", ylab = "")

# ACF do quadrado dos resíduos
acf(resid(fit)^2, xlab = "", ylab = "")

O teste Portmanteau rejeita a hipótese nula de que os resíduos se comportem como ruído branco multivariado, a níveis convencionais de significância. De qualquer forma, reporto o teste para diversas ordens. Fica claro que os resíduos do modelo não estão bons, isto é, de que o modelo não está bem ajustado.

# Teste Portmanteau para diversas ordens
p <- c(); chi <- c(); graus <- c(); ordem <- c()
# Loop para calcular o teste em diversas ordens
for (i in 1:8){
  p[i] <- as.numeric(serial.test(fit, lags.pt = i+3, type = "PT.asymptotic")$serial$p.value)
  chi[i] <- as.numeric(serial.test(fit, lags.pt = i+3, type = "PT.asymptotic")$serial$statistic)
  graus[i] <- as.numeric(serial.test(fit, lags.pt = i+3, type = "PT.asymptotic")$serial$parameter)
  ordem[i] <- i+3
}
# Agrupa as informações num data.frame
PT <- data.frame(ordem = ordem, graus_liberdade = graus, estatística = chi, p_valor = p)
kable(PT, align = 'c')

O teste Breusch-Godfrey possibilita outra forma de verificar a presença de correlação entre os resíduos e as suas defasagens (autocorrelação serial). Novamente, rejeitamos a hipótese nula de que os resíduos não tenham correlação com suas defasagens. A interpretação é a mesma que a apresentada acima, para o teste Portmanteau, os resíduos do modelo não se comportam como ruído branco multivariado.

# Teste Breusch-Godfrey para diversas ordens
p <- c(); LM <- c(); graus <- c(); ordem <- c()
# Loop para calcular o teste para diversas ordens
for (i in 1:10){
  p[i] <- as.numeric(serial.test(fit, lags.bg = i, type = "BG")$serial$p.value)
  LM[i] <- as.numeric(serial.test(fit, lags.bg = i, type = "BG")$serial$statistic)
  graus[i] <- as.numeric(serial.test(fit, lags.bg = i, type = "BG")$serial$parameter)
  ordem[i] <- i
}

BG <- data.frame(ordem = ordem, graus_liberdade = graus, estatística = LM, p_valor = p)
kable(BG, align = 'c', col.names = c("Ordem", "Graus de liberdade", "Estatística", "P-valor"), caption = "Resultados do teste Breusch-Godfrey")

# Parece não haver correlação (linear) entre os resíduos das equações e as variáveis defasadas e os resíduos defasados

Por fim, reporto o resultado do teste ARCH, para verificar a presença de autocorrelação no quadrado dos resíduos do modelo. Como se havia suspeitado pela inspeção visual dos correlogramas, rejeitamos a hipótese nula de que não há autocorrelação no quadrado dos resíduos. Abaixo reporto o teste para diversas ordens. Note que faço a versão multivariada do teste.

# Teste ARCH (multivariado) para diversas ordens
p <- c(); chi <- c(); graus <- c(); ordem <- c()
# Loop para calcular o teste para diversas ordens
for (i in 1:12){
  p[i] <- as.numeric(arch.test(fit, lags.multi = i)$arch.mul$p.value)
  chi[i] <- as.numeric(arch.test(fit, lags.multi = i)$arch.mul$statistic)
  graus[i] <- as.numeric(arch.test(fit, lags.multi = i)$arch.mul$parameter)
  ordem[i] <- i
}

ARCH <- data.frame(ordem = ordem, graus_liberdade = graus, estatística = chi, p_valor = p)
kable(ARCH, align = 'c', col.names = c("Ordem", "Graus de liberdade", "Estatística", "P-valor"), caption = "Resultados do teste ARCH-LM")

Como mencionado acima, os resultados não são bons e o modelo deveria ser revisto. Uma especificação alternativa, sugerida pelo critério AIC, é de um VAR(6). Abaixo reporto os resultados deste modelo.

Estimação VAR(6)

Como os resíduos do VAR(1) sugerem que este modelo não está bem ajustado aos dados, estimo um VAR(6). Novamente, estimo primeiro o modelo com constante e tendência. Como encontro que ambos estes parâmetros são não significantes, tento um modelo com constante e outro com tendência. Por fim, chego num modelo sem constante e sem tendência. Os resultados da estimação são reportados abaixo. Note que apenas as defasagens da própria série são significantes nas equações abaixo. Isto é, \(\Phi^{(3)}_{11}\) e \(\Phi^{(6)}_{11}\) são significantes para \(\hat{\text{Ibov}}_{t}\); e \(\Phi^{(1)}_{22}\), \(\Phi^{(2)}_{22}\), \(\Phi^{(4)}_{22}\) e \(\Phi^{(5)}_{22}\) são significantes para \(\hat{\text{Cemig}}_{t}\). Isto pode levantar suspeita de que não há granger-causalidade entre as séries.

As estimativas são resumidas nas matrizes abaixo, da esquerda para direita, lê-se \(\Phi_{1},\Phi_{2}, \Phi_{3}\) na primeira linha, seguida dos erros-padrão; na terceira linha, \(\Phi_{4},\Phi_{5}, \Phi_{6}\) com os respectivos erros-padrão em baixo:

\[\begin{align*} \begin{bmatrix} -0.001383 & 0.038963\\ 0.038564 & 0.103014 \end{bmatrix} \begin{bmatrix} 0.034991 & -0.055340 \\ 0.070705 & -0.095163 \end{bmatrix} \begin{bmatrix} -0.102461 & 0.048365 \\ -0.075297 & 0.006767 \end{bmatrix} \\ \begin{pmatrix} 0.043635 & 0.034813\\ 0.054705 & 0.043645 \end{pmatrix} \begin{pmatrix} 0.043410 & 0.034710 \\ 0.054423 & 0.043516 \end{pmatrix} \begin{pmatrix} 0.043361 & 0.034712 \\ 0.054362 & 0.043519 \end{pmatrix} \\ \end{align*}\]

\[\begin{align*} \begin{bmatrix} -0.011346 & -0.043614\\ 0.051028 & -0.099660 \end{bmatrix} & \begin{bmatrix} -0.044283 & -0.021964\\ 0.020575 & -0.088661 \end{bmatrix} & \begin{bmatrix} -0.115764 & 0.028714\\ -0.069425 & -0.046751 \end{bmatrix} \\ \begin{pmatrix} 0.043386 & 0.034719\\ 0.054393 & 0.043527 \end{pmatrix} & \begin{pmatrix} 0.043385 & 0.034725\\ 0.054391 & 0.043535 \end{pmatrix} & \begin{pmatrix} 0.043272 & 0.034399\\ 0.054250 & 0.043126 \end{pmatrix} \end{align*}\]

summary(fit2 <- VAR(x, p = 6, type = "none"))

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: ribov, rcemig 
## Deterministic variables: none 
## Sample size: 1492 
## Log Likelihood: 6869.173 
## Roots of the characteristic polynomial:
## 0.758 0.758 0.704 0.704 0.689 0.689 0.6292 0.6292 0.6135 0.6135 0.6063 0.6063
## Call:
## VAR(y = x, p = 6, type = "none")
## 
## 
## Estimation results for equation ribov: 
## ====================================== 
## ribov = ribov.l1 + rcemig.l1 + ribov.l2 + rcemig.l2 + ribov.l3 + rcemig.l3 + ribov.l4 + rcemig.l4 + ribov.l5 + rcemig.l5 + ribov.l6 + rcemig.l6 
## 
##            Estimate Std. Error t value Pr(>|t|)   
## ribov.l1  -0.001383   0.043635  -0.032  0.97473   
## rcemig.l1  0.038963   0.034813   1.119  0.26324   
## ribov.l2   0.034991   0.043410   0.806  0.42034   
## rcemig.l2 -0.055340   0.034710  -1.594  0.11107   
## ribov.l3  -0.102461   0.043361  -2.363  0.01826 * 
## rcemig.l3  0.048365   0.034712   1.393  0.16374   
## ribov.l4  -0.011346   0.043386  -0.262  0.79373   
## rcemig.l4 -0.043614   0.034719  -1.256  0.20925   
## ribov.l5  -0.044283   0.043385  -1.021  0.30756   
## rcemig.l5 -0.021964   0.034725  -0.633  0.52715   
## ribov.l6  -0.115764   0.043272  -2.675  0.00755 **
## rcemig.l6  0.028714   0.034399   0.835  0.40401   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.02818 on 1480 degrees of freedom
## Multiple R-Squared: 0.0276,  Adjusted R-squared: 0.01971 
## F-statistic:   3.5 on 12 and 1480 DF,  p-value: 3.955e-05 
## 
## 
## Estimation results for equation rcemig: 
## ======================================= 
## rcemig = ribov.l1 + rcemig.l1 + ribov.l2 + rcemig.l2 + ribov.l3 + rcemig.l3 + ribov.l4 + rcemig.l4 + ribov.l5 + rcemig.l5 + ribov.l6 + rcemig.l6 
## 
##            Estimate Std. Error t value Pr(>|t|)  
## ribov.l1   0.038564   0.054705   0.705   0.4810  
## rcemig.l1  0.103014   0.043645   2.360   0.0184 *
## ribov.l2   0.070705   0.054423   1.299   0.1941  
## rcemig.l2 -0.095163   0.043516  -2.187   0.0289 *
## ribov.l3  -0.075297   0.054362  -1.385   0.1662  
## rcemig.l3  0.006767   0.043519   0.156   0.8764  
## ribov.l4   0.051028   0.054393   0.938   0.3483  
## rcemig.l4 -0.099660   0.043527  -2.290   0.0222 *
## ribov.l5   0.020575   0.054391   0.378   0.7053  
## rcemig.l5 -0.088661   0.043535  -2.037   0.0419 *
## ribov.l6  -0.069425   0.054250  -1.280   0.2008  
## rcemig.l6 -0.046751   0.043126  -1.084   0.2785  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.03533 on 1480 degrees of freedom
## Multiple R-Squared: 0.05006, Adjusted R-squared: 0.04236 
## F-statistic:   6.5 on 12 and 1480 DF,  p-value: 2.115e-11 
## 
## 
## 
## Covariance matrix of residuals:
##            ribov    rcemig
## ribov  0.0007927 0.0008005
## rcemig 0.0008005 0.0012478
## 
## Correlation matrix of residuals:
##         ribov rcemig
## ribov  1.0000 0.8049
## rcemig 0.8049 1.0000

Diagnósitico de resíduos, VAR(6)

A análise visual do correlograma dos resíduos indica que o VAR(6) pode ser um melhor ajuste do que o VAR(1); contudo, parece ainda haver correlação em algumas ordens mais elevadas. Novamente, há evidência de que há autocorrelação no quadrado dos resíduos.

# Segundo modelo: VAR(6)
# ACF dos resíduos
acf(resid(fit2), xlab = "", ylab = "")

# ACF do quadrado dos resíduos
acf(resid(fit2)^2, xlab = "", ylab = "")

O teste Portmanteau mostra que parece haver autocorrelação serial nos resíduos para ordens maiores que 8.

# Teste Portmanteau de autocorrelação serial nos resíduos
serial.test(fit2, lags.pt=8, type="PT.asymptotic")

## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object fit2
## Chi-squared = 11.893, df = 8, p-value = 0.1561

# Teste Portmanteau para diversas ordens
p <- c(); chi <- c(); graus <- c(); ordem <- c()
# Loop para calcular o teste em diversas ordens
for (i in 1:8){
  p[i] <- as.numeric(serial.test(fit2, lags.pt = 7+i, type = "PT.asymptotic")$serial$p.value)
  chi[i] <- as.numeric(serial.test(fit2, lags.pt = 7+i, type = "PT.asymptotic")$serial$statistic)
  graus[i] <- as.numeric(serial.test(fit2, lags.pt = 7+i, type = "PT.asymptotic")$serial$parameter)
  ordem[i] <- 7+i
}
# Agrupa as informações num data.frame
PT <- data.frame(ordem = ordem, graus_liberdade = graus, estatística = chi, p_valor = p)
kable(PT, align = 'c', col.names = c("Ordem", "Graus de liberdade", "Estatística", "P-valor"), caption = "Resultados do teste Portmanteau")

O teste Breusch-Godfrey vai no mesmo sentido e indica que há autocorrelação nos resíduos do modelo. Abaixo reporto o teste para diversas ordens.

# Teste Breusch-Godfrey para diversas ordens
p <- c(); LM <- c(); graus <- c(); ordem <- c()
# Loop para calcular o teste para diversas ordens
for (i in 1:10){
  p[i] <- as.numeric(serial.test(fit2, lags.bg = i, type = "BG")$serial$p.value)
  LM[i] <- as.numeric(serial.test(fit2, lags.bg = i, type = "BG")$serial$statistic)
  graus[i] <- as.numeric(serial.test(fit2, lags.bg = i, type = "BG")$serial$parameter)
  ordem[i] <- i
}
BG <- data.frame(ordem = ordem, graus_liberdade = graus, estatística = LM, p_valor = p)
kable(BG, align = 'c', col.names = c("Ordem", "Graus de liberdade", "Estatística", "P-valor"), caption = "Resultados do teste Breusch-Godfrey")

Por fim, o teste ARCH indica, como já estava patente no correlograma, que há autocorrelação no quadrado dos resíduos.

# Teste ARCH (multivariado) para diversas ordens
p <- c(); chi <- c(); graus <- c(); ordem <- c()
# Loop para calcular o teste para diversas ordens
for (i in 1:12){
  p[i] <- as.numeric(arch.test(fit2, lags.multi = i)$arch.mul$p.value)
  chi[i] <- as.numeric(arch.test(fit2, lags.multi = i)$arch.mul$statistic)
  graus[i] <- as.numeric(arch.test(fit2, lags.multi = i)$arch.mul$parameter)
  ordem[i] <- i
}

ARCH <- data.frame(ordem = ordem, graus_liberdade = graus, estatística = chi, p_valor = p)
kable(ARCH, align = 'c', col.names = c("Ordem", "Graus de liberdade", "Estatística", "P-valor"), caption = "Resultados do teste ARCH-LM")

Previsão

Como argumentado acima, os modelos propostos não se ajustaram bem aos dados. Em ambos os casos, parece haver autocorrelações de ordens elevadas que não foram capturadas pelos modelos. Além disso, há claramente problemas de heterocedasticidade condicional na série. Por uma questão de completude apresento as previsões, 20 passos a frente, para o modelo VAR(6) .

plot(predict(fit2, n.ahead = 20), xlim = c(1450, 1520), ylim = c(-.1, .1))

Questão 2

# Abre os arquivos
ibov <- read_excel("ibv2.xls")
ind <- read_excel("indices.xls")

# Transforma em objeto time.series
ibov <- ts(ind$Ibv)
merv <- ts(ind$Merv)
# Gráfico das séries
plot.ts(ibov)

plot.ts(merv)

# Computa o log-retorno das séries
ribov <- log(ibov) - log(lag(ibov, -1))
rmerv <- log(merv) - log(lag(merv, -1))
# Gráfico dos log-retornos
plot.ts(ribov, main = "Log-retorno do índice Ibovespa", xlab = "", ylab = "")

plot.ts(rmerv, main = "Log-retorno do índice Merv", xlab = "", ylab = "")

# Junta as duas séries
x <- ts.intersect(ribov, rmerv)

# ACF das séries (correlograma cruzado)
acf(x, xlab ="", ylab = "")

# Visualmente, não parece haver muita correlação entre as séries

# Grava os valores estimados das correlações numa matriz
a <- acf(x, plot = FALSE)
correl <- matrix(a$acf, ncol = 4)
# Lista os valores estimados para as 6 primeiras defasagens, numa lista
ro <- list()
for (i in 2:7){# o loop começa em 2, pois o valor 1 da lista corresponde às correlações de ordem 0
  
  ro_atual <- list(matrix(correl[i,], ncol = 2))
  ro <- c(ro, ro_atual)
}
ro

## [[1]]
##            [,1]       [,2]
## [1,] 0.03119903 0.02309126
## [2,] 0.03903730 0.03667003
## 
## [[2]]
##             [,1]        [,2]
## [1,] -0.02557302 -0.05245255
## [2,] -0.01429808  0.02909192
## 
## [[3]]
##             [,1]       [,2]
## [1,] -0.01591490 0.02080691
## [2,]  0.01871022 0.01156290
## 
## [[4]]
##              [,1]         [,2]
## [1,] -0.028748048 -0.011622372
## [2,] -0.005776735  0.005257297
## 
## [[5]]
##              [,1]        [,2]
## [1,] -0.052424961 -0.04258838
## [2,] -0.003654279  0.02527557
## 
## [[6]]
##            [,1]       [,2]
## [1,] 0.01126966 0.01403940
## [2,] 0.03600329 0.03135183

# Seleciona a ordem do modelo VAR
VARselect(x, lag.max = 8, type = "const")$selection

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      1      1      2

VARselect(x, lag.max = 8, type = "trend")$selection

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      1      1      2

VARselect(x, lag.max = 8, type = "both")$selection

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      1      1      2

VARselect(x, lag.max = 8, type = "none")$selection

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      1      1      2

# O crtitério AIC escolhe duas defasagem para as quatro especificações acima

summary(fit3 <- VAR(x, p = 2, type = "none"))

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: ribov, rmerv 
## Deterministic variables: none 
## Sample size: 1978 
## Log Likelihood: 9106.327 
## Roots of the characteristic polynomial:
## 0.2697 0.2651 0.1659 0.1659
## Call:
## VAR(y = x, p = 2, type = "none")
## 
## 
## Estimation results for equation ribov: 
## ====================================== 
## ribov = ribov.l1 + rmerv.l1 + ribov.l2 + rmerv.l2 
## 
##          Estimate Std. Error t value Pr(>|t|)  
## ribov.l1  0.02766    0.02546   1.086   0.2774  
## rmerv.l1  0.01258    0.02511   0.501   0.6164  
## ribov.l2 -0.00120    0.02546  -0.047   0.9624  
## rmerv.l2 -0.05151    0.02510  -2.053   0.0403 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.02563 on 1974 degrees of freedom
## Multiple R-Squared: 0.003928,    Adjusted R-squared: 0.00191 
## F-statistic: 1.946 on 4 and 1974 DF,  p-value: 0.1003 
## 
## 
## Estimation results for equation rmerv: 
## ====================================== 
## rmerv = ribov.l1 + rmerv.l1 + ribov.l2 + rmerv.l2 
## 
##          Estimate Std. Error t value Pr(>|t|)  
## ribov.l1  0.02924    0.02582   1.133   0.2575  
## rmerv.l1  0.02345    0.02547   0.921   0.3573  
## ribov.l2 -0.03716    0.02582  -1.439   0.1502  
## rmerv.l2  0.04522    0.02545   1.777   0.0758 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.026 on 1974 degrees of freedom
## Multiple R-Squared: 0.003852,    Adjusted R-squared: 0.001833 
## F-statistic: 1.908 on 4 and 1974 DF,  p-value: 0.1065 
## 
## 
## 
## Covariance matrix of residuals:
##           ribov     rmerv
## ribov 0.0006563 0.0003141
## rmerv 0.0003141 0.0006755
## 
## Correlation matrix of residuals:
##        ribov  rmerv
## ribov 1.0000 0.4717
## rmerv 0.4717 1.0000

## Análise de resíduos
# Gráfico dos resíduos
plot.ts(resid(fit3), main = "Gráfico dos resíduos do VAR(2)", xlab = "")

# ACF dos resíduos
acf(resid(fit3), main = "ACF dos resíduos do VAR(2)", xlab = "", ylab = "")

# ACF do quadrado dos resíduos
acf(resid(fit3)^2, main = "ACF do quadrado dos resíduos do VAR(2)", xlab = "", ylab = "")

# Teste Portmanteau de autocorrelação serial nos resíduos
serial.test(fit3, lags.pt=10, type="PT.asymptotic")

## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object fit3
## Chi-squared = 41.323, df = 32, p-value = 0.1251

# Teste Portmanteau para diversas ordens
p <- c(); chi <- c(); graus <- c(); ordem <- c()
# Loop para calcular o teste em diversas ordens
for (i in 1:8){
  p[i] <- as.numeric(serial.test(fit3, lags.pt = i+3, type = "PT.asymptotic")$serial$p.value)
  chi[i] <- as.numeric(serial.test(fit3, lags.pt = i+3, type = "PT.asymptotic")$serial$statistic)
  graus[i] <- as.numeric(serial.test(fit3, lags.pt = i+3, type = "PT.asymptotic")$serial$parameter)
  ordem[i] <- i+3
}
# Agrupa as informações num data.frame
PT <- data.frame(ordem = ordem, graus_liberdade = graus, estatística = chi, p_valor = p)
kable(PT, align = 'c', col.names = c("Ordem", "Graus de liberdade", "Estatística", "P-valor"), caption = "Resultados do teste Portmanteau")

# Teste Breusch-Godfrey para diversas ordens
p <- c(); LM <- c(); graus <- c(); ordem <- c()
# Loop para calcular o teste para diversas ordens
for (i in 1:10){
  p[i] <- as.numeric(serial.test(fit3, lags.bg = i, type = "BG")$serial$p.value)
  LM[i] <- as.numeric(serial.test(fit3, lags.bg = i, type = "BG")$serial$statistic)
  graus[i] <- as.numeric(serial.test(fit3, lags.bg = i, type = "BG")$serial$parameter)
  ordem[i] <- i
}

BG <- data.frame(ordem = ordem, graus_liberdade = graus, estatística = LM, p_valor = p)
kable(BG, align = 'c', col.names = c("Ordem", "Graus de liberdade", "Estatística", "P-valor"), caption = "Resultados do teste Breusch-Godfrey")

# Teste ARCH (multivariado) para diversas ordens
p <- c(); chi <- c(); graus <- c(); ordem <- c()
# Loop para calcular o teste para diversas ordens
for (i in 1:12){
  p[i] <- as.numeric(arch.test(fit3, lags.multi = i)$arch.mul$p.value)
  chi[i] <- as.numeric(arch.test(fit3, lags.multi = i)$arch.mul$statistic)
  graus[i] <- as.numeric(arch.test(fit3, lags.multi = i)$arch.mul$parameter)
  ordem[i] <- i
}

ARCH <- data.frame(ordem = ordem, graus_liberdade = graus, estatística = chi, p_valor = p)
kable(ARCH, align = 'c', col.names = c("Ordem", "Graus de liberdade", "Estatística", "P-valor"), caption = "Resultados do teste ARCH-LM")

plot(predict(fit3, n.ahead = 20), xlim = c(1960, 2001), ylim = c(-.08, .08))

Questão 3

Uso a função causality para testar a hipótese de granger-causalidade. No primeiro modelo estimado para a série da Cemig e do Ibovespa, encontrei que apenas \(\Phi_{22}\) era significante, isto é, que apenas a primeira defasagem de \(Cemig_{t}\) era significante para a equação de \(Cemig_{t}\). É possível suspeitar que as séries sejam não-acopladas. De fato, não encontramos evidência para rejeitar, a níveis convencionais, tanto a hipótese de que \(Cemig_{t-1} \nrightarrow Ibov_{t}\) como de que \(Cemig_{t-1} \nrightarrow Ibov_{t}\). Note, contudo, que rejeitamos a hipótese nula de que não há causalidade contemporânea à la Granger, isto é, \(Cemig_{t} \leftrightarrow Ibov_{t}\).

# Como nas equações estimadas verificamos que apernas defasagens de ribov eram significantes para ribov
# e que apenas defasagens de rcemig eram significantes para rcemig, podemos suspeitar, a princípio que 
# as séries não sejam acopladas

causality(fit, cause = "ribov")$`Granger`

## 
##  Granger causality H0: ribov do not Granger-cause rcemig
## 
## data:  VAR object fit
## F-Test = 1.1918, df1 = 1, df2 = 2990, p-value = 0.2751

causality(fit, cause = "rcemig")$`Granger`

## 
##  Granger causality H0: rcemig do not Granger-cause ribov
## 
## data:  VAR object fit
## F-Test = 0.46349, df1 = 1, df2 = 2990, p-value = 0.496

causality(fit, cause = "rcemig")$Instant

## 
##  H0: No instantaneous causality between: rcemig and ribov
## 
## data:  VAR object fit
## Chi-squared = 590.36, df = 1, p-value < 2.2e-16

# A níveis convencionais verificamos que não há evidência de que Ibov granger-causa Cemig e tampouco que 
# Cemig granger-causa Ibov
# Rejeitamos a níveis convencionais que não há causalidade contemporânea à la Granger entre as séries

Na estimação do segundo modelo para as séries da Cemig e do Ibovespa encontrei, analogamente ao caso apresetado acima, que apenas defasagens da própria série eram significantes para a série; isto é, apenas defasagens de \(Cemig_{t}\) eram significantes para \(Cemig_{t}\), e apenas defasagens de \(Ibov_{t}\) eram significantes para \(Ibov_{t}\). Novamente não encontramos evidência para concluir que há granger causalidade, em qualquer direção, entre as séries de retorno. Note que

# Como nas equações estimadas verificamos que apernas defasagens de ribov eram significantes para ribov
# e que apenas defasagens de rcemig eram significantes para rcemig, podemos suspeitar, a princípio que 
# as séries não sejam acopladas

causality(fit2, cause = "ribov")$`Granger`

## 
##  Granger causality H0: ribov do not Granger-cause rcemig
## 
## data:  VAR object fit2
## F-Test = 1.1299, df1 = 6, df2 = 2960, p-value = 0.342

causality(fit2, cause = "rcemig")$`Granger`

## 
##  Granger causality H0: rcemig do not Granger-cause ribov
## 
## data:  VAR object fit2
## F-Test = 1.2778, df1 = 6, df2 = 2960, p-value = 0.2639

causality(fit2, cause = "rcemig")$Instant

## 
##  H0: No instantaneous causality between: rcemig and ribov
## 
## data:  VAR object fit2
## Chi-squared = 586.56, df = 1, p-value < 2.2e-16

# A níveis convencionais verificamos que não há evidência de que Ibov granger-causa Cemig e tampouco que 
# Cemig granger-causa Ibov
# Rejeitamos a níveis convencionais que não há causalidade contemporânea à la Granger entre as séries

Na questão 2 encontramos um modelo em que a segunda defasagem da série de retornos da Merval é significante na equação dos retornos do Ibovespa. Isto pode ser um indício de que há uma relação de granger causalidade no sentido \(Merv_{t-2} \rightarrow Ibov_{t}\). Não encontramos, contudo, evidência para concluir que haja granger causalidade, a níveis convencionais de significância, entre as séries. Note que rejeitamos a hipótese de que não há causalidade contemporânea, à la Granger, entre as séries, isto é, \(Merv_{t} \leftrightarrow Ibov_{t}\).

# Como rmerv.l2 é significante para ibov pode haver uma relação de granger-causalidade

causality(fit3, cause = "rmerv")$`Granger`

## 
##  Granger causality H0: rmerv do not Granger-cause ribov
## 
## data:  VAR object fit3
## F-Test = 2.2136, df1 = 2, df2 = 3948, p-value = 0.1094

causality(fit3, cause = "ribov")$`Granger`

## 
##  Granger causality H0: ribov do not Granger-cause rmerv
## 
## data:  VAR object fit3
## F-Test = 1.6552, df1 = 2, df2 = 3948, p-value = 0.1912

causality(fit3, cause = "ribov")$Instant

## 
##  H0: No instantaneous causality between: ribov and rmerv
## 
## data:  VAR object fit3
## Chi-squared = 360.41, df = 1, p-value < 2.2e-16

# Marginalmente rejeitamos a 11% a hipótese de que rmerv não granger-causa ribov
# Rejeitamos a níveis convencionais que não há causalidade contemporânea à la Granger

Questão 4

# Pacotes adicionais
suppressWarnings(library(rugarch))

## Loading required package: parallel

## 
## Attaching package: 'rugarch'

## The following object is masked from 'package:stats':
## 
##     sigma

suppressWarnings(library(fExtremes))

## Loading required package: timeDate

## Loading required package: timeSeries

## 
## Attaching package: 'timeSeries'

## The following object is masked from 'package:zoo':
## 
##     time<-

## Loading required package: fBasics

## Loading required package: fGarch

suppressWarnings(library(astsa))

## 
## Attaching package: 'astsa'

## The following object is masked from 'package:fBasics':
## 
##     nyse

suppressWarnings(library(aTSA))

## 
## Attaching package: 'aTSA'

## The following object is masked from 'package:vars':
## 
##     arch.test

## The following object is masked from 'package:graphics':
## 
##     identify

# Importar os dados
IBM <- read.table("d-ibmln.txt", quote="\"", comment.char="")

# Transforma em ts()
ibm <- ts(IBM)
# A série já está em log-retorno
head(ibm)

## [1]  0.428 -0.428 -1.439  1.510  2.398  4.263

##         ACF  PACF
##  [1,]  0.00  0.00
##  [2,] -0.03 -0.03
##  [3,]  0.01  0.01
##  [4,] -0.01 -0.01
##  [5,]  0.03  0.03
##  [6,]  0.00 -0.01
##  [7,]  0.01  0.01
##  [8,]  0.00  0.00
##  [9,]  0.00  0.00
## [10,] -0.01 -0.01
## [11,]  0.00  0.00
## [12,]  0.02  0.02
## [13,]  0.00  0.00
## [14,]  0.01  0.02
## [15,]  0.01  0.01
## [16,]  0.01  0.01
## [17,]  0.01  0.01
## [18,] -0.01 -0.01
## [19,] -0.02 -0.02
## [20,]  0.03  0.03
## [21,]  0.01  0.01
## [22,] -0.02 -0.01
## [23,] -0.02 -0.02
## [24,]  0.01  0.01
## [25,] -0.01 -0.01
## [26,] -0.02 -0.02
## [27,]  0.01  0.01
## [28,]  0.02  0.02
## [29,]  0.01  0.01
## [30,] -0.01  0.00
## [31,] -0.02 -0.02
## [32,]  0.02  0.02
## [33,]  0.01  0.01
## [34,]  0.00  0.00
## [35,] -0.02 -0.02
## [36,] -0.02 -0.02

mean.model <- arima(ibm, order = c(1,0,2))
mean.model

## 
## Call:
## arima(x = ibm, order = c(1, 0, 2))
## 
## Coefficients:
##           ar1     ma1      ma2  intercept
##       -0.6723  0.6745  -0.0268     0.0389
## s.e.   0.1304  0.1307   0.0118     0.0154
## 
## sigma^2 estimated as 2.189:  log likelihood = -16184.4,  aic = 32378.81

ts.diag(mean.model)

mean.model2 <- arima(ibm, order = c(20,0,2), fixed = c(NA, rep(0, 18), NA, NA, NA, NA))

## Warning in arima(ibm, order = c(20, 0, 2), fixed = c(NA, rep(0, 18), NA, :
## some AR parameters were fixed: setting transform.pars = FALSE

mean.model2

## 
## Call:
## arima(x = ibm, order = c(20, 0, 2), fixed = c(NA, rep(0, 18), NA, NA, NA, NA))
## 
## Coefficients:
##          ar1  ar2  ar3  ar4  ar5  ar6  ar7  ar8  ar9  ar10  ar11  ar12
##       0.0353    0    0    0    0    0    0    0    0     0     0     0
## s.e.  0.2191    0    0    0    0    0    0    0    0     0     0     0
##       ar13  ar14  ar15  ar16  ar17  ar18  ar19    ar20      ma1     ma2
##          0     0     0     0     0     0     0  0.0337  -0.0329  -0.030
## s.e.     0     0     0     0     0     0     0  0.0107   0.2192   0.011
##       intercept
##          0.0397
## s.e.     0.0157
## 
## sigma^2 estimated as 2.188:  log likelihood = -16181.68,  aic = 32375.36

ts.diag(mean.model2)

# Estima o modelo
spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,3)),
                   mean.model = list(armaOrder = c(1,2), include.mean = TRUE),
                   fixed.pars = list(beta2 = 0),
                   distribution.model = "std")

fit <- ugarchfit(spec = spec, data = ibm)
fit

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,3)
## Mean Model   : ARFIMA(1,0,2)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.028094    0.012489   2.2496 0.024475
## ar1    -0.695279    0.142834  -4.8677 0.000001
## ma1     0.697011    0.143258   4.8654 0.000001
## ma2    -0.023459    0.012113  -1.9366 0.052790
## omega   0.038366    0.008470   4.5297 0.000006
## alpha1  0.077751    0.009147   8.5004 0.000000
## beta1   0.460590    0.050854   9.0572 0.000000
## beta2   0.000000          NA       NA       NA
## beta3   0.443131    0.047696   9.2907 0.000000
## shape   6.521700    0.407805  15.9922 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.028094    0.012534   2.2414 0.025001
## ar1    -0.695279    0.133882  -5.1932 0.000000
## ma1     0.697011    0.134788   5.1712 0.000000
## ma2    -0.023459    0.012453  -1.8838 0.059591
## omega   0.038366    0.009782   3.9221 0.000088
## alpha1  0.077751    0.010595   7.3383 0.000000
## beta1   0.460590    0.026472  17.3992 0.000000
## beta2   0.000000          NA       NA       NA
## beta3   0.443131    0.026973  16.4288 0.000000
## shape   6.521700    0.467932  13.9373 0.000000
## 
## LogLikelihood : -15200.26 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       3.4033
## Bayes        3.4104
## Shibata      3.4033
## Hannan-Quinn 3.4057
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       3.392 0.06550
## Lag[2*(p+q)+(p+q)-1][8]      5.797 0.02345
## Lag[4*(p+q)+(p+q)-1][14]     7.830 0.40183
## d.o.f=3
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       1.818  0.1776
## Lag[2*(p+q)+(p+q)-1][11]     5.258  0.5417
## Lag[4*(p+q)+(p+q)-1][19]     6.819  0.8011
## d.o.f=4
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[5]   0.02487 0.500 2.000  0.8747
## ARCH Lag[7]   0.30236 1.473 1.746  0.9477
## ARCH Lag[9]   1.08745 2.402 1.619  0.9190
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  5.2814
## Individual Statistics:             
## mu     0.2218
## ar1    0.9004
## ma1    0.8647
## ma2    0.1208
## omega  0.7550
## alpha1 0.1677
## beta1  0.4162
## beta3  0.4389
## shape  0.5626
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.1 2.32 2.82
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias           1.7159 0.0862118   *
## Negative Sign Bias  3.5549 0.0003801 ***
## Positive Sign Bias  0.6037 0.5460820    
## Joint Effect       13.2923 0.0040453 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     38.55     0.005049
## 2    30     54.03     0.003214
## 3    40     61.71     0.011703
## 4    50     77.20     0.006223
## 
## 
## Elapsed time : 2.837802

## 
## please wait...calculating quantiles...

# Computa as previões do modelo
pred <- ugarchforecast(fit)
# Grava o valor dos retornos previstos
(r1 <- pred@forecast$seriesFor[1])

## [1] 0.03496415

(r5 <- pred@forecast$seriesFor[5])

## [1] 0.04100543

# Grava o valor das variâncias previstas
var1 <- pred@forecast$sigmaFor[1]
var5 <- pred@forecast$sigmaFor[5]
# Computa o erro-padrão
(sd1 <- sqrt(var1))

## [1] 1.470796

(sd5 <- sqrt(var5))

## [1] 1.471969

# Computa o quantil 95 (uso a estimativa do modelo para o número de graus de liberdade da t)
(quantil <- qt(.95, 6.51520))

## [1] 1.916005

# VaR
VaR_1 <- r1 - quantil*sd1
VaR_5 <- r5 - quantil*sd5*sqrt(5)
# Em reais
r_VaR_1 <- abs(VaR_1*10^6)
r_VaR_5 <- abs(VaR_5*10^6)

# Organizo os resultados numa tabela
VaR <- data.frame(var = c(VaR_1, VaR_5), var_reais = c(r_VaR_1, r_VaR_5))

# Para facilitar a visualização
class(VaR[,2]) <- c("money", class(VaR[,2]))
print.money <- function(x, ...) {
  print.default(paste0("R$", formatC(as.numeric(x), format="f", digits=0, big.mark=",")))
}
format.money  <- function(x, ...) {
  paste0("R$", formatC(as.numeric(x), format="f", digits=0, big.mark=","))
}
rownames(VaR) <- c("k = 1", "k = 5")
kable(VaR, align = 'c', col.names = c("VaR", "VaR (em reais)"), caption = "Estimação do VaR usando o moddelo ARMA-GARCH")

#### Histograma ####
hist(ibm, freq = F, breaks = 40, xlim = c(-10,10), ylim = c(0,.4), xlab = "", ylab = "Densidade", main = "Distribuição dos retornos")
# Define o 'eixo-x'
X <- seq(-10,10,.05)
# Linha vermelha com a densidade da normal
lines(X, dnorm(X), col = 2)
# Linha azul com a densidade da t-Student(3 g.l.)
lines(X, dt(X, 3), col = 4)
# Legenda
legend("topleft", legend = c("Normal", "t-Student"), col = c(2,4), lty=1)

# Minimiza a distância entre p_i e o valor 0.05
fn <- function(x){abs((x-1/2)/length(ibm) - 0.05)}
# Encontra o valor de x que minimiza a distância acima
optim(par = 1, fn = fn, method = "Brent", lower = 1, upper = length(ibm))$par

## [1] 447.4

# Calcula os valores dos quantis
(p_447 <- round((447-1/2)/length(ibm),5))

## [1] 0.04996

(p_448 <- round((448-1/2)/length(ibm),5))

## [1] 0.05007

# Calcula o peso
(f_447 <- (0.05 - p_447)/(p_448 - p_447))

## [1] 0.3636364

# Ordena a série dos retornos
r <- sort(as.numeric(ibm))
# Computa a combinação linear
(q_05 <- (1-f_447)*r[447] + f_447*r[448])

## [1] -2.148636

# Uso a função do R para conferir o valor encontrado acima
quantil <- quantile(ibm, probs = seq(0.01,1,0.01), type = 8)
quantil[5]

##       5% 
## -2.14875

# As estimativas coincidem até a 3ª casa decimal

# Pacotes adicionais
library(zoo)
library(fExtremes)
# Importo os dados
IBM <- read.table("d-ibmln.txt", quote="\"", comment.char="")
# Transformo em zoo
zibm <- zoo(IBM)

# Encontra o número n de grupos de 21 elementos
n <- floor(length(zibm)/21)
# Remove algumas observações inicias, de maneira a compatibilizar a série com a escolha dos grupos
sibm <- ibm[(length(zibm)-n*21+1):length(zibm)]
# Divide a série em n colunas (cada coluna é um grupo)
cut <- matrix(sibm, ncol = n)
# Encontra o valor máximo em cada coluna (grupo) e grava o resultado num vetor
maximo <- apply(cut, 2, FUN = max)

out <- gevFit(maximo, type = "mle", method = "BFGS")
summary(out)

## 
## Title:
##  GEV Parameter Estimation 
## 
## Call:
##  gevFit(x = maximo, type = "mle", method = "BFGS")
## 
## Estimation Type:
##   gev mle 
## 
## Estimated Parameters:
##        xi        mu      beta 
## 0.1700723 2.1791886 0.9187072 
## 
## Standard Deviations:
##          xi         mu       beta 
## 0.03671146 0.05007071 0.03892154 
## 
## Log-Likelihood Value:
##   675.7119 
## 
## Type of Convergence:
##   0 
## 
## Description
##   Mon Nov 12 17:08:45 2018

# Gravo os valores estimados dos parâmetros
xi <- as.numeric(out@fit$par.ests[1])
mu <- as.numeric(out@fit$par.ests[2])
sigma <- as.numeric(out@fit$par.ests[3])
ses <- out@fit$par.ses

p = 0.05

VaR_tve <- mu - sigma/xi*(1-(-n*log(1-p))^(-xi))

VaR_tve_reais <- 10^6*VaR_tve


VaR_tve_5 <- 5^(xi)*VaR_tve

VaR_tve_reais_5 <- 10^6*VaR_tve_5

# Organizo os resultados numa tabela
VaR <- data.frame(var = c(VaR_tve, VaR_tve_5), var_reais = c(VaR_quant, VaR_quant5))
rownames(VaR) <- c("k = 1", "k = 5")