Introdução

Neste curso discutiremos sobre análise de séries temporais e, basicamente, estaremos interessados em fazer boas previsões. Mas, o que é uma boa previsão? Qual modelo devemos usar? Como muito bem disse Box, “todos os modelos estão errados, mas alguns são úteis”. Esse pensamento de Box nunca esteve tão atual. Em uma época em que as mudanças estão ocorrendo em uma velocidade imensa, acertar as previsões está cada dia mais difícil. O World Economic Forum publicou recentemente um texto interessante [We’re moving fast. But nobody knows where we’re going], no qual os autores fazem uma pergunta bem provocativa: “The end of forecasting?”. A ideia do artigo é simples e pertinente: em uma era onde as coisas se modificam em uma velocidade incrível, qual seria a utilidade e a veracidade das previsões?


The sooner we realize that long-term forecasting is becoming obsolete, the better we’ll be able to cope with the new reality. Businesses are already responding by hiring **Chief Digital Officers**, promoting entrepreneurial culture, and elevating the innovation imperative to the level of existential priority.[World Economic Forum]


Hal Varian (Varian 2014), economista chefe da Alphabet, escreveu um artigo super interessante, e que eu recomendo fortemente a leitura para meus alunos, sobre como o estudo da econometria deve evoluir na nova era do “Big Data”. Em suma, mais que boas previsões, é preciso modificar a cultura das empresas e capacitar pessoas para que estejam apitas a viver a nova realidade.


Nesse curso, nós vamos ajudá-lo a adquirir uma importante ferramenta de um cientista de dados, que é a capacidade de criar cenários e estar apto a responder à questão: “E se isso acontecer…”. Vamos explorar os modelos univariados e multivariados; calcular elasticidades, fazer previsões dentro e fora da amostra e criar cenários baseados em diferentes conjunturas econômicas.


É claro que este é um curso introdutório, há avanços importantes no mundo de séries de tempo que não veremos nesse curso, mas que temos pesquisado e implementado para a nossa realidade. Por exemplo, a economista Lucrezia Reichlin lidera um time de pesquisadores na empresa Now-casting.com onde o objetivo é fazer previsões diárias do PIB de diversos países. Estes modelos, além nos permitir utilizar séries temporais com múltiplas frequências (e.g. mensais, trimestrais, diárias), permite também a “extração do sentimento” de texto jornalísticos para atualizar as previsões, neste caso, usando machine learning.


Por fim, no meu site é possível acompanhar algumas das pesquisas que estamos fazendo e que utilizam modelos de séries de tempo um pouco mais complexos do que os modelos que veremos ao longo desse curso. Além disso, nesse site, disponibilizo uma versão preliminar do livro de Séries Temporais que fiz com o time do NMEC|FGV e que será muito útil para o acompanhamento desse curso.


O curso é dividido, basicamente, em três grandes blocos: (i) introdução a séries temporais [material disponibilizado nesse arquivo]; (ii) Modelos Univariados [Box & Jenkins] e (iii) Modelos multivariados [Autoregressive Distributed Lag]. Nessa primeira parte introdutória, abordaremos os seguintes tópicos:

  1. Explorar os softwares R e RStudio

    Vamos discutir nossas motivações para utilizar o software R e fazer uma breve introdução.

  2. Análise de Séries Temporais

    Nessa seção iremos iniciar nossa discussão sobre séries temporais. A ideia é explorar alguns pacotes do R utilizados para lidar com séries de tempo e ainda, discutir conceitos importantes como, sazonalidade, ciclo e tendência. Vamos também, mostrar uma estratégia de trading utilizando pacotes e conceitos de séries temporais.

  3. Estudo de caso: trading Itaú e Bradesco

    vamos aplicar a teoria da cointegração para definir uma estratégia de operação na bolsa de valores. A ideia é simples: se duas séries tem dependência de longo prazo e, no curto prazo uma série temporal se descola da outra, espera-se que, após algum tempo, haverá uma reversão para a outra série temporal.

  4. Estudo de caso: venda de sorvete

    Quais são as variáveis importantes para determinar a venda de sorvete? O que aconteceria com a venda de sorvetes se a temperatura cair? Por que é importante prever a venda de sorvetes? São essas e outras questões que discutiremos nesse tópico.

Explorando os softwares R e RStudio

Por que escolher o R?

  1. Programadores em R estão entre os maiores salários pagos em 2014

  2. R is Still Hot - and Getting Hotter

  3. R is Ready for Business

  4. R #5 in IEEE 2016 Top Programming Languages

IEEE 2016 Top Programming Languages

IEEE 2016 Top Programming Languages

  1. SWIRL
  2. SWIRL

R package makes it fun and easy to learn R programming and data science.

install.packages("swirl")
library("swirl")
swirl()


R software

O R não é apenas uma linguagem de programação, é também um ambiente integrado de desenvolvimento. É um projeto GNU, que foi desenvolvido em Bell Laboratories (antiga AT & T, agora Lucent Technologies) por John Chambers e seus colegas. (Team (2017))

De acordo com seus criadores, o R é um conjunto integrado de instalações de software para manipulação de dados, cálculo e exibição gráfica, e esse ambiente inclui:

  • Instalação eficaz de tratamento e armazenamento de dados;
  • Conjunto de operadores para cálculos em numéricos, em vetores e matrizes;
  • Grande coleção coerente e integrada de ferramentas intermediárias para análise de dados;
  • Instalações gráficas para análise de dados e exibição na tela ou em console;
  • Linguagem de programação bem desenvolvida, simples e eficaz que inclui condicionais, loops, funções recursivas definidas pelo usuário e instalações de entrada e saída.

Além disso, o R conta com uma ampla comunidade acadêmica que oferece suporte e ajuda aos usuários (CRAN) além de oferecer um extensão às funções usuais do R, com os pacotes.

Isso faz do R um projeto aberto, o que sigfinica ser continuamente melhorado, atualizado e expandido pela comunidade global de desenvolvedor e seus uários incrivelmente apaixonados. (Bowles 2015)

RStudio software

O RStudio é um ambiente de desenvolvimento integrado (IDE) para R. Nele estão inclusos ferramentas que facilitam o uso da linguagem R, como um console e uma janela que suporta a execução do código direto, e instrumentos que auxiliam na análise histórica do script e debugação de possíveis erros.

O RStudio está disponível em open source e edições comerciais e é executado no desktop (Windows, Mac e Linux) ou em um navegador conectado ao RStudio Server ou ao RStudio Server Pro (Debian / Ubuntu, RedHat / CentOS e SUSE Linux).

Instalando os softwares

É muito fácil instalar os softwares R e RStudio. Abaixo nós mostramos os passos para a instalação, mas, caso tenha alguma dúvida, assista a este vídeo no youtube: Instalação dos softwares R e RStudio

Instalando o R

  1. Vá no site cran.r-project.org
  2. Faça download do arquivo
  3. Clique em em executar

Instalando o RStudio

  1. Vá no site rstudio.com
  2. Clique em dowload
  3. Escolha o free
  4. Faça o download do Arquivo
  5. Clique em executar

