SAUDAÇÕES!

Seja bem vindo à página do professor Pedro Albuquerque. Para saber mais sobre meu currículo, disciplinas ministradas e interesses de pesquisa, navegue no menu disponível no topo da página.
Mostrando postagens com marcador Finanças. Mostrar todas as postagens
Mostrando postagens com marcador Finanças. Mostrar todas as postagens

sexta-feira, 15 de setembro de 2017

CAViaR - Conditional Autoregressive Value at Risk by Regression Quantiles

O modelo CAViaR - Conditional Autoregressive Value at Risk by Regression Quantiles é um modelo utilizado para se mensurar o Value at Risk para períodos diários e, assim, monitorar o risco de maneira autogregressiva.

Engle e Manganelli (2004) sugerem quatro possíveis funções para a estimação do CAViaR:

  • Adaptive: $VaR_{t}=VaR_{t-1}+\beta_{1}\left\{\left[1+\exp\left(G\left[y_{t-1}-VaR_{t-1}\right]\right)^{-1}\right]-\tau\right\}$.
  • Symmetric absolute value: $VaR_{t}=\beta_{1}+\beta_{2}VaR_{t-1}+\beta_{3}|y_{t-1}|$.
  • Asymmetric slope: $VaR_{t}=\beta_{1}+\beta_{2}VaR_{t-1}+\beta_{3}(y_{t-1})^{+}+\beta_{4}(y_{t-1})^{-}$.
  • Indirect GARCH(1, 1): $VaR_{t}=\left(\beta_{1}+\beta_{2}VaR_{t-1}^{2}+\beta_{3}y_{t-1}^{2}\right)^{1/2}$

onde $y_{t}$ é o retorno do ativo no tempo $t$ com $t=1,\dots, T$. $G$ é uma constante positiva (Engle e Manganelli (2004) sugerem $G=10$), $\tau$ é o nível do VaR desejado (usualmente $\tau=0.01$ ou $\tau=0.05$) e $(x)^{+}=\max(x,0)$ e $(x)^{-}=-\min(x,0)$

Como as funções são recursvivas é necessário uma primeira estimativa para o $VaR_{t=1}$. Considerando os mesmos dados de Engle e Manganelli (2004) os quais trabalharam com os retornos dos preços das açoes da GM, IBM e SP500 para Abril 7, 1986 até Abril 7, 1999, podemos obter a primeira estimativa de $VaR_{t=1}$ como:

#Limpa os workspace
rm(list=ls())
#Lê os dados
stocks<-read.table("dataCAViaR.txt")
colnames(stocks)<-c("GM","IBM","SP500")
#Nível do VaR desejado
tau<-0.05
#Divide a base em treinamento e validação
train<-stocks[1:2892,]
valid<-stocks[2893:3392,]
#Criar o vetor de VaR
VaR<-rep(NA,nrow(train))
#Primeira estimativa
VaR[1]<- -qnorm(tau)*sd(train$GM)
Em seguida, Engle e Manganelli (2004) sugerem estimar os parâmetros dos modelos por meio da estrutura de Regressão Quantílica: $\min_{\boldsymbol\beta}\frac{1}{T}\sum_{t=1}^{T}\left[\tau-1\left(y_{t} \leq VaR_{t}\right)\right]\left[y_{t}-VaR_{t}\right]$ onde $1(\cdot)$ é uma função indicadora que assume valor igual a 1 quando seu argumento é verdadeiro $y_{t}\lt VaR_{t}$ e 0 caso contrário. Assim, a função objetivo para o modelo Adaptive é dada por:
#Sugestão Engle e Manganelli (2004)
G<-10
#Function
adpatative<-function(x){
  beta1<-x[1]
  for(i in 2:nrow(train)){
    VaR[i]<- VaR[i-1] + beta1*((1+(exp(G*(train[i-1,"GM"] - VaR[i-1]))^(-1)))-tau)
  }
  #Objective Function
  res<-sum((tau-(train$GM < VaR))*(train[,"GM"]-VaR))/nrow(train)
  return(res)
}
O próximo passo é otimizar essa função, isso pode ser feito utilizando-se a função optim:
#Realiza a otimização usando o métood de Brent
res<-optim(par=c(0), adpatative, method="Brent", lower=c(-10), upper=c(10))
#Resultado do parâmetro de interesse
beta<-res$par
A função optim tem como argumentos nesse exemplo: valor inicial para o(s) parâmetro(s) de interesse, função a ser minimizada, método utilizado e limites inferiores e superiores para o(s) parâmetro(s) de interesse. O último passo é a previsão do VaR:
VaR<-rep(NA,nrow(dados))
VaR[1]<- -qnorm(tau)*sd(dados[,"GM"])
#Previsão do CAViaR
adpatativeForecast<-function(x){
  beta1<-x[1]
  for(i in 2:nrow(dados)){
    VaR[i]<- VaR[i-1] + beta1*((1+(exp(G*(dados[i-1,"GM"]-VaR[i-1]))^(-1)))-tau)
  }
  #Objective Function
  return(VaR)
}
onde o objeto dados pode ser as bases de treinamento e validação:
#Base de treinamento
dados<-train
#Faz a previsão
VaR.train<-adpatativeForecast(beta)
#Medida de acuracia hits
mean(dados$GM< -VaR.train)-tau

O mesmo processo pode ser realizado para as demais funções.

quarta-feira, 15 de junho de 2016

Dados financeiros para o R.


Existe uma série de dados livres em finanças e economia que estão disponíveis no R. As principais fontes de dados são listadas a a seguir:

Fonte Pacote Acesso livre Disponível CRAN Nome/Página
Yahoo, FRED, Oanda, Google Quantmod Sim Sim http://www.quantmod.com/
Quandl Quandl Sim Sim http://www.quandl.com/help/packages/r
TrueFX TFX Sim Sim http://rpubs.com/gsee/TFX
Bloomberg Rbbg Não Não http://findata.org/rbloomberg/
Interactive Broker IBrokers Não Sim https://www.interactivebrokers.com/en/main.php
Datastream rdatastream Não Não https://github.com/fcocquemas/rdatastream
Penn World Table pwt Sim Sim https://pwt.sas.upenn.edu/
Yahoo, FRED, Oanda fImport Sim Sim http://www.rmetrics.org/
ThinkNum Thinknum Sim Sim http://thinknum.com/

A disponibilidade desses dados possibilita a execução de modelos econométricos em finanças e economia, visando o ensino como também a pesquisa. Outras publicações que valem a pena conferir, são:

1) Quandl Package – 5,000,000 free datasets at the tip of your fingers!
2) Financial Data Accessible from R - Parte 1.
2) Financial Data Accessible from R - Parte 2.
2) Financial Data Accessible from R - Parte 3.
2) Financial Data Accessible from R - Parte 4.

terça-feira, 15 de setembro de 2015

Credit Score usando o R.


Credit Scoring é definido como sendo um modelo estatístico/econométrico o qual atribui uma medida de risco aos clientes de uma instituição financeira ou aos futuros clientes.

Usualmente, a avaliação dos clientes é realizada por meio de um Credit Scorecard o qual é um modelo estatístico para avaliação do risco estruturado de maneira a facilitar a tomada de decisão quanto a liberação do crédito ou não.

O uso de credit scorecard é muito popular principalmente para as organizações que lidam empréstimos como bancos.

Dentre as vantagens da utilização de um credit scorecard podemos listar:

  1. Credit Scorecard é implementado facilmente e pode ser monitorado ao longo do tempo.
  2. Pessoas sem o conhecimento técnico em estatística ou econometria podem utilizar facilmente o credit scorecard para tomar decisões.

As principais questões em Credit Score são:

  1. Quem receberá o crédito ?
  2. Quanto deverá ser esse crédito ?
  3. Quais as estratégias para a distribuição e cobrança do crédito ?

Em outras palavras, Credit Score é um conjunto de modelos de decisão e técnicas estatísticas subjacentes que auxiliam os credores na tomada de decisão quanto a concessão de crédito ao consumidor.

