Introduction to R for Quantitative Finance - Chapter 1
Marcus Marinho
30/04/2021
R Markdown
Instalando a biblioteca
#install.packages("zoo")Limpando as variáveis:
rm(list = ls())Setando o espaço de trabalho e instalando as bibliotecas:
setwd("C:/Users/onlym/OneDrive/Documentos/ProgetosR/Chapter_1")
library("zoo")##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericCarregando o arquivo e plotando o gráfico:
aapl <- read.zoo("aapl.csv", sep=",", header=TRUE,format= "%Y-%m-%d")plot(aapl, main = "APPLE Closing Prices", ylab = "Price USD", xlab = "Date")
Observando o início e o final da time série:
head(aapl)## 2000-01-03 2000-01-04 2000-01-05 2000-01-06 2000-01-07 2000-01-10
## 27.58 25.25 25.62 23.40 24.51 24.08tail(aapl)## 2013-04-17 2013-04-18 2013-04-19 2013-04-22 2013-04-23 2013-04-24
## 402.80 392.05 390.53 398.67 406.13 405.46Exibindo o maior preço histórico:
aapl[which.max(aapl)]## 2012-09-19
## 694.86Ao lidar com séries temporais, normalmente estamos mais interessados em retornos do que preços. Isso ocorre porque os retornos são estacionários. Calcularemos retornos simples e retornos compostos continuamente (em termos percentuais)
ret_simples <- diff(aapl) / lag(aapl, k = -1) * 100
ret_cont <- diff(log(aapl)) * 100Resumo de estatísticas do retorno simples também podem ser obtidas. Nos usamos o método coredatapara indicar que nós queremos os preços da ação e não o indice (data)
summary(coredata(ret_simples))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -51.85888 -1.32475 0.07901 0.12527 1.55332 13.91131O maior dia de prejuízo foi de -51,86%. A data em que isso ocorreu pode ser obtida assim:
ret_simples[which.min(ret_simples)]## 2000-09-29
## -51.85888Para obter uma melhor compreensão das frequencias relativas dos retornos diários, podemos traçar o histograma. O número de células usadas para agrupar os dados de retorno pode ser especificado usando o argumento break
hist(ret_simples, breaks = 100, main = "Histograma do retorno simples", xlab = "%")
Podemos restringir nossa análise a uma janela dos dados
aapl_2013 <- window(aapl, start = '2013-01-01', end = '2013-12-31')
aapl_2013[which.max(aapl_2013)]## 2013-01-02
## 545.85aapl_2013[which.min(aapl_2013)]## 2013-04-19
## 390.53Cálculo do VAR( Value at Risk). Trata-se da medida de perda percentual de uma carteira de investimentos sujeita aos riscos de mercado. Com o cálculo do VAR é possível obter o valor esperado da máxima perda(ou pior perda) em um horizonte de tempo com um intervalo de confiança.
quantile(ret_simples, probs = 0.01)## 1%
## -7.042678Isto significa que existe 1% de probabilidade de que o pior retorno de um dado dia seja negativo em -7%. Para um nível de significância de 5%, basta alterar a fórmula.
quantile(ret_simples, probs = 0.05)## 5%
## -4.183982Ao nível de significância de 5%, podemos dizer que em um dado dia o pior retorno seria de -4%.
ARIMA -Autoregressive Integrated Moving Average
Primeiro vamos instalar o pacote forecaste depois fazer upload dos arquivo com os preços mensais das casas UK
#install.packages("forecast")Carregandoo o forecast
library('forecast')## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoohp <- read.zoo("UKHP.csv", sep = ",",header= TRUE,format = "%Y-%m", FUN = as.yearmon)O argumento FUN aplicado a função(as.yearmon, que representa pontos mensais de dados) a coluna data. Para ter certeza que foi guardado os dados mensais (12 subperíodos por período) ao especificar as.yearmonnós pesquisamos a frequencia da série de dados.
frequency(hp)## [1] 12O resultado significa que temos doze subperíodos(meses) por período (anos). Nós novamente usamos retorno simples em nossa análise.
hp_ret <- diff(hp) / lag(hp, k = -1) * 100 MODEL IDENTIFICATION AND ESTIMATION
Usamos auto.arimado pacote forecastpara identificar o modelo ótimo e estimar os coeficientes em um passo. A função recebe alguns argumentos além do retorno das séries (hp_ret). Especificando stationary=TRUE, nós restringimos a pesquisa por modelos estacionários. De modo similar, seasonal=FALSE restringe a bisca por modelos não sazonais.
mod <- auto.arima(hp_ret, stationary = TRUE, seasonal = FALSE, ic="aic")Para determinar o coeficiente de valores ajustados, nós chamamos o output
mod## Series: hp_ret
## ARIMA(2,0,0) with non-zero mean
##
## Coefficients:
## ar1 ar2 mean
## 0.2299 0.3491 0.4345
## s.e. 0.0573 0.0575 0.1519
##
## sigma^2 estimated as 1.118: log likelihood=-390.97
## AIC=789.94 AICc=790.1 BIC=804.28Um processo ar2 parece se adequar melhor aos dados, de acordo com os critérios de informação de akike. Para confirmação visual, podemos plotar a função de autocorrelação parcial usando o comando pacf. Ele mostra autocorrelações parciais diferentes de zero até o atraso dois, portanto, um processo AR de ordem dois parece ser apropriado. Os dois coeficientes AR, a interceptação(que é realmente a média dse o modelo contiver um termo AR) e a respectivos erros padrão são fornecidos. No exemplo a seguir, eles são todos significativos no nível de 5%, uma vez que os respectivos intervalos de confiança não contém zero.
confint(mod)## 2.5 % 97.5 %
## ar1 0.1174881 0.3422486
## ar2 0.2364347 0.4617421
## intercept 0.1368785 0.7321623Um rápido caminho para validar o modelo é plotar o diagnóstico da TS usando o seguinte comando:
tsdiag(mod)
#### Nosso modelo parece bom, uma vez que os resíduos padronizados não mostram clusters de volatilidade, nenhuma autocorrelação significativa entre os resíduos de acordo com o gráfico ACF, e o teste Ljung-Box para autocorrelação mostra altos valores de p, de modo que a hipótese nula de resíduos independentes não pode ser rejeitada.Para avaliar o quão bem o modelo representa os dados na amostra, podemos plotar os retornos mensais brutos (a linha sólida preta fina) versus os valores ajustados (a linha pontilhada vermelha espessa).
plot(mod$x, lty = 1, main = "UK house prices: raw data vs. fitted values", ylab = "Return in percent", xlab = "Date")
lines(fitted(mod), lty = 2, lwd = 2, col = "red")
Além disso, podemos calcular medidas comuns de precisão
accuracy(mod)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.001203096 1.051422 0.8059929 -Inf Inf 0.7154503 -0.004848183This command returns the mean error, root mean squared error, mean absolute error, mean percentage error, mean absolute percentage error, and mean absolute scaled error.
Para a predição dos retornos dos próximos três meses:
predict(mod, n.ahead = 3)## $pred
## Apr May Jun
## 2013 0.5490544 0.7367277 0.5439708
##
## $se
## Apr May Jun
## 2013 1.057401 1.084978 1.165247Portanto, esperamos um ligeiro aumento nos preços médios das residências nos próximos três meses, mas com um erro padrão elevado de cerca de 1%. Para plotar a previsão com erros padrão, podemos usar o seguinte comando:
plot(forecast(mod))
COINTEGRAÇÃO
A ideia por trás da cointegração, um conceito introduzido por Granger (1981) e formalizado por Engle e Granger (1987), é encontrar uma combinação linear entre séries temporais não estacionárias que resultam em séries temporais estacionárias. Portanto, é possível detectar relações estáveis de longo prazo entre séries temporais não estacionárias (por exemplo, preços).
Primeiro carregamos as bibliotecas necessárias. A biblioteca urca tem alguns métodos úteis para testes de raiz unitária e para estimar relações de cointegração
#install.packages("urca")library("urca")Vamos importar o preço mensal do combustível de avião e heating oil (usd por galão)
prices <- read.zoo("JetFuelHedging.csv", sep = ",",FUN = as.yearmon, format = "%Y-%m", header = TRUE)Levando em consideração apenas o comportamento de curto prazo (variações mensais de preços) das duas commodities, pode-se derivar o hedge de variância mínima ajustando um modelo linear que explica as variações nos preços do combustível de aviação pelas variações nos preços do óleo para aquecimento. O coeficiente beta dessa regressão é a proporção de hedge ideal.
simple_mod <- lm(diff(prices$JetFuel) ~ diff(prices$HeatingOil)+0)
summary(simple_mod)##
## Call:
## lm(formula = diff(prices$JetFuel) ~ diff(prices$HeatingOil) +
## 0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52503 -0.02968 0.00131 0.03237 0.39602
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## diff(prices$HeatingOil) 0.89059 0.03983 22.36 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0846 on 189 degrees of freedom
## Multiple R-squared: 0.7257, Adjusted R-squared: 0.7242
## F-statistic: 499.9 on 1 and 189 DF, p-value: < 2.2e-16Obtemos uma razão de hedge de 0,89059 e um erro padrão residual de 0,0846. A sebe cruzada não é perfeita; a carteira protegida resultante ainda é arriscada.
Agora tentamos melhorar essa relação de hedge usando uma relação de longo prazo existente entre os níveis de combustível de aviação e os preços futuros do óleo para aquecimento. Você já pode adivinhar a existência de tal relação traçando as duas séries de preços (os preços do óleo para aquecimento estarão em vermelho) usando o seguinte comando:
plot(prices$JetFuel, main = "Jet Fuel and Heating Oil Prices",xlab = "Date", ylab = "USD")
lines(prices$HeatingOil, col = "red")
Usamos a técnica de estimativa em duas etapas de Engle e Granger. Em primeiro lugar, ambas as séries temporais são testadas para uma raiz unitária (não tacionalidade) usando o teste Dickey-Fuller aumentado.
jf_adf <- ur.df(prices$JetFuel, type = "drift")
summary(jf_adf)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.06212 -0.05015 0.00566 0.07922 0.38086
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.03050 0.02177 1.401 0.16283
## z.lag.1 -0.01441 0.01271 -1.134 0.25845
## z.diff.lag 0.19471 0.07250 2.686 0.00789 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.159 on 186 degrees of freedom
## Multiple R-squared: 0.04099, Adjusted R-squared: 0.03067
## F-statistic: 3.975 on 2 and 186 DF, p-value: 0.0204
##
##
## Value of test-statistic is: -1.1335 0.9865
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81A hipótese nula de não estacionariedade (a série temporal do combustível de aviação contém uma raiz unitária) não pode ser rejeitada ao nível de significância de 1%, pois a estatística de teste de -1,1335 não é mais negativa do que o valor crítico de -3,46. O mesmo se aplica aos preços do óleo para aquecimento (a estatística de teste é -1,041).
ho_adf <- ur.df(prices$HeatingOil, type = "drift")
summary(ho_adf)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.78839 -0.06344 -0.00128 0.07418 0.49186
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.02815 0.02115 1.331 0.1847
## z.lag.1 -0.01314 0.01263 -1.041 0.2992
## z.diff.lag 0.14419 0.07296 1.976 0.0496 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1534 on 186 degrees of freedom
## Multiple R-squared: 0.02418, Adjusted R-squared: 0.01369
## F-statistic: 2.304 on 2 and 186 DF, p-value: 0.1027
##
##
## Value of test-statistic is: -1.041 0.9002
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81Podemos agora prosseguir para estimar o modelo de equilíbrio estático e testar os resíduos para uma série temporal estacionária usando um teste Dickey-Fuller aumentado. Observe que diferentes valores críticos [por exemplo, de Engle e Yoo (1987)] devem ser usados agora, uma vez que a série sob investigação é estimada.
mod_static <- summary(lm(prices$JetFuel ~ prices$HeatingOil))
error <- residuals(mod_static)
error_cadf <- ur.df(error, type = "none")
summary(error_cadf)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19162 -0.03385 -0.00516 0.02879 0.47426
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.74476 0.08357 -8.912 4.46e-16 ***
## z.diff.lag 0.12304 0.07264 1.694 0.092 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07107 on 187 degrees of freedom
## Multiple R-squared: 0.3417, Adjusted R-squared: 0.3347
## F-statistic: 48.53 on 2 and 187 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -8.912
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62A estatística de teste obtida é -8,912 e o valor crítico para um tamanho de amostra de 200 ao nível de 1% é -4,00; portanto, rejeitamos a hipótese nula de não estacionariedade. Assim, descobrimos duas variáveis cointegradas e podemos prosseguir com a segunda etapa; ou seja, a especificação de um modelo de correção de erros (ECM). O ECM representa um modelo dinâmico de como (e com que rapidez) o sistema retorna ao equilíbrio estático estimado anteriormente e é armazenado na variável mod_static.
djf <- diff(prices$JetFuel)
dho <- diff(prices$HeatingOil)
error_lag <- lag(error, k = -1)
mod_ecm <- lm(djf ~ dho + error_lag)
summary(mod_ecm)##
## Call:
## lm(formula = djf ~ dho + error_lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19332 -0.03391 -0.00101 0.02128 0.44942
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.001511 0.005013 0.301 0.763
## dho 0.899489 0.032547 27.637 <2e-16 ***
## error_lag -0.655301 0.066303 -9.883 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06892 on 187 degrees of freedom
## Multiple R-squared: 0.8189, Adjusted R-squared: 0.817
## F-statistic: 422.9 on 2 and 187 DF, p-value: < 2.2e-16