Primeiros passos

Vamos alguns exemplos de como fazer operações usando o R. Para maiores detalhes, consulte o primeiro capítulo do livro Análise de Séries Temporais em R que está disponível no meu site.

Exemplo 1: Qual é o resultado de \(\sqrt[2]{25}\)?

sqrt(25)
## [1] 5


Exemplo 2: Qual é o resto da divisão de \(\frac{5}{2}\)? 

5%%2
## [1] 1


Exemplo 3: Tome \(a=25\) e \(b=45\) analise se “a” é maior:

a=25
b=45

a>b
## [1] FALSE


Exemplo 4: Tome \(a=15\) e \(b=25\) faça \(c=a+b\) e \(d=\left ( \sqrt[3]{c}\right )\):

a=15
b=25
(c=a+b)
## [1] 40
d=c^(1/3)
d
## [1] 3.419952


Exemplo 5: Tome \(nome_1 = Big\) e \(nome_2 = Data\). Faça \(nomecompleto= nome_1 +nome_2\):

Para resolver essa questão é necessário fazer uso da função paste(), essa função é usada para concatenar duas ou mais variáveis character.

nome1="big"
nome2='data'
nomecompleto=paste(nome1,nome2)
nomecompleto
## [1] "big data"
## Caso não seja útil o espaço entre as palavras, basta usar a função paste0().


Exemplo 6: Transforme 06/02/2017 em data:

data="06/02/2017"
formatodata = as.Date(x = data,format="%d/%m/%Y")
formatodata
## [1] "2017-02-06"


Exemplo 7: Concatene as listas com os nomes dos alunos:

nomes1<-c("Anna","Pedro","Carlos","Bruno","Vanessa","Paula","Italo")
nomes2<-c("Jorge","Davi","Mariana","Carolina","Alice")

(nomes_completo <- c(nomes1,nomes2))
##  [1] "Anna"     "Pedro"    "Carlos"   "Bruno"    "Vanessa"  "Paula"   
##  [7] "Italo"    "Jorge"    "Davi"     "Mariana"  "Carolina" "Alice"


Exemplo 8: Veja a tabela criada abaixo. a) Acesse o preço do Produto D na tabela produtos; b) Acesse os preço do Produto D e Produto E na tabela produtos; c) Acesse a coluna produtos; d) Crie uma coluna quantidade, na tabela preços, que receba os valores \(quantidade<-c(50,100,120,150,200)\):

Produto<-c("Produto A"," Produto B","Produto C", "Produto D","Produto E")
Preco<-c(5,15,4,6,8)
tabela_preco_produto<-data.frame(Produto,Preco)
tabela_preco_produto
##      Produto Preco
## 1  Produto A     5
## 2  Produto B    15
## 3  Produto C     4
## 4  Produto D     6
## 5  Produto E     8
## Acesse o preço do *Produto D* na tabela produtos: 
tabela_preco_produto[4,2]
## [1] 6
## Acesse os preço do *Produto D* e *Produto E* na tabela produtos:
tabela_preco_produto[4:5,2]
## [1] 6 8
## Acesse a coluna produtos
tabela_preco_produto[,"Produto"]
## [1] Produto A   Produto B Produto C  Produto D  Produto E 
## Levels:  Produto B Produto A Produto C Produto D Produto E
## Crie uma coluna quantidade, na tabela preços, que receba os valores:

tabela_preco_produto$quantidade<-c(50,100,120,150,200)
tabela_preco_produto
##      Produto Preco quantidade
## 1  Produto A     5         50
## 2  Produto B    15        100
## 3  Produto C     4        120
## 4  Produto D     6        150
## 5  Produto E     8        200


Exemplo 9: Seja a matriz A. Qual é o valor na primeira linha e segunda coluna da matriz A?

A<-matrix(c(0,7,0,1,4,5,0,7,0),ncol=3,nrow = 3,byrow = TRUE)
A
##      [,1] [,2] [,3]
## [1,]    0    7    0
## [2,]    1    4    5
## [3,]    0    7    0
A[1,2]
## [1] 7


Exemplo 10: Crie uma lista com os parâmetros:

  • Nome do Curso: Curso de Formação Executiva em Big Data e Data Science
  • Turma: 3
  • Aula: Introdução ao R
minha_lista<-list("Curso de Formação Executiva em Big Data e Data Science",3, "Introdução ao R")
minha_lista
## [[1]]
## [1] "Curso de Formação Executiva em Big Data e Data Science"
## 
## [[2]]
## [1] 3
## 
## [[3]]
## [1] "Introdução ao R"


Exemplo 11: Tome \(a=31\) e \(b=30\), analise qual dos dois objetos é o maior.

a=31
b=30
if(a<b){
   cat("a é maior que b")
 }else{
   cat("b é maior que a")
 }
## b é maior que a


Exemplo 12: Tome o vetor \(v=[16,18,59,35,27,37,38]\), faça uma estrutura que imprima os valores do vetor.

v=c(16,18,59,35,27,37,38)
for (i in 1:length(v)){
  print(v[i])
  }
## [1] 16
## [1] 18
## [1] 59
## [1] 35
## [1] 27
## [1] 37
## [1] 38


Exemplo 13: Construa uma função que retornará a média de dois números.

media<-function(a,b){
   m<-(a+b)/2
   return(m)
 }

media(3,5)
## [1] 4


—

Análise de Séries Temporais


Vamos introduzir a ideia de análise de séries temporais e abordar alguns conceitos importantes sobre o tema, como: sazonalidade, estacionariedade, heterocedásticidade e ciclo.


Exemplos de séries de tempo


1. Séries Temporais em diferentes frequências

## Daily morning gold prices
library(ggfortify)
library(fpp2)
autoplot(gold)

## Quarterly production of woollen yarn in Australia: tonnes
autoplot(woolyrnq)

## Australian monthly gas production
autoplot(gas)

## Half-hourly electricity demand
autoplot(taylor)


2. Ondas de rádio


3. Eletrocardiograma


4. curvas de Carga


5. Sazonalidade ou Ciclo

  • Exemplo de um evento cíclico

O ciclo de de vida da população de lebres e de seu principal predador, o lince-canadense, foi revelado nos registros do número de peles vendidas por caçadores para a Hudson Bay Company.

## Número de armadilhas para capturar Lince no Canadá
library(ggfortify)
data(lynx); trapped <- lynx
autoplot(lynx)


A ST apresenta tendência?

trapped.time <- time(trapped)
trapped.lm <- lm(trapped ~ trapped.time)
summary(trapped.lm)
## 
## Call:
## lm(formula = trapped ~ trapped.time)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -1594  -1211   -755   1032   5366 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -4630.034   8493.112  -0.545    0.587
## trapped.time     3.285      4.523   0.726    0.469
## 
## Residual standard error: 1589 on 112 degrees of freedom
## Multiple R-squared:  0.004689,   Adjusted R-squared:  -0.004198 
## F-statistic: 0.5276 on 1 and 112 DF,  p-value: 0.4691

De acordo com o modelo, houve um aumento de 3,39 armadilhas por ano, contudo esse parâmetro não foi significante, indicando que não há tendência.