Estas técnicas informam quem receberá o crédito, quanto crédito deverá ser fornecido e quais estratégias operacionais melhorarão a rentabilidade dos devedores para os credores (Thomas, Eldelman e Crook, 2002).

Credit Score usando o R.

Para demonstrar como o Credit Score pode ser formulado usando o R, iremos trabalhar com os dados German.csv o qual representa um conjunto de dados de crédito para uma instituição financeira Alemã.

Os dados são compostos por 300 empréstimos "ruins" (por exemplo, ausência de pagamento ou atraso) e 700 empréstimos "bons" (por exemplo, pagamentos sem atraso). O objetivo é fornecer insumos para a tomada de decisão quanto aos futuros empréstimos com base nos padrões anteriormente observados.

É comum em Credit Score classificar as contas "ruins" como aquelas contas que em algum período de tempo apresentaram inadimplência por 60 dias ou mais (em empréstimos hipotecários, 90 dias ou mais é por vezes utilizado).

Vamos importar os dados para o R:

#Limpa o Workspace
rm(list=ls())

#Importa os dados German.csv 
dados.df<-read.csv("https://dl.dropboxusercontent.com/u/36068691/Blog/german.csv")

#Apresenta as variáveis do DataFrame
names(dados.df)

#Apresenta a estrutura do DataFrame
str(dados.df)
No R as variáveis categóricas ou binárias são usualmente tratadas como fatores enquanto as variáveis contínuas ou discretas são tratadas como valores numéricos. Nesse caso, algumas conversões são necessárias:
#Transforma em fatores as variáveis categóricas e "dummies"
dados.df[,"CHK_ACCT"]     <-as.factor(dados.df[,"CHK_ACCT"])
dados.df[,"HISTORY"]      <-as.factor(dados.df[,"HISTORY"])
dados.df[,"NEW_CAR"]      <-as.factor(dados.df[,"NEW_CAR"])
dados.df[,"USED_CAR"]     <-as.factor(dados.df[,"USED_CAR"])
dados.df[,"FURNITURE"]    <-as.factor(dados.df[,"FURNITURE"])
dados.df[,"RADIO.TV"]     <-as.factor(dados.df[,"RADIO.TV"])
dados.df[,"EDUCATION"]    <-as.factor(dados.df[,"EDUCATION"])
dados.df[,"RETRAINING"]   <-as.factor(dados.df[,"RETRAINING"])
dados.df[,"SAV_ACCT"]     <-as.factor(dados.df[,"SAV_ACCT"])
dados.df[,"EMPLOYMENT"]   <-as.factor(dados.df[,"EMPLOYMENT"])
dados.df[,"MALE_DIV"]     <-as.factor(dados.df[,"MALE_DIV"])
dados.df[,"MALE_SINGLE"]  <-as.factor(dados.df[,"MALE_SINGLE"])
dados.df[,"MALE_MAR"]     <-as.factor(dados.df[,"MALE_MAR"])
dados.df[,"CO.APPLICANT"] <-as.factor(dados.df[,"CO.APPLICANT"])
dados.df[,"GUARANTOR"]    <-as.factor(dados.df[,"GUARANTOR"])
dados.df[,"TIME_RES"]     <-as.factor(dados.df[,"TIME_RES"])
dados.df[,"REAL_ESTATE"]  <-as.factor(dados.df[,"REAL_ESTATE"])
dados.df[,"PROP_NONE"]    <-as.factor(dados.df[,"PROP_NONE"])
dados.df[,"OTHER_INSTALL"]<-as.factor(dados.df[,"OTHER_INSTALL"])
dados.df[,"RENT"]         <-as.factor(dados.df[,"RENT"])
dados.df[,"OWN_RES"]      <-as.factor(dados.df[,"OWN_RES"])
dados.df[,"NUM_CREDITS"]  <-as.factor(dados.df[,"NUM_CREDITS"])
dados.df[,"JOB"]          <-as.factor(dados.df[,"JOB"])
dados.df[,"NUM_DEPEND"]   <-as.factor(dados.df[,"NUM_DEPEND"])
dados.df[,"TELEPHONE"]    <-as.factor(dados.df[,"TELEPHONE"])
dados.df[,"FOREIGN"]      <-as.factor(dados.df[,"FOREIGN"])

#Variável dependente
dados.df[,"RESPONSE"]     <-as.factor(dados.df[,"RESPONSE"])

#Transforma em numeric
dados.df[,"AMOUNT"]       <-as.numeric(dados.df[,"AMOUNT"])
dados.df[,"INSTALL_RATE"] <-as.numeric(dados.df[,"INSTALL_RATE"])
dados.df[,"AGE"]          <-as.numeric(dados.df[,"AGE"])
dados.df[,"NUM_DEPEND"]   <-as.numeric(dados.df[,"NUM_DEPEND"])
dados.df[,"DURATION"]     <-as.numeric(dados.df[,"DURATION"])
O próximo passo é separar os dados em dois grupos:
  1. Dados para estimação. (Treinamento)
  2. Dados para teste. (Validação)
Uma sugestão é separar a base de dados (aleatoriamente) da seguinte forma: 60% das observações deverão compor a base de treinamento e 40% a base de validação:
#Índices obtidos após a aleatorização
ordena <- sort(sample(nrow(dados.df), nrow(dados.df)*.6))

#Dados para o treinamento
treinamento<-dados.df[ordena,]

#Dados para a validação
validacao<-dados.df[-ordena,]
A ideia é construir o(s) modelo(s) de Credit Scoring com o DataFrame "treinamento" e em seguida avaliar o ajuste com o DataFrame "validacao". Uma das formas mais simples para modelar os dados de crédito é por meio da regressão logística:
#Índices obtidos após a aleatorização
ordena = sort(sample(nrow(dados.df), nrow(dados.df)*.6))

#Dados para o treinamento
treinamento<-dados.df[ordena,]

#Dados para a validação
validacao<-dados.df[-ordena,]
Como há muitas possíveis variáveis para o modelo, podemos proceder com a abordagem de Stepwise para selecionar o modelo com a "melhor" combinação de variáveis explicativas:
#Regressão Logística
modelo.completo <- glm(RESPONSE ~ . ,family=binomial,data=treinamento)

#Abordagem Stepwise para seleção de variáveis
stepwise <- step(modelo.completo,direction="both") 
Após algumas iterações, observa-se que o conjunto de variáveis com o menor valor para o Critério de Informação de Akaike é:
#Modelo com as variáveis indicadas pelo Stepwise
stepwise <- glm(RESPONSE ~  JOB+NUM_CREDITS+EMPLOYMENT+RETRAINING+NEW_CAR+TELEPHONE+MALE_DIV+
                  FURNITURE+PROP_NONE+MALE_MAR+RENT+NUM_DEPEND+REAL_ESTATE+EDUCATION+FOREIGN+
                  TIME_RES, family=binomial,data=treinamento)

#Resume os resultados do modelo
summary(stepwise)
Percebemos que nem todas as variáveis são significantes:
Uma medida interessante para interpretar o modelo é a medida de Razão de chances (Odds Ratio):
#Calcula a razão de chances
exp(cbind(OR = coef(stepwise), confint(stepwise)))
Obtemos como resultados as seguintes Razões de chance:
E como dito anteriormente, a interpretação é bem interessante. Veja por exemplo a variável NUM_DEPEND, nesse caso, para cada dependente a mais que um proponente possui isso aumenta a sua chance de ser considerado inadimplente em aproximadamente 12%. Finalmente, vamos testar a qualidade do modelo aplicando o modelo estimado na base de validação para termos uma ideia do grau de acerto desse modelo:
#Faz a previsão para a base de validação (probabilidade)
predito<-predict(stepwise,validacao,type="response")

#Escolhe quem vai ser "1" e quem vai ser "0"
predito<-ifelse(predito>=0.8,1,0)
  
#Compara os resultados
table(predito,validacao$RESPONSE)
Obtemos assim a seguinte matriz de confusão:

Logo a nossa taxa de acerto (acurácia) nesse modelo é dada por:

$
Tx.Acerto=\frac{75+166}{400}=\frac{241}{400} \approx 60\%
$

Podemos melhorar a taxa de acerto refinando o modelo por meio da exclusão de variáveis não significantes e pela inclusão de componentes estatisticamente significantes.

quarta-feira, 15 de abril de 2015

Text Mining com kernels string no R - Parte 3.


Nesse post daremos continuidade ao que já vimos em posts anteriores como Text Mining com kernels string no R - Parte 1 e Text Mining com kernels string no R - Parte 2. Nessa etapa, vamos tentar reconhecer o padrão das informações prestadas pela Petrobras de 2006 até o primeiro trimestre 2014, para isso, considere os seguintes pacotes:

#Habilita as bibliotecas necessárias
library(tm)
library(wordcloud)
library(Rstem)

Vamos trabalhar com os textos informados pela Petrobras, os quais podem ser obtidos nesse link. Salve esses arquivos em alguma pasta do seu computador, como por exemplo, C:\Text Mining. Em seguida, precisamos transformar os arquivos em PDF para o formato TXT, o qual é possível de ser lido pelo R:

#Seleciona o arquivo de texto a ser minerado
dest<- "C:\\Text Mining\\2006 - 1T.pdf"

#Executa o programa que converte pdf em txt
exe<-"C:\\Program Files\\xpdfbin-win-3.03\\bin64\\pdftotext.exe"
system(paste("\"", exe, "\" \"", dest,"\"", sep= ""), wait= F)

#Cria o arquivo no formato txt
filetxt<- sub(".pdf", ".txt", dest)

#Lê os dados do arquivo
txt<- readLines(filetxt, warn=FALSE)

#Deixa tudo minúsculo
txt<- tolower(txt)

#Remove algumas expressões regulares
txt<-removeWords(txt, c("\\f", stopwords("portuguese")))

#Cria os objetos no formato corpus
corpus<-Corpus(VectorSource(txt))

#Remove a pontuação
corpus<-tm_map(corpus,removePunctuation)

#Remove os números
corpus<-tm_map(corpus,removeNumbers)
Note que o código acima só é possível de ser executado se você tiver instalado o arquivo que converte PDF para TXT, qual seja, Xpdf. Na linha número 5 deve estar o endereço do arquivo executável pdftotext.exe o qual é obtido ao se instalar o programa Xpdf. Cada texto em PDF no arquivo é transformado em TXT e pode ser lido. Para cada texto, podemos associar um retorno obtido no trimestre após a divulgação da informação, por exemplo, para o primeiro trimestre de 2006 temos:
#Coloca os textos na forma de Matriz
tdm<- TermDocumentMatrix(corpus)
m<- as.matrix(tdm)

#Cria um data.frame com a frequência das palavras
d<- data.frame(freq= sort(rowSums(m), decreasing =TRUE))

#Coloca as palavras nas colunas
d$word <- row.names(d)
d$stem<- wordStem(d$word, language="portuguese")

#Remove as palavras muito grandes
d<- d[nchar(row.names(d))<20,]

#Agrega as frequências das palavras
agg_freq<- aggregate(freq ~stem, data = d, sum)
agg_word<- aggregate(word ~stem, data =d, function(x) x[1])

#Une os dados de frequência e palavras
d<- cbind(freq = agg_freq[,2], agg_word)

#Ordena as palavras pela frequência
d<- d[order(d$freq,decreasing =TRUE),]

#1 Trimestre 2006
d$Trimestre<-"1t2006"
d$Retorno<-0.03789927 #Retorno do 2t2006
Uma vez criado um DataFrame para cada texto com o seu respectivo retorno trimestral, podemos usar as frequências das palavras como uma covariável na previsão do retorno trimestral. Essa proposta não é nova e tem sido usada em finanças na busca por padrões. Um bom texto para se começar a estudar esse assunto é o artigo Text mining for market prediction: A systematic review.

sábado, 5 de janeiro de 2013

Redes neurais no R: Aplicações em finanças.


Uma rede neural artificial, muitas vezes chamada apenas de rede neural, é um modelo matemático inspirado pelas redes neurais biológicas.

Uma rede neural consiste em um grupo de neurônios artificiais interligados e processa a informação através de uma abordagem conexionista à computação. Na maioria dos casos, uma rede neural é um sistema adaptável, que muda a sua estrutura durante uma fase de aprendizagem.

As redes neurais são usadas para modelar relações complexas entre entradas e saídas ou para encontrar padrões em dados.

Considere, por exemplo, a órbita dos planetas em torno do sol ou o calendário das mares.

A ideia central do aprendizado de máquina é:

"Pode-se usar o computador para descobrir e descrever padrões baseados em dados ?"

De maneira formal, a hipótese básica do aprendizado de máquina é a Hipótese de Aprendizagem Indutiva:

Uma função adequada encontrada para modelar a função alvo para um conjunto suficientemente grande de dados também irá funcionar para modelar adequadamente exemplos não observados.

Isso quer dizer que se encontrarmos alguma "boa fórmula" para o movimento dos planetas em torno do sol, por exemplo, e que essa fórmula tenha sido construída com base em uma amostra suficientemente grande, espera-se que essa "boa fórmula" funcione bem "out-of-sample".

Nesse post mostrarei como usar o pacote neuralnet com uma aplicação a finanças.

Pacote neuralnet.

O pacote neuralnet foi construído de forma a ser possível trabalhar com Perceptron Multi-camadas (multi-layer perceptrons) no contexto da análise de regressão, isto é, para aproximar as relações funcionais entre variáveis independentes e variáveis resposta.

Assim, as redes neurais podem ser utilizadas como extensões do modelo linear generalizado.

Uma das vantagens do pacote neuralnet é a possibilidade de poder se trabalhar com um número arbitrário de covariáveis e​​ também de variáveis ​​de resposta, assim como o número de camadas ocultas.

Redes neurais.

Em muitas situações, a relação funcional entre as covariáveis ​​(também conhecidas como variáveis ​​de entrada ou variáveis independentes) e as variáveis ​​de resposta (também conhecidas como variáveis ​​de saída ou variáveis dependentes) é de grande interesse.

As redes neurais artificiais podem ser aplicadas como aproximação a qualquer relação funcional complexa.

Ao contrário dos modelos lineares generalizados (McCullagh e Nelder, 1989), não é necessário que o tipo de relação entre variáveis ​​dependentes e variáveis ​​resposta seja, por exemplo, uma combinação linear.

Isso faz das redes neurais artificiais uma valiosa ferramenta quantitativa. Esses modelos são, particularmente, extensões diretas dos modelos lineares generalizados e podem ser aplicados de maneira semelhante.

Dados observados são usados ​​para treinar a rede neural, fazendo assim com que a rede neural "aprenda" uma aproximação da relação entre as variáveis independentes e dependentes de forma iterativa por meio da adaptação contínua dos seus parâmetros. De fato, usando a nomenclatura estatística, os parâmetros do modelo são estimados por meio da amostra.

Perceptron Multi-Camadas.

O pacote neuralnet se concentra nos modelos Perceptron Multi-Camadas (Multi-Layer Perceptrons (MLP))(Bishop, 1995), os quais são úteis na modelagem por meio de relações funcionais entre as variáveis.

A estrutura subjacente de um MLP é um grafo orientado, isto é, consiste de vértices (neurônios) e arestas (sinapses).

Os neurônios são organizados em camadas, que são normalmente ligadas por sinapses. No pacote neuralnet, uma sinapse só pode se conectar a camadas posteriores.

A camada de entrada é constituída por todas as covariáveis (um neurônio para cada covariável) separadas por neurônios (camadas ocultas) até as variáveis ​​resposta.

Essas camadas intermédias são denominadas camadas ocultas (ou variáveis latentes), por não serem diretamente observáveis​​.

As camadas de entrada e as camadas ocultas incluem um neurônio constante, o qual estará associado ao intercepto em modelos lineares, ou seja, não é diretamente influenciado por qualquer covariável.