Outra informação que gostaríamos de saber é o tamanho do ciclo. Como observamos graficamente, parece que o ciclo de vida desses animais é de aproximadamente 5 anos.

ma10 <- filter(x=trapped, filter=rep(x=1/10,times=10), sides=2)
ma5 <- filter(x=trapped, filter=rep(x=1/5,times=5), sides=2)
plot(trapped,col="grey", ylab = "Lynx Trapped", xlab = "Year")
lines(ma10,col="red",lwd=2)
lines(ma5,col="purple",lwd=2)
legend("topright", c(" 5-Year Moving Average", "10-Year Moving Average"), bty = "n", col = c("purple","red"), pch = 15)
abline(trapped.lm, col="black",lwd=2, lty="dashed")


* Exemplo de um evento sazonal

library(BETS)
# Vendas de autoveículos no mercado interno 
veiculos <- BETS.get(code = 1379)
library(fpp2)
autoplot(veiculos)

monthplot(veiculos)

ggseasonplot(veiculos, polar = T)


Alguns conceitos importantes


1. Caracterização, Modelagem e Previsão de uma série temporal

  • Uma série temporal (ST) é um conjunto de observações, geralmente equiespaçadas, obtidas a partir da medição/observação de uma variável ao longo do tempo;
  • A frequência de observação/medição de uma ST pode variar dependendo do fenômeno observado: minuto, diária, semanal, mensal, anual etc;
  • A partir da análise de ST, é possível obter subsídios para a escolha de um modelo adequado para modelar a série, escolhido dentro de uma classe de modelos pré-existentes;
  • Uma vez construído, um modelo de ST pode ser utilizado para efetuar previsões probabilísticas sobre o futuro da série e/ou analisar a sensibilidade de algumas variáveis à variável objetivo;
  • A capacidade de realizar previsões é fundamental no processo de tomada de decisões em diversos contextos e em diversos lugares como orgãos públicos e empresas.


2. Função de autocorrelação

  • É uma medida de “memória” do processo;
  • A FAC é o procedimento padrão para investigar a dependência linear implicada por um modelo de ST;
  • A FAC ou correlograma de um determinado PE é obtida através do cálculo da seguinte expressão:
    \[\rho(k)= \frac{\gamma_k}{\gamma_0}\]

Exemplo de uma FAC gerada por um processo AR(1)

require(tseries)
require(BETS)

# Simulando um processo AR(1) com parâmetro 0.5
set.seed(10)
sim.AR<-arima.sim(model=list(ar=0.6),n=100)
library(fpp2)
autoplot(sim.AR)

BETS.corrgram(sim.AR,type = "correlation",lag.max = 36)
BETS.corrgram(sim.AR,type = "partial",lag.max = 36)

3. Processo Ruído Branco

  • RB – é um PE estacionário onde as observações são descorrelatadas, isto é, o processo não possui “memória” linear;
  • Geralmente representado por \(\varepsilon_t\) . É importante na construção de modelos para processos com dependência;

Exemplo de uma FAC gerada por um processo ruído branco

require(tseries)

# Simulando um processo Ruído Branco
set.seed(11)
sim.wn<-arima.sim(model=list(),n=100)
library(fpp2)
autoplot(sim.wn)

acf(sim.wn)

require(BETS)
BETS.corrgram(sim.wn,type = "correlation",lag.max = 36)
BETS.corrgram(sim.wn,type = "partial",lag.max = 36)


4. Intervalo de confiança

Para observarmos o intervalo de confiança vamos utilizar dois modelos bem simples: Naive e SNaive.

library(fpp2)
library(forecast)
# Daily closing stock prices of Google Inc
fcgoog <- naive(goog, h =  20)
autoplot(fcgoog)

summary(fcgoog)
## 
## Forecast method: Naive method
## 
## Model Information:
## Call: naive(y = goog, h = 20) 
## 
## Residual sd: 8.9145 
## 
## Error measures:
##                     ME     RMSE      MAE        MPE      MAPE MASE
## Training set 0.4436236 8.921089 6.008889 0.06493981 0.9815741    1
##                    ACF1
## Training set 0.04680557
## 
## Forecasts:
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1001         838.96 827.5272 850.3928 821.4750 856.4450
## 1002         838.96 822.7915 855.1285 814.2325 863.6875
## 1003         838.96 819.1577 858.7623 808.6751 869.2449
## 1004         838.96 816.0943 861.8257 803.9900 873.9300
## 1005         838.96 813.3954 864.5246 799.8623 878.0577
## 1006         838.96 810.9554 866.9646 796.1306 881.7894
## 1007         838.96 808.7116 869.2084 792.6990 885.2210
## 1008         838.96 806.6231 871.2969 789.5049 888.4151
## 1009         838.96 804.6615 873.2585 786.5050 891.4150
## 1010         838.96 802.8062 875.1138 783.6675 894.2525
## 1011         838.96 801.0416 876.8784 780.9688 896.9512
## 1012         838.96 799.3555 878.5645 778.3901 899.5299
## 1013         838.96 797.7383 880.1817 775.9169 902.0031
## 1014         838.96 796.1822 881.7378 773.5371 904.3829
## 1015         838.96 794.6808 883.2392 771.2408 906.6792
## 1016         838.96 793.2287 884.6913 769.0200 908.9000
## 1017         838.96 791.8212 886.0988 766.8674 911.0526
## 1018         838.96 790.4546 887.4654 764.7774 913.1426
## 1019         838.96 789.1254 888.7946 762.7446 915.1754
## 1020         838.96 787.8308 890.0892 760.7646 917.1554
# Quarterly Australian Beer production
fcbeer <- snaive(ausbeer,h = 16)
autoplot(fcbeer)

summary(fcbeer)
## 
## Forecast method: Seasonal naive method
## 
## Model Information:
## Call: snaive(y = ausbeer, h = 16) 
## 
## Residual sd: 19.1207 
## 
## Error measures:
##                    ME     RMSE      MAE      MPE    MAPE MASE       ACF1
## Training set 3.098131 19.32591 15.50935 0.838741 3.69567    1 0.01093868
## 
## Forecasts:
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2010 Q3            419 394.2329 443.7671 381.1219 456.8781
## 2010 Q4            488 463.2329 512.7671 450.1219 525.8781
## 2011 Q1            414 389.2329 438.7671 376.1219 451.8781
## 2011 Q2            374 349.2329 398.7671 336.1219 411.8781
## 2011 Q3            419 383.9740 454.0260 365.4323 472.5677
## 2011 Q4            488 452.9740 523.0260 434.4323 541.5677
## 2012 Q1            414 378.9740 449.0260 360.4323 467.5677
## 2012 Q2            374 338.9740 409.0260 320.4323 427.5677
## 2012 Q3            419 376.1020 461.8980 353.3932 484.6068
## 2012 Q4            488 445.1020 530.8980 422.3932 553.6068
## 2013 Q1            414 371.1020 456.8980 348.3932 479.6068
## 2013 Q2            374 331.1020 416.8980 308.3932 439.6068
## 2013 Q3            419 369.4657 468.5343 343.2438 494.7562
## 2013 Q4            488 438.4657 537.5343 412.2438 563.7562
## 2014 Q1            414 364.4657 463.5343 338.2438 489.7562
## 2014 Q2            374 324.4657 423.5343 298.2438 449.7562


5. Valores ajustados e resíduos

O valor ajustado é a previsão de uma observação usando todas as observações passadas, isto é:

  1. são previsões 1-passo-a-frente
  2. não é a previsão verdadeira, uma vez que os parâmetros são estimados a cada novo dado
fc <- naive(oil)
autoplot(oil, series = "Dados reais") + xlab("Year") + autolayer(fitted(fc), series = "valores ajustados") + ggtitle("Produção de petróleo na Saudi Arabia")


O resíduo é a diferença entre o valor ajustado/previsto e o valor real [são erros 1-passo-a-frente]

autoplot(residuals(fc))


Os resíduos devem se parecer com um ruído branco:

  1. Características essenciais
    • eles devem ser não correlacionados;
    • devem ter média zero.
  2. Propriedades úteis (para previsão e intervalo de confiança dos parâmetros)
    • devem ter variância constante;
    • normalmente distribuidos.


Uma maneira fácil de testar essas propriedades é usando a função “checkresiduals()”.

checkresiduals(fc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Naive method
## Q* = 12.59, df = 10, p-value = 0.2475
## 
## Model df: 0.   Total lags used: 10


6. Treinamento e Teste de previsão


observações:

  1. O conjunto de testes não deve ser usado para fazer as previsões;
  2. As previsões devem ser construídas usando o conjunto de treinamento;
  3. Um bom modelo para os dados de treinamento, não necessariamente, é um bom modelo para conjunto de teste.
training <- window(oil, end = 2003)
test <- window(oil, start = 2004)
fc <- naive(training, h = 10)
autoplot(fc) + autolayer(test, series = "conjunto de teste") + autolayer(fitted(fc), series = "valores ajustados")


!! Cuidado
Erro de previsão é a diferença entre o valor observado e o valor previsto no conjunto de teste. Resíduos são os erros do conjunto de treinamento, normalmente 1-passo-a-frente (o erro de previsão pode ser múltiplos passos-a-frente)


7. Acurácia dos modelos de previsão

require(forecast)
accuracy(fc,test)
##                    ME     RMSE      MAE      MPE      MAPE      MASE
## Training set  9.87358 52.56156 39.42504 2.506565 12.570647 1.0000000
## Test set     21.60250 35.09832 29.97666 3.963914  5.777875 0.7603458
##                   ACF1 Theil's U
## Training set 0.1801528        NA
## Test set     0.4029519  1.184862

8. Out-of-sample rolling evaluation

com os mesmos parâmetros estimados para os dados passados, move-se a origem de previsão no período out-of-sample, fazendo-se previsões para cada origem;



Estudo de caso: venda de Sorvete

Imagine que você trabalhe na Kibon e que um dos C levels da empresa lhe fizesse as seguintes demandas:


a) Estimativa de vendas do produto para os próximos 6 meses;
b) Ele diz que o país está em um momento de crise e pode ser que a renda das famílias caiam, mas ele não tem certeza sobre essa afirmação. O que você propõe?;
c) Ele diz que a empresa está pensando em aumentar o preço do Chicabon em 10% e lhe pergunta se isso irá afetar as vendas;
d) O presidente lhe informa que o carnaval é a época do ano que mais vende Chicabon, mas ele está preocupado com as notícias de que o El Niño fará as temperaturas caírem aproximadamente 10% em todo Brasil e ele pergunta como serão as vendas se isso realmente acontecer;
e) Ele lembra também de uma campanha de marketing que ocorrerá no final do ano e que pretende fazer um contrato com a empresa de marketing baseado no aumento das vendas, o que você propõe?

Imagem meramente ilustrativa. Os números usados no exercício não refletem a venda deste produto

Imagem meramente ilustrativa. Os números usados no exercício não refletem a venda deste produto


Bem, vamos responde-lo seguindo as seguintes etapas:

Carregue o conjunto de dados e trace as variáveis cons (consumo de sorvete), temp (temperatura), renda e preço;

require(ggplot2)
require(gridExtra)

df <- read.csv("Icecream.csv")

p1 <- ggplot(df, aes(x = X, y = cons)) +
  ylab("consumo") +
  xlab("") +
  geom_line() +
  expand_limits(x = 0, y = 0)
p2 <- ggplot(df, aes(x = X, y = temp)) +
  ylab("Temperatura") +
  xlab("") +
  geom_line() +
  expand_limits(x = 0, y = 0)
p3 <- ggplot(df, aes(x = X, y = income)) +
  ylab("renda") +
  xlab("Period") +
  geom_line() +
  expand_limits(x = 0, y = 70)
p4 <- ggplot(df, aes(x = X, y = price)) +
  ylab("preço") +
  xlab("Period") +
  geom_line() +
  expand_limits(x = 0, y = 0.25)
grid.arrange(p1, p2, p3, p4, ncol=1, nrow=4)


Agora vamos estimar um modelo ARIMA para os dados sobre consumo de sorvete usando a função auto.arima. Em seguida, passe o modelo como entrada para a função forecast para obter uma previsão para os próximos 6 períodos (ambas as funções são do pacote forecast).

require(forecast)
fit_cons <- auto.arima(df$cons)
fcast_cons <- forecast(fit_cons, h = 6)
autoplot(fcast_cons)


Bem, parece que atendemos o primeiro pedido do C level, mas será que podemos melhorar esse modelo?! Outra dúvida é: quanto estamos errando?


Vamos usar a função accuracy do pacote forecast para encontrar o erro médio de escala absoluta (MASE) do modelo ARIMA ajustado.

require(forecast)
accuracy(fit_cons)
##                        ME       RMSE        MAE        MPE     MAPE
## Training set 0.0001020514 0.03525274 0.02692065 -0.9289035 7.203075
##                   MASE       ACF1
## Training set 0.8200619 -0.1002901


Objetivando melhorar nosso modelo vamos estimar um modelo ARIMA estendido para os dados de consumo com a variável de temperatura como um regressor adicional (usando a função auto.arima). Em seguida, faremos uma previsão para os próximos 6 períodos (note que esta previsão requer uma suposição sobre a temperatura esperada; suponha que a temperatura para os próximos 6 períodos será representada pelo seguinte vetor: fcast_temp <- c (70,5, 66, 60,5, 45,5, 36, 28)).

require(forecast)
fit_cons_temp <- auto.arima(df$cons, xreg = df$temp)
fcast_temp <- c(70.5, 66, 60.5, 45.5, 36, 28)
fcast_cons_temp <- forecast(fit_cons_temp, xreg = fcast_temp, h = 6)
autoplot(fcast_cons_temp)