A figura anterior retirada de Günther e Fritsch (2010) representa um exemplo de uma rede neural com dois neurônios de entrada (A e B) e um neurônio de saída (Y), além de uma camada oculta composta de três neurônios.

Um peso (parâmetro) está associado a cada uma das sinapses, representando o efeito correspondente do neurônio e de todos os dados passam pela rede neural como sinais.

Os sinais são processados ​​inicialmente pela função de integração combinando todos os sinais de entrada e, em seguida, pela função de ativação transformando os resultados do neurônio.

O modelo perceptron multicamada mais simples consiste de uma única camada de entrada com $n$ covariáveis ​e uma camada de saída com um único neurônio de saída. Esse modelo calcula a seguinte função:


onde $w_{0}$ denota o intercepto, $\mathbf{w} = (w_{1},\dots, w_{n})$ o vetor de todos os demais pesos (parâmetros), e $\mathbf{x} = (x_{1},\dots, x_{n})$ o vetor de todas as covariáveis.

A função é matematicamente equivalente à estrutura padrão do modelo linear generalizado com função de ligação $f^{-1}(.)$. Portanto, todos os pesos calculados são, neste caso, equivalentes aos parâmetros da regressão via MLG.

Para aumentar a flexibilidade da modelagem mais camadas ocultas podem ser incluídas, aumentando assim a "não-linearidade" do modelo. No entanto, Hornik et al. (1989) mostraram que uma única camada oculta é suficiente para modelar qualquer função contínua por partes.

O modelo perceptron multicamada com uma camada oculta consistindo de $J$ neurônios calcula a seguinte função:


onde $w_{0}$ denota o intercepto do neurônio de saída e $w_{0j}$ representa o intercepto do $j$-ésimo neurônio oculto. Além disso, $w_{j}$ denota o peso sináptico correspondente à sinapse começando no $j$-ésimo neurônio oculto e que conduz ao neurônio de saída. $\mathbf{w_{j}} = (w_{1j},\dots, w_{nj})$ o vetor de todos os pesos sinápticos correspondentes às sinapses que conduzem ao $j$-ésimo neurônio oculto, e $\mathbf{x} = (x_{1},\dots, x_{n})$ é o vetor de todas as covariáveis.

Apesar das redes neurais serem extensões diretas dos MLG, os parâmetros não podem ser interpretados da mesma maneira.

De maneira formal, todos os neurônios ocultos e os neurônios de saída calculam a seguinte função: $f(g(z_{0},z_{1},\dots, z_{k})) = f(g(\mathbf{z}))$ a partir das saídas de todos os neurônios precedentes $z_{0}, z_{1},\dots, z_{k}$, onde $g:\mathbb{R}^{k+1}\rightarrow \mathbb{R}$ representa a função de integração e $f:\mathbb{R}\rightarrow \mathbb{R}$ é a função de ativação. O neurônio unitário $z_{0}$ é uma constante o qual está relacionado com o conceito de intercepto em modelos de regressão.