summary(fcast_cons_temp)
## 
## Forecast method: Regression with ARIMA(0,1,0) errors
## 
## Model Information:
## Series: df$cons 
## Regression with ARIMA(0,1,0) errors 
## 
## Coefficients:
##         xreg
##       0.0028
## s.e.  0.0007
## 
## sigma^2 estimated as 0.001108:  log likelihood=58.03
## AIC=-112.06   AICc=-111.6   BIC=-109.32
## 
## Error measures:
##                       ME       RMSE        MAE      MPE     MAPE      MASE
## Training set 0.002563685 0.03216453 0.02414157 0.564013 6.478971 0.7354048
##                    ACF1
## Training set -0.1457977
## 
## Forecasts:
##    Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
## 31      0.5465774 0.5039101 0.5892446 0.4813234 0.6118313
## 32      0.5337735 0.4734329 0.5941142 0.4414905 0.6260566
## 33      0.5181244 0.4442225 0.5920263 0.4051012 0.6311476
## 34      0.4754450 0.3901105 0.5607796 0.3449371 0.6059529
## 35      0.4484147 0.3530078 0.5438217 0.3025024 0.5943270
## 36      0.4256524 0.3211393 0.5301654 0.2658135 0.5854913


Ao observar o MAPE e o MASE percebemos que apenas inserindo a variável temperatura no nosso modelo já conseguimos melhorar o nosso erro de previsão. Agora nós vamos ver se a variável temperatura é estatisticamente significante no nível de 5% para o nosso modelo;

require(lmtest)
coeftest(fit_cons_temp)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## xreg 0.0028453  0.0007302  3.8966 9.756e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Para tentar melhorar ainda mais nosso modelo vamos inserir as seguintes variáveis:

  1. Valores da variável de temperatura,
  2. Valores da variável renda com lag igual a 1 (estamos supondo que a queda da renda só irá afetar o consumo no mês seguinte),
  3. Valores da variável de preço.

Para mostrar o que eu fiz, estou imprimindo a matriz com as variáveis que entrarão no modelo.

price_column <- matrix(df$price, ncol = 1)
temp_column <- matrix(df$temp, ncol = 1)
income <- c(NA, df$income)
income_matrix <- embed(income, 2)
vars_matrix <- cbind(temp_column, income_matrix[,2],price_column)
print(vars_matrix)
##       [,1] [,2]  [,3]
##  [1,]   41   NA 0.270
##  [2,]   56   78 0.282
##  [3,]   63   79 0.277
##  [4,]   68   81 0.280
##  [5,]   69   80 0.272
##  [6,]   65   76 0.262
##  [7,]   61   78 0.275
##  [8,]   47   82 0.267
##  [9,]   32   79 0.265
## [10,]   24   76 0.277
## [11,]   28   79 0.282
## [12,]   26   82 0.270
## [13,]   32   85 0.272
## [14,]   40   86 0.287
## [15,]   55   83 0.277
## [16,]   63   84 0.287
## [17,]   72   82 0.280
## [18,]   72   80 0.277
## [19,]   67   78 0.277
## [20,]   60   84 0.277
## [21,]   44   86 0.292
## [22,]   40   85 0.287
## [23,]   32   87 0.277
## [24,]   27   94 0.285
## [25,]   28   92 0.282
## [26,]   33   95 0.265
## [27,]   41   96 0.265
## [28,]   52   94 0.265
## [29,]   64   96 0.268
## [30,]   71   91 0.260



Agora, estimaremos um modelo onde essas variáveis são levadas em consideração.

require(BETS)
fit_vars <- auto.arima(df$cons, xreg = vars_matrix,lambda = 0)
summary(fit_vars)
## Series: df$cons 
## Regression with ARIMA(1,0,0) errors 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1  intercept   xreg1   xreg2    xreg3
##       0.4191     -1.967  0.0097  0.0110  -1.7383
## s.e.  0.1942      0.599  0.0012  0.0033   1.7208
## 
## sigma^2 estimated as 0.005138:  log likelihood=37.34
## AIC=-62.67   AICc=-59.02   BIC=-54.27
## 
## Training set error measures:
##                        ME      RMSE        MAE        MPE     MAPE
## Training set 0.0002430829 0.0261767 0.02000142 -0.3454377 5.434105
##                   MASE        ACF1
## Training set 0.6092868 -0.09847473
BETS.t_test(fit_vars)
##                 Coeffs  Std.Errors        t Crit.Values Rej.H0
## ar1        0.419134114 0.194161630 2.158687     2.04523   TRUE
## intercept -1.967007334 0.598993997 3.283851     2.04523   TRUE
## xreg1      0.009732735 0.001165196 8.352875     2.04523   TRUE
## xreg2      0.010998589 0.003258491 3.375363     2.04523   TRUE
## xreg3     -1.738334523 1.720785749 1.010198     2.04523  FALSE

Observando os resultados do modelo, percebemos que apenas a variável preço não é significante, contudo, como o C level perguntou especificamente sobre essa variável decidimos mante-la. Uma pergunta que precisamos fazer é se o nosso novo modelo está errando menos, vamos ver isso agora?!

require(forecast)
accuracy(fit_vars)
##                        ME      RMSE        MAE        MPE     MAPE
## Training set 0.0002430829 0.0261767 0.02000142 -0.3454377 5.434105
##                   MASE        ACF1
## Training set 0.6092868 -0.09847473

Como podemos observar, tanto o MAPE quanto o MASE melhoraram significativamente. Agora, nós vamos usar este modelo para fazer uma previsão para os próximos 6 períodos. Observe que a previsão requer uma matriz da temperatura esperada, renda e preço para os próximos 6 períodos. Vamos criar uma matriz usando a variável fcast_temp e os seguintes valores para a renda esperada: 91, 91, 93, 96, 96, 96 e do preço esperado: 0.282, 0.265, 0.265, 0.265, 0.268, 0.267.

require(forecast)
expected_temp_income_price <- matrix(c(fcast_temp, 91, 91, 93, 96, 96, 96,0.282, 0.265, 0.265, 0.265, 0.268, 0.267),
                               ncol = 3, nrow = 6)
fcast_cons_temp_income_price <- forecast(fit_vars,
                                   xreg = expected_temp_income_price,
                                   h = 6)
autoplot(fcast_cons_temp_income_price)

summary(fcast_cons_temp_income_price)
## 
## Forecast method: Regression with ARIMA(1,0,0) errors
## 
## Model Information:
## Series: df$cons 
## Regression with ARIMA(1,0,0) errors 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1  intercept   xreg1   xreg2    xreg3
##       0.4191     -1.967  0.0097  0.0110  -1.7383
## s.e.  0.1942      0.599  0.0012  0.0033   1.7208
## 
## sigma^2 estimated as 0.005138:  log likelihood=37.34
## AIC=-62.67   AICc=-59.02   BIC=-54.27
## 
## Error measures:
##                        ME      RMSE        MAE        MPE     MAPE
## Training set 0.0002430829 0.0261767 0.02000142 -0.3454377 5.434105
##                   MASE        ACF1
## Training set 0.6092868 -0.09847473
## 
## Forecasts:
##    Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
## 31      0.4879484 0.4451220 0.5348952 0.4239943 0.5615492
## 32      0.4665596 0.4223279 0.5154238 0.4006367 0.5433298
## 33      0.4463201 0.4034828 0.4937054 0.3824964 0.5207935
## 34      0.3964974 0.3583608 0.4386924 0.3396806 0.4628175
## 35      0.3587917 0.3242689 0.3969898 0.3073594 0.4188304
## 36      0.3321792 0.3002150 0.3675467 0.2845588 0.3877689
accuracy(fcast_cons_temp_income_price)
##                        ME      RMSE        MAE        MPE     MAPE
## Training set 0.0002430829 0.0261767 0.02000142 -0.3454377 5.434105
##                   MASE        ACF1
## Training set 0.6092868 -0.09847473


Bem, agora que já temos o nosso modelo definido, vamos tentar responder todas as perguntas do C level.

a) Estimativa de vendas do produto para os próximos 6 meses:
O modelo “fcast_cons_temp_income_price” nos dá essa previsão. É importante observar o que estamos esperando das demais variáveis que entraram no modelo.


b) Ele diz que o país está em um momento de crise e pode ser que a renda das famílias caiam, mas ele não tem certeza sobre essa afirmação. O que você propõe?

Neste caso, nossa proposta é fazer uma previsão das vendas do Chicabon em que a renda caia 5%, neste caso o C level terá um cenário pessimista para as vendas, uma vez que o cenário “provável” já foi dado na resposta anterior. Vamos fazer isso:
* Observe que reduzimos a nossa previsão para a renda em 5%;

require(forecast)
expected_temp_income_5_price <- matrix(c(fcast_temp, 86.45, 86.45, 88.35, 91.2, 91.2, 91.2,0.282, 0.265, 0.265, 0.265, 0.268, 0.267),
                               ncol = 3, nrow = 6)
fcast_expected_temp_income_5_price <- forecast(fit_vars,
                                   xreg = expected_temp_income_5_price,
                                   h = 6)
autoplot(fcast_expected_temp_income_5_price)

summary(fcast_expected_temp_income_5_price)
## 
## Forecast method: Regression with ARIMA(1,0,0) errors
## 
## Model Information:
## Series: df$cons 
## Regression with ARIMA(1,0,0) errors 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1  intercept   xreg1   xreg2    xreg3
##       0.4191     -1.967  0.0097  0.0110  -1.7383
## s.e.  0.1942      0.599  0.0012  0.0033   1.7208
## 
## sigma^2 estimated as 0.005138:  log likelihood=37.34
## AIC=-62.67   AICc=-59.02   BIC=-54.27
## 
## Error measures:
##                        ME      RMSE        MAE        MPE     MAPE
## Training set 0.0002430829 0.0261767 0.02000142 -0.3454377 5.434105
##                   MASE        ACF1
## Training set 0.6092868 -0.09847473
## 
## Forecasts:
##    Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
## 31      0.4641306 0.4233947 0.5087859 0.4032983 0.5341388
## 32      0.4437859 0.4017132 0.4902649 0.3810808 0.5168088
## 33      0.4240677 0.3833661 0.4690904 0.3634260 0.4948280
## 34      0.3761079 0.3399325 0.4161331 0.3222130 0.4390177
## 35      0.3403412 0.3075937 0.3765751 0.2915538 0.3972926
## 36      0.3150973 0.2847768 0.3486460 0.2699256 0.3678283
accuracy(fcast_expected_temp_income_5_price)
##                        ME      RMSE        MAE        MPE     MAPE
## Training set 0.0002430829 0.0261767 0.02000142 -0.3454377 5.434105
##                   MASE        ACF1
## Training set 0.6092868 -0.09847473
vendas_Nrenda <- fcast_expected_temp_income_5_price$mean*100000
vendas_Arenda <- fcast_cons_temp_income_price$mean*100000
comp_renda <- cbind(vendas_Arenda,vendas_Nrenda)

ts.plot(comp_renda,col= 1:2,lty=1:2)
legend("topright", legend = c("Renda Antiga","Renda Nova"), col = 1:2, lty = 1:2,cex = 0.8)

barplot(comp_renda,horiz = T)

summary(comp_renda)
##  vendas_Arenda   vendas_Nrenda  
##  Min.   :33218   Min.   :31510  
##  1st Qu.:36822   1st Qu.:34928  
##  Median :42141   Median :40009  
##  Mean   :41472   Mean   :39392  
##  3rd Qu.:46150   3rd Qu.:43886  
##  Max.   :48795   Max.   :46413


c) Ele diz que a empresa está pensando em aumentar o preço do Chicabon em 10% e lhe pergunta se isso irá afetar as vendas

require(forecast)
expected_temp_income_Nprice <- matrix(c(fcast_temp, 91, 91, 93, 96, 96, 96,0.3102, 0.2915, 0.2915, 0.2915, 0.2948, 0.286),
                               ncol = 3, nrow = 6)
fcast_cons_temp_income_Nprice <- forecast(fit_vars,
                                   xreg = expected_temp_income_Nprice,
                                   h = 6)
autoplot(fcast_cons_temp_income_price)

summary(fcast_cons_temp_income_Nprice)
## 
## Forecast method: Regression with ARIMA(1,0,0) errors
## 
## Model Information:
## Series: df$cons 
## Regression with ARIMA(1,0,0) errors 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1  intercept   xreg1   xreg2    xreg3
##       0.4191     -1.967  0.0097  0.0110  -1.7383
## s.e.  0.1942      0.599  0.0012  0.0033   1.7208
## 
## sigma^2 estimated as 0.005138:  log likelihood=37.34
## AIC=-62.67   AICc=-59.02   BIC=-54.27
## 
## Error measures:
##                        ME      RMSE        MAE        MPE     MAPE
## Training set 0.0002430829 0.0261767 0.02000142 -0.3454377 5.434105
##                   MASE        ACF1
## Training set 0.6092868 -0.09847473
## 
## Forecasts:
##    Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
## 31      0.4646055 0.4238279 0.5093064 0.4037109 0.5346853
## 32      0.4455546 0.4033143 0.4922189 0.3825996 0.5188686
## 33      0.4262264 0.3853176 0.4714783 0.3652761 0.4973469
## 34      0.3786467 0.3422271 0.4189420 0.3243879 0.4419810
## 35      0.3424599 0.3095086 0.3789193 0.2933687 0.3997658
## 36      0.3213871 0.2904614 0.3556055 0.2753138 0.3751707
accuracy(fcast_cons_temp_income_Nprice)
##                        ME      RMSE        MAE        MPE     MAPE
## Training set 0.0002430829 0.0261767 0.02000142 -0.3454377 5.434105
##                   MASE        ACF1
## Training set 0.6092868 -0.09847473


Vamos ver a nova quantidade vendida

vendas_Nprice <- fcast_cons_temp_income_Nprice$mean*100000
vendas_price <- fcast_cons_temp_income_price$mean*100000
comp <- cbind(vendas_price,vendas_Nprice)

ts.plot(comp,col= 1:2,lty=1:2)
legend("topright", legend = c("Preço Antigo","Novo Preço"), col = 1:2, lty = 1:2,cex = 0.8)

barplot(comp,horiz = T)