A função de integração é, muitas vezes, definida como $g(\mathbf{z})= w_{0}z_{0} + \sum_{i=1}^{n}w_{i}z_{i}= w0 + \mathbf{w}^{'}\mathbf{z}$. A função de ativação $f$ é geralmente uma função não-linear, não-decrescente, limitada e também diferenciável tal como o função logística $f(u) = 1/(1+\exp^{-u})$ ou a tangente hiperbólica.

Essa função deve ser escolhida em relação à variável de resposta, assim como nos modelos lineares generalizados. A função logística é, por exemplo, apropriada para variáveis ​​resposta binárias, uma vez que mapeia a saída de cada neurônio para o intervalo $[0,1]$. O pacote neuralnet usa a mesma função de integração, bem como função de ativação para todos os neurônios.

Aprendizado supervisionado.


As redes neurais são estimadas por meio de um processo de treinamento da rede que usa os dados para "aprender". Especificamente, o pacote neuralnet concentra-se em algoritmos de aprendizado supervisionado.

Estes algoritmos de aprendizagem são caracterizados pela utilização de saídas (outputs, resultados ou ainda variáveis dependentes), as quais são comparadas com o valor predito pela rede neural, adaptando dinamicamente os valores dos parâmetros de modo que o "erro" seja minimizado.

Os parâmetros de uma rede neural são os seus pesos. Todos os pesos são geralmente iniciados com valores aleatórios provenientes de uma distribuição normal padrão. Durante um processo iterativo de formação da rede, as seguintes etapas são repetidas:


  • A rede neural calcula uma saída de $o(\mathbf{x})$ para as entradas dadas $\mathbf{x}$ e para os parâmetros correntes (pesos atuais). Se o processo de formação ainda não estiver concluído, os resultados previstos serão diferentes da saída $\mathbf{y}$ observada.
  • Uma função de erro, como a Soma dos Quadrados dos Erros (SSE - Sum of Squared Errors):
    ou a entropia cruzada:

    mede a diferença entre a saída prevista e o resultado observado, tal que $l=1,\dots,L$ representam os índices para as observações, ou seja, é o par ordenado com os dados de entrada e saída, e $H=1,\dots,H$ representam os nós de saída.
  • Todos os pesos são adaptados segundo alguma regra de aprendizagem definida a priori.

O processo termina se um critério pré-determinado é atingido, por exemplo, se todas as derivadas parciais da função de erro com respeito aos pesos $(\partial E/\partial \mathbf{w})$ são menores do que um dado limiar. Um algoritmo de aprendizagem amplamente utilizado é o algoritmo resilient backpropagation.

Resilient backpropagation.


O algoritmo resilient backpropagation é baseado no algoritmo de retropropagação tradicional, o qual modifica os pesos de uma rede neural, a fim de encontrar um mínimo local para a função de erro estipulada.

Em outras palavras, o gradiente da função de erro $(\partial E /\partial\mathbf{w})$ é calculado em relação aos pesos, a fim de encontrar uma raiz. Particularmente, os pesos são modificadas para irem na direção oposta das derivadas parciais até que um mínimo local seja atingido.

Esta ideia básica é aproximadamente ilustrada na figura abaixo retirada do texto de Günther e Fritsch (2010), para uma função de erro univariada:


Caso a derivada parcial seja negativa, o peso é aumentado (lado esquerdo da figura) e se a derivada parcial for positiva, o peso é reduzido (parte da direita da figura). Isto garante que um mínimo local para a função de erro seja atingido.

Todas as derivadas parciais são calculadas usando a regra da cadeia, desde que a função do cálculo de uma rede neural seja basicamente uma composição de funções de integração e ativação. Uma explicação detalhada é dada em Rojas (1996).

O pacote neuralnet possibilita a escolha do método entre o clássico backpropagation e o resilient backpropagation, com retrocesso do peso (Riedmiller, 1994) ou sem retrocesso de peso (Riedmiller e Braun, 1993) e também a versão modificada globalmente proposta por Anastasiadis et al. (2005).

Aplicação em finanças.

A aplicação de redes neurais em finanças não é nova. Wonga e Selvib (1998) apresentam uma ampla revisão da utilização de redes neurais em finanças entre os anos 1990 a 1996. Os autores investigaram a tendência de aplicações produzidas no período 1990-1996. Wonga e Selvib (1998) analisaram a literatura de acordo com os seguintes critérios: (1) ano de publicação, (2) área de aplicação, (3) domínio do problema, (4) fase do processo de decisão, (5) nível de gestão (6), nível de interdependência de tarefas, (7) meios de desenvolvimento, (8) de interação corporativa / acadêmica, (9) Tecnologia / técnica estatística, e (10) estudo comparativo.

Yoon e Swales (1991) apresentam uma discussão sobre a previsibilidade dos preços de ações por meio de redes neurais. Segundo os autores, a previsão do desempenho do preço de ações é um problema difícil e complexo. Técnicas multivariadas quantitativas e qualitativas têm sido repetidamente usadas para auxiliar na formação da expectativas dos investidores quanto aos preços de ações. Yoon e Swales (1991) analisaram a capacidade das redes neurais e seu poder de previsão em contraste com o método de análise discriminante.

Nesse post farei uma breve aplicação do método de redes neurais com o intuito de avaliar o poder de previsibilidade das ações por meio dessa ferramenta.

Inicialmente, vamos obter a série temporal (01/01/2001 a 31/12/2012) para o ativo PETR4.SA por meio do pacote quantmod:

#Limpa o Workspace
rm(list=ls())

#Habilita o pacote quantmod
library(quantmod)

#Início do período de interesse
inicio = as.Date("2001-01-01") 

#Fim do período de interesse
fim = as.Date("2012-12-31") 

#Obtêm os dados da PETR4
getSymbols("PETR4.SA", src="yahoo",from=inicio,to=fim)

Nesse caso, queremos tentar construir um preditor na forma de um Modelo Autoregressivo de ordem 5, isto é:

$y_{t}=\phi_{0}+\phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\phi_{3}y_{t-3}+\phi_{4}y_{t-4}+\phi_{5}y_{t-5}$

Nesse caso, o valor da ação no tempo $t$ seria predita pelos $5$ valores (diários) que a antecedem:

#Dados para o neuralnet
t0<-as.numeric(Cl(PETR4.SA))            #Cinco dias antes
t0<-t0[-((length(t0)-4):length(t0))]

t1<-as.numeric(Cl(PETR4.SA)) [-1]       #Quatro dias antes
t1<-t1[-((length(t1)-3):length(t1))]

t2<-as.numeric(Cl(PETR4.SA)) [-c(1,2)]  #Três dias antes
t2<-t2[-((length(t2)-2):length(t2))]

t3<-as.numeric(Cl(PETR4.SA)) [-(1:3)]   #Dois dias antes
t3<-t3[-((length(t3)-1):length(t3))]

t4<-as.numeric(Cl(PETR4.SA)) [-(1:4)]   #Um dia antes
t4<-t4[-length(t0)]

t5<-as.numeric(Cl(PETR4.SA)) [-(1:5)]   #Variável dependente
Em seguida vamos dividir a série temporal em duas partes:
  1. Período de estimação (01/01/2001 até 31/12/2006).
  2. Período de validação (01/01/2007 até 31/12/2012).
#Cria a base
dados<-cbind(t1,t2,t3,t4,t5)

#Dados para estimação:
dados.Treina<-dados[1:1457,]

#Dados para validação
dados.Valida<-dados[1458:2917,]
Usando os dados para estimação e trabalhando com um modelo de redes neurais com 1 camada com 7 neurônios e admitindo o número máximo de iterações igual a 10000 e valor de corte (threshold) igual 1, a programação para estimar os parâmetros da rede é:
set.seed(12345)
library("neuralnet")
maxit<-as.integer(1000000)
nn <- neuralnet(dados.Treina[,5]~dados.Treina[,4]+
+dados.Treina[,3]+dados.Treina[,2]+dados.Treina[,1],
+data=dados.Treina, hidden=7,threshold =1,stepmax= maxit)
onde o comando fixa a semente para a geração de números aleatórios com o intuito de permitir a reprodução dos resultados obtidos.
CUIDADO!!
A programação acima pode demorar para ser concluída, por isso, use com parcimônia os valores do número máximo de iterações e parâmetro threshold.
Uma vez obtidos os valores para os parâmetros do modelo de redes neurais podemos analisar os resultados usando os comandos:
#Apresenta os valores para os pesos
nn$result.matrix

#Faz o gráfico do modelo
plot(nn)
O próximo passo é realizar a previsão para os dados de validação:
#Faz a previsão
previsao<-compute(nn,dados.Valida[,1:4])

#Valores da previsao
previsao.nn<-previsao$net.result
Outra hipótese bastante comum para as séries temporais financeiras é a Hipótese do Passeio Aleatório a qual pode ser representada matematicamente como: $y_{t}=\mu+y_{y-1}+\epsilon_{t}$ Assumindo que $\epsilon_{t}\sim N(0,\sigma^{2})$ podemos utilizar os dados de treinamento para estimar os parâmetros desse modelo:
#Gera as estimativas para o Random Walk
epsilon<-(dados.Treina[,5]-dados.Treina[,4])
mu<-mean(epsilon)
sigma2<-var(epsilon)

#Faz a previsao usando o Random Walk
previsao.rw<-dados.Valida[,4]+rnorm(nrow(dados.Valida),
+mu,sqrt(sigma2))
O gráfico comparativo entre a previsão via redes neurais e Passeio Aleatório contra os valores observados é elaborado usando o seguinte código:
#Monta a base com as previsões
Tempo<-seq(1458,2917)
previsao.todos<-as.data.frame(cbind(
+Tempo,dados.Valida[,5],previsao.nn,previsao.rw))
colnames(previsao.todos)<-c("Tempo","Obsevado",
+"Predito.NN", "Predito.RW")

#Faz o gráfico
#Faz o gráfico
library("ggplot2")
ggplot(previsao.todos, aes(Tempo)) + 
  geom_line(aes(y = Obsevado, colour = "Obsevado")) + 
  geom_line(aes(y = Predito.NN, colour = "Predito NN"))+
  geom_line(aes(y = Predito.RW, colour = "Predito RW"))+
  ggtitle("Séries Temporais")
É interessante notar que o modelo de Passeio Aleatório se ajustou melhor aos valores observados do que o modelo de redes neurais, especialmente quando a volatilidade apresentada era grande. Comparando a distribuição dos erros entre os dois métodos temos:
#Encontra os erros
erro.rw<-previsao.todos[,2]-previsao.todos[,4]
erro.nn<-previsao.todos[,2]-previsao.todos[,3]

#Une os dados
library("reshape")
df.m <- melt(cbind(erro.rw,erro.nn))

#Calcula as funções densidade via Kernel.
ggplot(df.m) + geom_density(aes(x = value,colour = X2)) 
+ labs(x = "Valores",y="Densidade") +
ggtitle("Densidade por meio do estimador Kernel.")
O gráfico de densidades demonstra que o modelo Passeio Aleatório apresenta menos variabilidade na previsão do que o modelo de redes neurais. Uma proposta seria incorporar ao modelo de redes neurais alguma covariável associada a volatilidade e comparar novamente os resultados ou ainda, adicionar mais camadas ocultas no modelo com o intuito de aumentar a não-linearidade.

quinta-feira, 8 de novembro de 2012

Otimização de portfólio.


O pacote do R denominado Tawny fornece uma maneira simples para otimizar carteiras de investimento de forma a minimizar o risco em uma carteira.

Este otimizador pode executar a otimização de carteiras para diversos períodos temporais simultaneamente.

A ideia é: "não colocar todos os seus ovos em uma mesma cesta".

Para fins de ilustração, considere um subconjunto de ativos do S&P 500. Abaixo estão os códigos em R para otimizar a carteira de investimentos:

#Habilita a biblioteca Tawny 
library(tawny)

#Lê o subconjunto de dados do S&P 500
data(sp500.subset)
dados <- create(TawnyPortfolio, sp500.subset, window=190)

# Otimizando por meio de uma janela temporal de tamanho 190
pesos <- optimizePortfolio(dados, create(SampleFilter) )

quarta-feira, 3 de outubro de 2012

Testando a hipótese de reação exagerada em mercados.


Até meados da década de 1980, pesquisas econômicas financeiras eram essencialmente sobre como investigar a interação entre o comportamento racional dos maximizadores de utilidade.

Muito sobre essa abordagem de pesquisa e investigação foi produzido mas pouco sobre o comportamento individual dos participantes do mercado foi explicitado. Com as influências da psicologia clínica e economia experimental impactando a disciplina de finanças a partir de ângulos diferentes, a investigação sobre os aspectos comportamentais em finanças começaram a surgir.

Um dos primeiros artigos empíricos na área de finanças comportamentais foi o estudo seminal de DeBondt e Thaler (1985) sobre a reação exagerada do mercado de ações.

Os resultados DeBondt e Thaler (1985), mostraram que quando os ativos eram mantidos por longos períodos de tempo, as ações apresentavam um comportamento significante para a reversão dos preços desses mesmos ativos.

A explicação para a quantidade substancial de reversão no preço dos ativos é baseada no comportamento. Os investidores, em geral, tendem a "super-ponderar" desempenhos mais recentes do ativo em questão e "sub-ponderar" informações de mais longo prazo em suas decisões. Esse comportamento é consistente com a heurística da disponibilidade apresentada por Tversky e Kahneman (1973).

Por exemplo, se uma ação sofreu recentemente uma diminuição nos ganhos devido a alguma mudança não estrutural temporária nas oportunidades econômicas, a queda dos preços pode ser exagerada devido a atualidade da notícia.

Uma vez que o preço da ação foi pressionado, o seu aumento gradual até a sua recuperação completa pode ser lento. Isso até que os investidores percebam sua "reação exagerada" à notícia temporariamente ruim, criando assim uma pressão compradora para "inverter" a queda dos preços.

DeBondt e Thaler (1985) encontraram uma persistência significante para esse fenômeno de reversão e sugeriram que os mercados tendem a sofrer de reação exagerada dos investidores.

A beleza da abordagem de DeBondt e Thaler (1985) está em demonstrar como uma hipótese de reação exagerada pode ser testada empiricamente.

Especificamente, o que DeBondt e Thaler (1985) fazem é calcular o retorno das carteiras para alguns períodos, por exemplo, o desempenho da carteira pode ser calculada para os 36 meses anteriores e, posteriormente, avaliamos como o portfólio se comporta, especialmente para um determinado período de avaliação ou seja, o desempenho nos 36 meses subsequentes.

Para testar a hipótese de reação exagerada, DeBondt e Thaler (1985) formaram portfólios "perdedores" e "vencedores" compostos pelas 30 ações com o pior desempenho e 30 ações com o melhor desempenho, respectivamente, durante o período de formação.

Os autores, então, seguem estas carteiras para os 36 meses seguintes e avaliam e comparam seus desempenhos. Com o tempo, a carteira "perdedora", composta das ações com o pior desempenho no período de formação, consistentemente supera a carteira de "vencedora" durante o período de avaliação. Esses resultados, de acordo com DeBondt e Thaler (1985), sugerem uma forte evidência da presente reação exagerada do mercado de ações, consistente com a teoria da tomada de decisão irracional dos agentes.

O objetivo deste post é mostrar como você pode realizar um simples estudo sobre a reação exagerada do Mercado de ações usando o R para isso.

Teste de reação exagerada usando o R.


DeBondt e Thaler (1985) acumularam em seu trabalho os retornos ajustados do mercado ao longo de um período de formação de três anos para cada ação do New York Stock Exchange (NYSE) com início em dezembro de 1932 e continuando até dezembro de 1977. As carteiras "perdedoras" e "vencedora" foram as carteiras avaliadas ao longo dos 17 períodos disjuntos de três anos de avaliação.

No nosso exemplo, trabalharemos com os ativos disponíveis de 2000 a 2012 e assim como DeBondt e Thaler (1985) utilizaremos três anos para o período de formação e três anos para o período de avaliação.

Conrad e Kaul (1993) introduziram a ideia de que os retornos acumulados em horizontes longos podem enviesar os resultados sobre a super-reação dos agentes. Uma solução para isso é usar períodos de retenção (holding periods) ao investigar horizontes longos.

Tal como acontece com a maior parte da investigação em economia financeira, o primeiro passo para a análise de dados é lê-los a partir de um conjunto de observações de dados externos. Esta tarefa é realizada pela seguinte sintax:

#Limpa o console
rm(list=ls(all=TRUE))

#Lista de ativos Bovespa
ativos<-read.csv("http://dl.dropbox.com/u/36068691/Ativos.csv")

#Obtêm os dados da Bovespa
lista<-as.character(ativos[,"Sigla"])
Em seguida, captura-se a data atual da máquina para a composição da série de interesse:
#Pega a data atual
d <- Sys.time()
current.year <- format(d, "%Y")
current.month <- format(d, "%m")
current.day <- format(d, "%d")
A fonte dos nossos dados é a API do Yahoo Finance, pior isso, precisamos especificar algumas parametrizações:
#Yahoo finance considera janeiro como 00 então:
month <- as.numeric(current.month) - 1
month <- ifelse(month < 10, paste("0",month, sep=""), m)

#Pega os prerídos diários para o ticket 1
tckr <- lista[1]                                   
fn <- sprintf("http://ichart.finance.yahoo.com/table.csv?s=%s&a=0&b=1&c=2000&d=%s&e=%s&f=%s&g=d&ignore=.csv",tckr, month, current.day, current.year)

#Guarda os dados
dados <- read.csv(fn, colClasses=c("Date", rep("numeric",6)))
dados<-dados[c("Date","Close")]
colnames(dados)<-c("Date",lista[1])
O código acima captura os dados do primeiro ativo listado e armazena em um data.frame denominado dados. Em seguida, precisamos repetir esse processo para cada um dos ativos disponíveis:
for(i in 2:length(lista))
{
  tckr <- lista[i]                                
  fn <- sprintf("http://ichart.finance.yahoo.com/table.csv?s=%s&a=0&b=1&c=2000&d=%s&e=%s&f=%s&g=d&ignore=.csv",tckr, month, current.day, current.year)
  
  #Guarda os dados
  data<-data.frame()
  data <- try(read.csv(fn, colClasses=c("Date", rep("numeric",6))))
  if(ncol(data)>0)
    {
      data<-data[c("Date","Close")]
      colnames(data)<-c("Date",tckr)
      dados<-merge(dados, data, by.x = "Date", by.y = "Date", all = TRUE)
    }
}
Uma função muito útil para se utilizar em análise de séries temporais é a função shift desenvolvida pelo R-Bloggers. Essa função é capaz de defasar a série temporal para trás (lag) e para frente (lead):
shift<-function(x,shift_by){
  stopifnot(is.numeric(shift_by))
  stopifnot(is.numeric(x))
  
  if (length(shift_by)>1)
    return(sapply(shift_by,shift, x=x))
  
  out<-NULL
  abs_shift_by=abs(shift_by)
  if (shift_by > 0 )
    out<-c(tail(x,-abs_shift_by),rep(NA,abs_shift_by))
  else if (shift_by < 0 )
    out<-c(rep(NA,abs_shift_by), head(x,-abs_shift_by))
  else
    out<-x
  out
}
Então, para cada ativo calculamos os retornos nos períodos de formação (três anos antes) e no período de avaliação (três anos depois):
#Para cada ativo calcular o retorno três anos antes e três anos depois
resultado<-data.frame()
for(i in 1:length(lista))
{
  precos<-dados[,c("Date",lista[i])]
  precos[,"pform"]<-shift(precos[,lista[i]],-3*360)
  precos[,"paval"]<-shift(precos[,lista[i]],+3*360)  
  precos[,"pform"]<-(precos[,lista[i]]-precos[,"pform"])/precos[,"pform"]
  precos[,"paval"]<-(precos[,"paval"]-precos[,lista[i]])/precos[,lista[i]]
  precos$Empr<-lista[i]
  precos<-na.omit(precos[,c("Date","Empr","pform","paval")])
  resultado<-rbind(resultado,precos)
}
Podemos em seguida eliminar as observações que são outliers:
#Remove retornos exagerados
resultado<-resultado[resultado[,"pform"]>-1 & resultado[,"paval"]>-1, ]

#Ordena os resultados
resultado <- resultado[order(resultado[,"Date"], resultado[,"Empr"]), ]

#Frequência dos ativos 
firmas <-  as.data.frame(table(resultado[,"Date"]))
names(firmas) <- c("Date", "freq")
firmas[,"Date"]<-as.Date(firmas[,"Date"])

#Une as frequências com a base resultado
dat <- merge(resultado, firmas, by="Date",all.x=T)
O código acima além de fazer a limpeza dos outliers calcula a frequência de dados para cada ativo e une esse resultado no data.frame dat. Os ativos são então ranqueados permitindo assim a definição dos "perdedores" e dos "vencedores":
#Ranqueia as empresas
firma <- unlist(tapply(rep(1, nrow(resultado)), resultado[,"Date"], cumsum))
dat <- cbind(dat, firma)
O teste T é então utilizado para se comparar o comportamento das carteiras "perdedoras" e dos "vencedoras" e consequentemente a possível existência de uma inversão das categorias "perdedoras" e "vencedoras":
#Carteira perdedora
perdedor <- dat[dat[,"firma"]<=5, ]
summary(perdedor[,"paval"])
t.test(perdedor[,"paval"])

#Carteira vencedora
vencedor <- dat[dat[,"firma"] >= dat$freq-4, ]
summary(vencedor[,"paval"])
t.test(vencedor[,"paval"])
A carteira "perdedoras" apresenta intervalo de confiança de 95% como [0.5463877, 0.6533983] e estimativa pontual de 0.599893. Já a carteira dita "vencedora" apresenta intervalo de confiança de 95% como [1.068938, 1.364774] e estimativa pontual de 1.216856, dessa forma não há evidências de reação exagerada para o conjunto de dados e período estudado.

terça-feira, 18 de setembro de 2012

Teste de razão de variâncias de Lo-MacKinlay.


Existem muitos textos que sugerem através de evidências empíricas que o retorno de ações contém componentes que podem ser preditos.

Por exemplo, Keim e Stambaugh (1986) encontraram um modelo estatisticamente significante para prever o preço de ações usando um modelo de previsão com variáveis pré-determinadas.

Fama e French (1987) mostraram que um longo período de retenção de retornos é negativamente correlacionado serialmente, implicando assim que 25% a 40% da variação de um longo período de retornos é previsível com base nos retornos passados.

Esse post, no entanto, tem como objetivo principal apresentar o Teste de Razão de Variâncias proposto por Lo e MacKinlay (1988). Os autores indicam que o modelo de passeio aleatório (ou passeio do bebâdo) é, em geral, inconsistente (para os ativos estudados) com o comportamento estocástico dos retornos semanais, especialmente para as ações com menor capitalização (Small Caps).

Lo e MacKinlay (1988) alertam que estes resultados não implicam necessariamente que o mercado de ações é ineficiente ou que os preços das ações não são avaliações racionais dos seus valores "fundamentais". O teste proposto pelos autores pode ser interpretado como um teste para a rejeição de "algum" modelo econômico de formação eficiente de preços.

Teste de razão de variâncias.

A eficiência do mercado de ações tem sido debatida por pesquisadores e profissionais do mercado financeiro. Uma forma particular de eficiência do mercado de ações é a eficiência informacional.

A eficiência informacional se baseia na premissa de que os preços dos ativos refletem as informações relevantes disponíveis instantaneamente aos investidores e ao público em geral. Como a chegada de informações é imprevisível, os preços dos ativos também se tornam imprevisíveis. A hipótese de que a chegada de informações fundamentais ao mercado é aleatória e, portanto, os movimentos dos preços dos ativos também serão aleatórios, é englobada pela teoria do passeio aleatório dos preços dos ativos.

De maneira simplista, a teoria do passeio aleatório indica que uma vez que a chegada de informações é imprevisível, o melhor preditor do preço de um ativo é seu valor atual. Esta ideia simples pode ser facilmente incorporada no modelo de passeio aleatório bem conhecido dos preços dos ativos, que pode ser expresso da seguinte forma: $P_{t}=P_{t-1}+\epsilon$.

$P_{t}$ é o preço atual do ativo em estudo, $P_{t-1}$ é o preço do período anterior, e $\epsilon$ é um termo de erro aleatório. Cada termo de erro aleatório representa a chegada de uma nova informação, o qual se assumido ser imprevisível deve ser independentes uns dos outros. Admitindo sobre a hipótese nula o termo de erro aleatório é independente e identicamente distribuído segundo uma distribuição normal, então a variância do termo de erro aleatório é linear no intervalo de tempo durante o qual os preços são observados. Simplesmente, a variância da variação de preços quinzenais deve ser o dobro da variação de preços semanais. Além do mais, a variância das alterações de preços mensais deve ser quatro vezes superior ao de alterações de preços semanais, e assim por diante.

A relação linear entre o intervalo de tempo observado para os preços do ativo de interesse e a sua variação é a essência da especificação do teste desenvolvido por Lo e MacKinlay (1988). Lo e MacKinlay (1988) desenvolveram distribuições limitantes para os estimadores de razão de variância, com e sem a existência de heterocedasticidade, e mostraram que os preços dos ativos não necessariamente seguem um passeio aleatório. Seus estimadores são definidos como se segue:

A equação (01) representa a média das alterações dos preços em T períodos de tempo:


Já a equação (02) representa um estimador de variância para as alterações dos preços no period de tempo especificado, isso é $t=1,\dots,T$:


Por fim, a equação (03) representa um estimador de variância para as mudanças q-temporais de preços:


e $m=q(T-q+1)(1-\frac{q}{T})$ é um ajuste executado no denominador da variância do estimador q-temporal para acomodar as observações que se sobrepõem, além de ajudar a aumentar a potência do teste de razão de variância.

Especificação do teste.


Assim, o teste da razão de variâncias é definido por:


Para acomodar a possível heterocedasticidade, utiliza-se a estatística do teste padronizado $z$ o qual é assintoticamente distribuído segundo uma normal padrão:



onde


e


Aplicação usando os dados da Bovespa


Utilizando o R, vamos calcular o teste de razão de variâncias de Lo e MacKinlay para as ações PETR4 em um período $q=1$ diário, isso é, com defasagem de um dia.

#Limpa o console
rm(list=ls(all=TRUE))

#Habilita o pacote quantmod
library(quantmod)

#Início do período de interesse
start = as.Date("2011-01-01") 

#Fim do período de interesse
end = as.Date("2011-12-31") 

#Obtêm os dados da PETR4
getSymbols("PETR4.SA", src="yahoo",from=start,to=end)

A seguir é necessário calcular os log-retornos de duas séries: a série original e a série defasada (nesse caso em um dia):

#Aplica o log-retorno nos preços
ret <- as.matrix(c(rep(NA, 0), diff(log(Cl(PETR4.SA)), lag=1)))

#Aplica o log-retorno nos preços (Defasagem diária)
ret2 <- as.matrix(c(rep(NA, 0), diff(log(Cl(PETR4.SA)), lag=2)))
#Plota a série temporal dos preços de fechamento
plot(ret, main="Petr4",type="l",
+ylab="Log-return Close price. (Lag 1)",xlab="Time")

plot(ret2,main="Petr4",type="l",
+ylab="Log-return Close price. (Lag 2)",xlab="Time")
Em seguida os parâmetros exigidos para o teste são calculados:
#Calcula muhat
muhat <- mean(ret, na.rm=TRUE)
nq <- nrow(ret)
#Parâmetros Sigma e Delta
sigatop <- (ret - muhat)^2
sigatop1 <- c(NA, sigatop[-length(sigatop)])
deltop <- sigatop * sigatop1
delbot <- sigatop
sigctop <- (ret2 - 2*muhat)^2
#Guarda os resultados
dados <- data.frame(ret, sigatop, sigatop1, deltop, delbot, sigctop)
Fazendo a soma, obtemos os valores desejados:
#Calcula a soma dos valores
sigatop <- sum(sigatop, na.rm=TRUE)
sigctop <- sum(sigctop, na.rm=TRUE)
deltop <- sum(deltop, na.rm=TRUE)
delbot <- sum(delbot, na.rm=TRUE)
Em seguida o cálculo dos testes é realizado:
nq<-nrow(ret)
q<-2
qm1<-q-1
theta<-0
j<-1
m <- q*(nq-q+1)*(1-q/nq)
siga <- sigatop/(nq-1)
sigc <- sigctop/m
varrat2 <- sigc/siga
delta <- nq*deltop/(delbot**2)
while(j <= qm1)
{
  theta <- theta + ((2*(q-j)/q)**2)*delta
  j<-j+1
}

z <- sqrt(nq)*(varrat2-1)/sqrt(theta)
print(nq, cat("Número de retornos:"))
print(varrat2, cat("Razão de variâncias para retronos com Lag=2:"))
print(z, cat("Teste Robusto a Heterocedasticidade de Lo-MacKinlay:"))
Após realizar a análise obtêm-se:
  • Número de retornos:[1] 249
  • Razão de variâncias para retronos com Lag=2:[1] 1.081269
  • Teste Robusto a Heterocedasticidade de Lo-MacKinlay:[1] 1.052326
Os resultados indicam que a hipótese nula de presença de um comportamento como um passeio aleatório para a série PETR4 diária entre 01-01-2011 até 31-12-2011 não pode ser rejeitada, uma vez que o valor do teste não é maior do que o valor crítico sobre a hipótese nula, qual seja, aproximadamente 1.96 a um nível de confiança de 95%. Uma maneira direta de se atingir o mesmo resultado é por meio da biblioteca vrtest. A função que realiza o teste de Lo-MacKinlay é a função Lo.Mac(y=, kvec=) onde o parâmetro y= recebe o vetor de ativos e kvec= a defasagem utilizada no teste.

sábado, 8 de setembro de 2012

Calculando a função densidade de probabilidade da distribuição estável usando FFT.


Distribuições estáveis ​​paretianas têm propriedades atraentes para modelagem empírica em finanças, porque incluem a distribuição normal como um caso especial, mas também pode permitir caudas mais pesadas e assimetria.

Uma razão principal para a pouca utilização dessa distribuição em trabalhos acadêmicos aplicados é devido ao fato de que, em geral, não há expressão de "forma fechada" para a a função de densidade de probabilidade, e que as aproximações numéricas computacionais são não-triviais e computacionalmente extensivas.

Nesse post vou mostrar como é possível calcular a função densidade de probabilidade via Fast-Fourier Transform (FFT).

O trabalho original sobre esse assunto foi produzido por Mittnik, Doganoglu e Chenyao (1999).

A distribuição alfa-estável.


A distribuição alfa-estável, em geral, não possui expressão analítica para sua função densidade de probabilidade (f.d.p) ou ainda para a sua função distribuição acumulada (f.d.a), mas pode ser escrita por meio de sua função característica (Rachev e Mittnik, 2000 ):


onde $0<\alpha<2$ é o expoente da distribuição ou índice de cauda, $-1\le \beta\le 1$ é o parâmetro de assimetria, $\sigma>0$ é o parâmetro de escala e $\delta \in \Re$ é o parâmetro de locação, a função $\omega(.,.)$ é dada por:




A distribuição alfa-estável representada acima e definida pela notação $S_{\alpha}^{0}(\beta,\delta,\sigma)$ é denominada parametrização $S_{0}$ segundo Nolan (2010).

A função densidade de probabilidade pode ser aproximada utilizando o método FFT (Fast Fourier Transform) o qual é computacionalmente eficiente e permite um processo de aproximação mais rápido do que expansão por séries (Bergström, 1952) ou integração direta (Nolan, J. P., 2001. Maximum likelihood estimation of stable parameters. Manuscrito não publicado.).

Segundo Durrett (2010) página 106 uma função densidade de probabilidade pode ser escrita pela Transformada de Fourier da função característica, em outras palavras:


A integral acima pode ser calculada para $n$ pontos igualmente espaçados com distância $h$ e soma resultante pode ser computada por meio do método FFT (Fast Fourier Transform). Mittnik e Doganoglu (1999) sugerem que os valores de $n$ e $h$ devem ser respectivamente $2^{13}$ e $h=0.01$ para que uma boa aproximação seja possível.

Implementação no R.


Podemos então implementar o método descrito anteriormente no R.

Inicialmente, definimos os valores dos parâmetros para os quais desejamos calcular a função densidade de probabilidade:

#Limpa o console
rm(list=ls(all=TRUE))

#Parâmetros para o calculo da PDF
alpha<-0.5
beta<-0.3
sigma<-1
delta<-0
O limite superior para o suporte da distribuição $xmax$ e o número de pontos utilizados para a aproximação $n$ deve ser definido:
#Valor máximo em módulo da PDF
xmax<-20

#2^N espaços iguais entre [-xmax,xmax]
n<-13
Devido à natureza do método FFT, valores distantes do centro da distribuição podem ser subestimados. Por esta razão o método aqui sugerido calcula a função densidade de probabilidade da distribuição alfa-estável para o intervalo $[-xmax, xmax] \times 2^{mult}$ onde o multiplicador $mult$ é utilizado para aumentar (temporariamente) o suporte da distribuição, em seguida, esse intervalo aumentado é truncado em seu intervalo original. Seja o valor $mult=4$ como uma primeira tentativa, no entanto, para melhor precisão deve-se utilizar $mult>4$.
#Define o multiplicador
mult<-4

#Calcula o intervalo do suporte aumentado
xmax <- xmax*(2^mult)

#Aumenta n para cobrir o intervalo maior com os mesmos pontos do grid
n <- n + mult
M <- 2^n
R <- pi/xmax
dt <- 1/(R*M)
Como o método FFT admite a utilização de números complexos e devido a função característica da distribuição alfa-estável possuir número imaginário na sua estrutura precisamos definí-lo no R:
#Número imaginário
i<-complex(real = 0, imaginary = 1)

#Define o grid de valores com pontos uniformemente espaçados 
xx <-seq(from = -2^(n-1)+0.5, to = (2^(n-1)-0.5))/(2^n*dt)
Em seguida, a função característica é definida:
#Função característica da distribuição alfa-estável
piby2 <- pi/2
if(abs(alpha-1)<0.001)
  {
    yy <- exp( -sigma*(abs(xx))*( 1+i*beta*sign(xx)/piby2*log(sigma*abs(xx))) + i*delta*xx )
  }
else
{
    yy <- exp( -(sigma*abs(xx))^alpha*( 1+i*beta*sign(xx)*tan(alpha*piby2)*( (sigma*abs(xx))^(1-alpha)-1 ) ) + i*delta*xx )
}
Note que caso $\alpha$ esteja muito próximo de um, admitimos que $\alpha=1$ e utilizamos a estrutura da função característica nessa condição. O próximo passo é aplicar o método FFT:
#Utiliza o método FFT
yy1 <- c(yy[seq(from=(2^(n-1)+1),to=2^n)], yy[seq(from=1,to=2^(n-1))])
z <- Re( fft(yy1) )/(2*pi)*R
E após a utilização do método FFT pela função fft do R, a função densidade de probabilidade é calculada:
#Computa a função densidade de probabilidade
x <- (2*pi)*(seq(from=0,to=(M-1),by=1)/(M*R)-1/(2*R))
y <- c(z[seq(from=(2^(n-1)+1),to=2^n)], z[seq(from=1,to=2^(n-1))])
Por fim, o intervalo original do suporte da distribuição é obtido e o esboço de sua função densidade de probabilidade é apresentado:
#Encontra os pontos da f.d.p no intervalo original
T <- which((x<=xmax/(2^mult)) & (x>=-xmax/(2^mult)))
x <- x[T]
y <- y[T]

#Plota a função densidade de probabilidade
plot(x,y,type="l", main ="Distribuição alfa-estável"
+,xlab="Suporte",ylab="f.d.p",col="blue")
O valor da função densidade de probabilidade para um ponto arbitrário qualquer pode ser obtido por meio de interpolação entre os pontos que contém o valor desejado. Já a aproximação via FFT para outras parametrizações segue o mesmo processo, exceto pela definição da expressão da função característica na sintaxe apresentada anteriormente que deve ser alterada para a parametrização desejada.