summary(comp)
##   vendas_price   vendas_Nprice  
##  Min.   :33218   Min.   :32139  
##  1st Qu.:36822   1st Qu.:35151  
##  Median :42141   Median :40244  
##  Mean   :41472   Mean   :39648  
##  3rd Qu.:46150   3rd Qu.:44072  
##  Max.   :48795   Max.   :46461


d) O presidente lhe informa que o carnaval é a época do ano que mais vende Chicabon, mas ele está preocupado com as notícias de que o El Niño fará as temperaturas caírem aproximadamente 10% em todo Brasil e ele pergunta como serão as vendas se isso realmente acontecer


e) Ele lembra também de uma campanha de marketing que ocorrerá no final do ano e que pretende fazer um contrato com a empresa de marketing baseado no aumento das vendas, o que você propõe?


Trading Itau vs Bradesco

Neste exercício vamos aplicar a teoria da cointegração para definir uma estratégia de operação na bolsa de valores. A ideia é simples: se duas séries tem dependência de longo prazo e, no curto prazo uma série temporal se descola da outra, espera-se que, após algum tempo, haverá uma reversão para a outra série temporal.

Para alcançar nosso objetivo, as seguintes atividades serão desenvolvidas:

  1. Vamos trabalhar com as séries de ações de dois dos principais bancos privados brasileiros: Bradesco e Itaú. Antes de iniciar nossa estratégia de operação, vamos fazer uma análise exploratória dessas duas séries temporais;

  2. Criar uma estratégia de curto prazo para operar as ações;

  3. Long-Short através de Cointegração: um exemplo usando o pacote PairTrading.


Análise exploratória: Bradesco e Itaú

Primeiramente, vamos consultar o valor das ações desses dois bancos e observa-las graficamente.

#install.packages("quantmod")
require(quantmod)
## Banco Bradesco S.A.
getSymbols('BBDC4.SA',src='yahoo')
## [1] "BBDC4.SA"
chartSeries(BBDC4.SA)

## Itaúsa - Investimentos Itaú S.A.
getSymbols('ITUB4.SA',src='yahoo')
## [1] "ITUB4.SA"
chartSeries(ITUB4.SA)

Que tal agora observarmos algum período da série?

## últimos 4 meses
chartSeries(BBDC4.SA,subset = 'last 4 months')

chartSeries(ITUB4.SA,subset = 'last 4 months')

## podemos também observar os gráficos das duas ST ao mesmo tempo
head(as.xts(merge(BBDC4.SA,ITUB4.SA)))
##            BBDC4.SA.Open BBDC4.SA.High BBDC4.SA.Low BBDC4.SA.Close
## 2007-01-02            NA            NA           NA             NA
## 2007-01-03            NA            NA           NA             NA
## 2007-01-04            NA            NA           NA             NA
## 2007-01-05            NA            NA           NA             NA
## 2007-01-08            NA            NA           NA             NA
## 2007-01-09            NA            NA           NA             NA
##            BBDC4.SA.Volume BBDC4.SA.Adjusted ITUB4.SA.Open ITUB4.SA.High
## 2007-01-02              NA                NA        90.937        92.892
## 2007-01-03              NA                NA        91.750        93.293
## 2007-01-04              NA                NA        92.233        92.952
## 2007-01-05              NA                NA        91.821        92.516
## 2007-01-08              NA                NA        89.158        90.254
## 2007-01-09              NA                NA        89.877        90.784
##            ITUB4.SA.Low ITUB4.SA.Close ITUB4.SA.Volume ITUB4.SA.Adjusted
## 2007-01-02       90.714       78.28006         5921845          16.50530
## 2007-01-03       91.090       78.74992         7301247          16.60435
## 2007-01-04       91.467       77.99983         6424324          16.44620
## 2007-01-05       88.346       74.99985         8512351          15.81366
## 2007-01-08       88.688       76.62003         8683065          16.15527
## 2007-01-09       85.495       74.98011        10453820          15.80949
chartSeries(c(BBDC4.SA, ITUB4.SA))


Modelo básico de curto prazo para operar ações: bandas de Bollinger

Outra ferramenta interessante do quantmod são das bandas de Bollinger (https://en.wikipedia.org/wiki/Bollinger_Bands).

O conceito é simples, mas lógico. Procura usar uma ideia clássica e sólida. A ideia de que os preços não se afastam por muito tempo de uma média de valores, procurando retornar a essa média de tempos em tempos. Logo, a banda é construída projetando uma quantidade de desvios padrões acima e abaixo desta média.

Com isso, conseguimos observar situações de sobrevenda e de sobrecompra, conseguimos estabelecer suportes e resistencias, conseguimos estabelecer tendência de alta ou de baixa e ainda observar a expansão ou retração da volatilidade do ativo (http://lseducacao.com.br/como-operar-bandas-de-bollinger/).

Bollinger Bands consist of:

a. an N-period moving average (MA)

b. an upper band at K times an N-period standard deviation above the moving average (MA + Kσ)

c. a lower band at K times an N-period standard deviation below the moving average (MA − Kσ)


Typical values for N and K are 20 and 2, respectively.

Conforme observado, o default é 20,2, mas o quantmod me permite escolher as bandas. Foi isso que eu fiz para as ações do Bradesco.

chartSeries(BBDC4.SA, subset = 'last 12 months', theme="white",TA="addVo();addBBands(30,2);addCCI()") 

chartSeries(ITUB4.SA, subset = 'last 12 months', theme="white",TA="addVo();addBBands();addCCI()")


Long-Short através de Cointegração: um exemplo usando o pacote PairTrading

Estratégias de arbitragem estatística são baseadas em encontrar uma série temporal que possua a característica de estacionariedade ou reversão à média. Isto significa que é possível identificar situações em que a série divergiu de seu comportamento histórico e prever com alguma segurança que a série convergirá ou reverterá para um comportamento “médio”. O conceito de cointegração formaliza matematicamente este comportamento e permite a realização de testes estatísticos para detectar séries com este comportamento.

No contexto de operações com pares de ativos (pairs trading), a existência de uma relação de cointegração entre as séries de preços de dois ativos significa que pode ser possível realizar operações lucrativas de arbitragem. Por outro lado, se o par não for cointegrado, será impossível encontrar uma relação consistente para operar o par.

É necessário termos um teste para identificar quais pares de ações são cointegrados. Mesmo dentro do universos dos pares que são cointegrados, não há garantia de sucesso. É preciso que o par possua algumas características específicas para que uma estratégia de arbitragem seja consistentemente lucrativa:

a. Relação de cointegração estável ao longo do tempo
b. Reversão frequente do spread à média
c. Variabilidade razoavelmente grande nas divergências


Quais são os passos para operar com pares de ações?

Para operar com pares de ações, os seguintes passos devem ser seguidos:

a - selecionar duas ações que movem similarmente;

b - ideia básica: vamos vender as ações com preço elevado e comprar as ações com preço baixo;

c - Monitorar as diferenças entre as duas ações (isso no curto prazo);

d - No longo prazo eu posso usar a teoria da cointegração para fazer esse monitoramento;

e - Regra de Decisão:

\[ spread = log(Y_t) - (alpha + beta*log(X_t))\], onde \(Y_t\) e \(X_t\) são ações.

Se \(spread\) > muito alto: compra-se \(X_t\) e vende-se \(Y_t\)
Se \(spread\) < muito baixo: compra-se \(Y_t\) e vende-se \(X_t\) 


*Selecionar duas ações que movem similarmente*
BBDC4.SA_2013 <- BBDC4.SA['2013::'] 
ITUB4.SA_2013 <- ITUB4.SA['2013::'] 
pairs <- cbind(BBDC4.SA_2013$BBDC4.SA.High,ITUB4.SA_2013$ITUB4.SA.High)
ts.plot(pairs)


Estimando o spread

O pacote PairTrading não está mais disponível na versão mais nova do R, dessa forma, para trabalharmos com o pacote iremos utilizar o pacote que está disponível no Github. Utilizando este pacote vamos estimar o spread entra as ações do Itaú e do Bradesco. A primeira coisa que precisamos fazer é estimar a regressão abaixo:

\[ log(Y_t) = alpha + beta*log(X_t) + u_t\]

# install.packages("devtools")
# require(devtools)
# install_github("cran/PairTrading")
require(PairTrading)

## Estimando o spread
reg <- EstimateParameters(pairs, method = lm)


Estimando os parâmetros para o back-test

Para rodar o back-test nós precisamos estimar os parâmetros historicamente usando a função “EstimateParametersHistorically”. Esse função faz algo como uma “rolling regression” para estimar os parâmetros. Por isso ela é diferente da função “EstimateParameter”.


#estimate parameters for back test
params <- EstimateParametersHistorically(pairs, period = 180)

** Criando o sinal para a operação**

Agora, nós vamos criar o sinal para a operação usando o “spread”. A função “Simple” nos dá uma estratégia de operação bem simples: se o spread é maior (menor) que um valor específico, nós iremos comprar (vender).

Neste nosso caso, nós fixamos nosso “valor específico” em \(0.05\), então, baseado na nossa regressão:

\[ log(ITAU_t) = alpha + beta*log(BRAD_t) + u_t\]

onde o spread é igual ao erro (\(u_t = spread\))

\[ spread_t = log(ITAU_t) - (alpha + beta*log(BRAD_t))\]

Então, iremos comprar Bradesco, se o spread for maior que \(0.05\) e comprar Itaú, se o spread for menor que \(-0.05\). O gráfico abaixo mostra os pontos em que devemos comprar uma ou outra ação.

manchas azuis na margem inferior (spread maior que \(0.05\)) -> compro Itaú manchas azuis na margem superior (spread menor que \(-0.05\)) -> compro Bradesco

#create & plot trading signals
signal <- Simple(params$spread, 0.05)
barplot(signal,col="blue",space = 0, border = "blue",xaxt="n",yaxt="n",xlab="",ylab="")
par(new=TRUE)

plot(params$spread)


e - Back-test performance

A questão é a seguinte: fiquei rico?!!

#Performance of pair trading
return.pairtrading <- Return(pairs, lag(signal), lag(params$hedge.ratio))
plot(100 * cumprod(1 + return.pairtrading))


Para meus alunos de graduação

  1. A ementa do curso está no meu site

  2. Avaliação

    Provas (70%) + (2 Testes + 1 trabalho) (30%)}
    Nota final = 0.7 (0.5 P1 + 0.5 P2) + 0.3 (T1 + T2)

  3. Observações importantes:

  • A regra estabelecida acima NÃO muda ao longo do curso;
  • A menor nota entre os dois testes e o trabalho será eliminada;
  • Não será tratado por e-mail nenhum assunto de prova;
  • Não haverá “arredodamento” da nota antes da PS;
  • Não é possível trocar as notas das provas pelas notas dos testes e vice-versa;
  • A revisão só será permitida se a prova ou o teste estiverem a caneta;
  • A revisão deverá ser feita por escrito e, posteriormente será analisada pelo professor;
  • O prazo para revisão é de até 1 (uma) semana após a revisão feita em sala.


Horário das aulas: terças e quintas das 07:30 às 09:35 hs

Materiais de Apoio

  1. Introduction to R for Data Science: este é um curso introdutório e irá ajuda-lo a entender conceitos básicos de programação em R;

  2. Vídeos sobre o R (FGV/IBRE | NMEC): vídeos em português produzidos pelo nosso time da FGV que facilitarão o entendimento de conceitos básicos de programação em R;

  3. Khan Academy: ideal para aprender conceitos básicos de matemática e estatística;

  4. BETS (Brazilian Economic Time Series package) – este é o pacote R que estamos desenvolvendo, fiquem à vontade para instala-lo. Iremos utiliza-lo no nosso curso;

# Instalação do pacote pelo CRAN 
install.package("BETS")
                      
# Para obter sempre a versão em desenvolvimento e acompanhar o precessod de criação:
install.packages("devtools")
require(devtools)
install_github("pedrocostaferreira/BETS")
require(BETS)
  1. Forecasting using R: este é um excelente curso do Rob Hyndman que poderá lhe ajudar a compreender um pouco mais o mundo de séries temporais.

Prof. Dr. Pedro Costa Ferreira

Doutor em Engenharia Elétrica - (Decision Support Methods) e Mestre em Economia. Co-autor dos livros “Planejamento da Operação de Sistemas Hidrotérmicos no Brasil” e “Análise de Séries Temporais em R: um curso introdutório”. É o primeiro pesquisador da América Latina a ser recomendado pela empresa RStudio Inc. Atuou em projetos de Pesquisa e Desenvolvimento (P&D) no setor elétrico nas empresas Light S.A. (e.g. estudo de contingências judiciais), Cemig S.A, Duke Energy S.A, entre outras. Ministrou cursos de estatística e séries temporais na PUC-Rio e IBMEC e em empresas como o Operador Nacional do Setor Elétrico (ONS), Petrobras e CPFL S.A. Atualmente é professor de Econometria de Séries Temporais e Estatística, cientista chefe do Núcleo de Métodos Estatísticos e Computacionais (FGV|IBRE), coordenador dos cursos Economia Descomplicada (FGV|IDE) e Big Data e Data Science (FGV|IDE) e sócio-diretor da empresa Model Thinking Br (MTBr). É também revisor de importantes journals, como Energy Policy e Journal of Applied Statistics. Principais estudos são em modelos Econométricos, Incerteza Econômica, Preços, R software e Business Cycle.

contatos:
email: pedro.guilherme@fgv.br
website: pedrocostaferreira.github.io
GitHub: github.com/pedrocostaferreira
Linkedin: linkedin.com/pedro-costa-ferreira

Referências

Bowles, M. 2015. “Ensemble Packages in R. See: Http://Blog. Revolutionanalytics. Com/2014/04/Ensemble-Packages-in-R. Html.” Accessed.

Team, R Core. 2017. “R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. Available Online at: H Ttp.” Www. R-Project. Org.

Varian, Hal R. 2014. “Big Data: New Tricks for Econometrics.” The Journal of Economic Perspectives 28 (2). American Economic Association: 3–27.