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.

segunda-feira, 15 de janeiro de 2018

Leitura das bases do IBGE por meio do microdadosBrasil e dplyr.


Como pesquisadores uma fonte de informação rica e muito importante para a academia são as bases de dados do Instituto Brasileiro de Geografia e Estatística (IBGE).

A leitura das principais bases de dados nacionais pode ser realizada por meio do pacote microdadosBrasil.

O primeiro passo é a instalação do pacote no R, fazemos isso da seguinte forma:

#Habilita o pacote devtools
library(devtools)
#Instala o pacote microdadosBrasil
devtools::install_github("lucasmation/microdadosBrasil")
#Habilita o pacote microdadosBrasil
library(microdadosBrasil)

Nesse post, trataremos da Pesquisa Nacional por Amostra de Domicílios (PNAD) como exemplo, para baixar e ler a PNAD precisamos dos seguintes comandos:

#Fazemos somente uma vez caso os dados não existam localmente
download_sourceData("PNAD", 2015, root_path = "C:\\Dados")
#Lê os dados da PNAD
pnad <- read_PNAD("pessoas", 2015, root_path = "C:\\Dados")
#Salva a PNAD lida
save.image("C:\\Dados\\pnad.RData")
No bloco anterior os microdados, dicionários, arquivos de leitura da PNAD são armazenados na pasta definida pelo usuário (nesse exemplo C:\Dados). Em seguida esses dados são lidos pela função read_PNAD e então salva a imagem no arquivo pnad.RData. No código anterior definimos ainda que desejamos os dados da base de pessoas, para o ano de 2015. Maiores informações sobre os tipos de bases e pesquisas disponíveis os interessados podem consultar o repositório do microdadosBrasil. Uma vez salva a imagem, a leitura fica muito mais rápida no R e podemos chamar os dados por meio da função load:
#Chama os dados da pnad salvos no RData
load("C:\\Dados\\pnad.RData")
Como exemplo vamos calcular a renda média das pessoas na PNAD 2015. A variável rendimento de interesse é dada pela variável V9532 e a variável de peso (expansão amostral) é dada por V4729. O cômputo do rendimento médio é dado por?
#Chama a biblioteca dplyr
library(dplyr)
#Exemplo Média do Rendimento
pnad %>%
  group_by(UF) %>%
  filter(V9532<999999999)%>%
  summarise(weighted.mean(V9532 , 
                          w = V4729, 
                          na.rm = TRUE)
  ) -> agg.df
No bloco anterior salvamos no data.frame agg.df a renda média das pessoas por Unidade da Federação (UF) removendo-se os valores missing informados pelo IBGE como V9532>999999999. Outra estatística muito comum é o cálculo de proporções. Podemos calcular a proporção por cor nas Unidades da Federação (UF) como:
#Calcula a proporção
pnad %>%
  group_by(UF, V0404) %>%
  summarise (n.cor = sum(V4729))%>%
  group_by(UF) %>%
  mutate (n = sum(n.cor)) %>%
  mutate(freq = n.cor/n) %>%
  mutate(freq2 = paste0(round(freq,6)*100,"%")) -> cor.df
#Define os labels dos factors
cor.df[,"Cor"] = factor(cor.df$V0404,
                    levels=c(2,4,6,8,0,9),
                    labels=c("Branca","Preta","Amarela",
                             "Parda","Indígena","Sem declaração"),
                    ordered=FALSE)
Nesse exemplo, salvamos no objeto cor.df essas proporções e também definimos os labels das variáveis qualitativas. Na última parte desse exemplo, podemos estar interessados no cômputo do Índice de Gini por Unidade da Federação:
#Exemplo Gini
library(reldist)
pnad %>%
  group_by(UF) %>%
  filter(V9532<999999999)%>%
  summarise(gini.idx=gini(V9532,V4729))
Onde o pacote reldist apresenta a função gini que calcula o Índice de Gini ponderado pelo peso amostral.

sexta-feira, 15 de dezembro de 2017

Introdução à Inferência Variacional.


A Inferência Variacional é uma abordagem para se estimar distribuições a posteriori muito difíceis ou complexas.

Frequentemente, em Inferência Bayesiana estamos interessados em computar as medidas $p(z|x)$ para variáveis latentes $z_{1},\dots,z_{n}$ dado as observações $x_{1},\dots,x_{n}$. Quando essa distribuição é intratável podemos aproximar $p(z|x)$ por $q(z)$ onde a medida de proximidade é dada pela divergência de Kullback-Leibler.

Usualmente a Inferência Variacional é realizada através dos seguintes passos:

  1. Comece com um modelo para $p(z,x)$.
  2. Encontre a aproximação variacional $q(z|\nu)$
  3. onde $\nu$ são os parâmetros variacionais.
  4. Encontre a Evidence Lower Bound (ELBO) dada por $L(\nu)=E_{Q}[\log(p(z,x))-\log(q(z|\nu))]$
  5. Minimize a quantia $L(\nu)$
  6. Volte para o passo 3 e faça isso até haver convergência.

Com base no exposto, vamos incialmente considerar um modelo de Regressão Probit construído da seguinte forma:

#Define a semente
set.seed(8333)
#Gera os dados
x<-rnorm(1000)
#Gera a latente
z<- rnorm(1000,10+5*x,1)
#Gera a variável dependente:
y<-ifelse(z<0,0,1)
O objetivo é estmar os parâmetros $a=10,b=5$ do modelo de Regressão Probit. Note que o Processo Gerador de Dados temos que a distribuição a posteriori é dada por: $p(a,b,z|y,x)= \prod_{i=1}^{n}y_{i}^{I(z_{i}>0)}(1-y_{i})^{1-I(z_{i}>0)}\times \frac{1}{\sqrt{2\pi}}\exp\left[-\frac{(z_{i}-(a+bx_{i}))^{2}}{2}\right]$ o qual também pode ser escrito como: $p(a,b,z|y,x)\propto \prod_{i=1}^{n}y_{i}^{I(z_{i}>0)}(1-y_{i})^{I(z_{i}\leq 0)}\times \exp\left[-\frac{(z_{i}-(a+bx_{i}))^{2}}{2}\right]$ para distribuições a priori uniformes para $a,b$. Dado o modelo precisamos definir a aproximação variacional. Nesse caso Blei et. al. (2017) sugere como regra ótima para a aproximação variacional fazer: $q(\theta_{j})\propto \exp\left[E_{\theta_{-j}}(\log(p(\theta|y,x)))\right]$ onde $\theta=(\theta_{1},\dots,\theta_{p})$ são os parâmetros e variáveis latentes do modelo. Nesse sentido, temos: $\log(p(a,b,z|y,x))=\sum_{i=1}^{n}I(z_{i}>0)\log(y_{i})+I(z_{i}\leq 0)\log(1-y_{i})-\sum_{i=1}^{n}\frac{(z_{i}-(a+bx_{i}))^{2}}{2}$ A distribuição varicional ótima é dada por: $q_{j}(z_{j})\propto \exp\left[E_{a,b,z_{-j}}(\log(p(a,b,z|y,x)))\right]$ Note que: $E_{a,b,z_{-j}}(\log(p(a,b,z|y,x)))=I(z_{j}>0)\log(y_{j})+I(z_{j}\leq 0)\log(1-y_{j})-\frac{E_{a,b}(z_{j}-(a+bx_{j}))^{2}}{2}$ Aplicando a exponencial temos que a distribuição ótima para $z_{j}$ é uma distribuição Normal Truncada, isso é $q_{j}(z_{j})\sim N^{+}(E(a)+E(b)x_{j},1)$ se $y_{j}=1$ e $q_{j}(z_{j})\sim N^{-}(E(a)+E(b)x_{j},1)$ se $y_{j}=0$. Similarmente, para $a$ temos $E_{b,z}[\log(p(a,b,z|y,x))]\propto E_{b,z}\left(-\frac{\sum_{i=1}^{n}(z_{i}-(a+bx_{i}))^{2}}{2}\right)$ , removendo os termos que não dependem de $a$ e complentando quadrados temos: $q_{a}(a)\sim N(\frac{\sum_{i=1}^{n}[E(z_{i})-E(b)x_{i}]}{n},\frac{1}{n})$. Finalmente para $b$ temos $E_{a,z}[\log(p(a,b,z|y,x))]\propto E_{a,z}\left(-\frac{\sum_{i=1}^{n}(z_{i}-(a+bx_{i}))^{2}}{2}\right)$, novamente removendo os termos que não dependem de $b$ e complentando quadrados temos: $q_{b}(b)\sim N(\frac{\sum_{i=1}^{n}[E(z_{i})-E(a)]x_{i}}{\sum_{i=1}^{n}x_{i}^{2}},\frac{1}{\sum_{i=1}^{n}x_{i}^{2}})$. Nesse caso, a otimização ocorre atribuindo os valores dos parâmetros a sua respectiva média da distribuição variacional. As equações para atualização dessas distribuições são dadas por:
#Equações de Update
update.newZj = function(newA,newB,j) {
  mu = newA + newB*x[j]
  if (y[j] == 1) {
    return(mu + dnorm(-1*mu)/(1-pnorm(-1*mu)))
  } else {
    return(mu - dnorm(-1*mu)/(pnorm(-1*mu)))
  }
}
update.newA = function(newZ,newB) {
  return(sum(newZ-newB*x)/n)
}
update.newB = function(newZ,newA) {
  return(sum(x*(newZ-newA))/sum(x^2))
}
Finalmente, podemos construir um sistema iterativo baseado nos valores otimizados:
#Total de observações
n<-length(z)
#Chute inicial
newA<-3
newB<-9
#Inicializa vetores
a.vec<-rep(NA,n)
b.vec<-rep(NA,n)
newZ<-rep(NA,n)
#Inicia o processo iterativo
for(it in 1:1000){
  #Update Z
  for (i in 1:n) {
    newZ[i] = update.newZj(newA,newB,i)
  }
  #Correct Numerical
  newZ[is.infinite(newZ)]<-100
  #Update A
  newA = update.newA(newZ,newB)
  newB = update.newB(newZ,newA)
  a.vec[it] = newA
  b.vec[it] = newB
}


plot(a.vec,type="l")
summary(a.vec)
plot(b.vec,type="l")
summary(b.vec)

segunda-feira, 20 de novembro de 2017

Usando a função DEoptim em paralelo.


Um dos pacotes mais interessantes para realizar otimização é o pacote DEoptim. Essa abordagem permite a otimização de funções dos mais diversos tipos de uma maneira fácil.

Entretanto, por usar Otimização Evolucionária o processo de otimização pode ser lento, caso utilize-se muitas populações, iterações ou níveis tolerância extremos.

Aqui mostrarei como usar a opção parallelType da função DEoptim para que esse processo seja acelerado por meio da programação em paralelo.

Considere a seguinte função com múltiplos máximos e mínimos locais:

f<-function(x,y){
  res<-(268/((3*x+18)^2+(3*y+7)^2+(18+7)))-(2*(268)/((3*x-21)^2+(3*y+21)^2+(21+21)))-(3*(268)/((3*x+18)^2+(3*y-21)^2+(18+21)))+(4*(268)/((3*x-21)^2+(3*y-7)^2+(21+7))) 
  return(res)
}
Podemos fazer o gráfico dessa função através do seguinte bloco de comandos:
#Chama a biblioteca
library(plot3D)
#Cria a grade de valores na função
z<-f(M[,"x"], M[,"y"])
#Faz o gráfico
persp3D(z = z, x = M[,"x"], y = M[,"y"],
        expand = 0.3, main = "Função complicada", facets = FALSE, scale = FALSE,
        clab = "Valor da função",zlab="Valor", colkey = list(side = 1, length = 0.5))
O próximo passo para realizar a otimização é escrever a função no formato adequado para o uso do método DEoptim. Para isso, basta trocar o argumento da função que ao invés de serem dois parâmetros, será agora um vetor contendo esses dois argumentos, em outras palavras:
fEval<-function(parms){
  x<-parms[1]
  y<-parms[2]
  res<-(268/((3*x+18)^2+(3*y+7)^2+(18+7)))-(2*(268)/((3*x-21)^2+(3*y+21)^2+(21+21)))-(3*(268)/((3*x+18)^2+(3*y-21)^2+(18+21)))+(4*(268)/((3*x-21)^2+(3*y-7)^2+(21+7))) 
  return(res)
}
Note que só as três primeiras linhas foram alteradas. A seguir, precisamos configurar os nossos clusters para que a função possa ser otimizada em paralelo, a função busca pelo menor valor da função:
#Chama o pacote DEoptim e doSNOW
library(DEoptim)
library(doSNOW)
#Guarda o número de clusters disponíveis
numSlaves<- detectCores()
#Cria os clusters
cl <− makeCluster(numSlaves) 
#Passa todas as bibliotecas que são usadas para o cálculo da função
#clusterEvalQ(cl , library(pracma, ... ))
#Passa todos os objetos necessários para o cômputo da função
#clusterExport(cl , ... ) 
#Registra os clusters
registerDoSNOW (cl)
#Realiza a otimização
outDEoptim <− DEoptim (fEval, c(-15,-15), c(15,15), DEoptim.control(list(trace=FALSE, NP=50, itermax=50, parallelType = 2)))
#Fecha os clusters
stopCluster(cl)
#Guarda os melhores resultados.
fim <- outDEoptim $ optim $ bestmem
fim
Alguns comentários são importantes: as linhas 8 e 10 são utilizadas para se passar bibliotecas e objetos do R que são utilizados no cômputo da função que deve ser otimizada. Nesse caso, como não usamos nenhum pacote especial para o cômputo da função e nem bibliotecas esses comandos estão comentados. Caso deseje-se passar esses valores, basta atribuí-los na forma de lista em cada comando. Na função DEoptim o primeiro argumento é a função a ser minimizada, seguido pelos limites inferiores e superiores para os argumentos. O último argumento é uma lista que controla as condições da função de otimização: NP é o número de populações a serem geradas no algoritmo de Otimização Evolucionária, itermax é o número total de iterações e parallelType = 2 define a forma de paralelismo a ser considerada. Observe ainda que se o valor de NP e/ou itermax e o paralelismo não for adotado, esse processo pode ser muito lento computacionalmente, daí a vantagem de se utilizar o processamento em paralelo que pode reduzir esse tempo computacional significativamente. Mote ainda que o comando trace=FALSE é para suprimir os resultados em cada iteração, e caso deseje-se maximizar a função ao invés de minimizar, basta trabalhar com o valor do negativo da função, ou ainda com seu recíproco.

segunda-feira, 16 de outubro de 2017

Orientação a objetos no R - S4 Classes e Métodos.


Uma das principais atividades e, também uma boa prática na programação usando o R é a utilização da abordagem de programação orientada a objetos. Essa proposta é útil na construção de pacotes no R bem como na organização e estabilidade das bibliotecas criadas. Nesse tutorial iremos trabalhar com a classe S4 de objetos do R.

Suponha que desejamos criar um pacote que nos ajude a calcular a menção de alunos de uma turma. Nesse caso, o sistema de classes S4 pode ser definido criando-se o objeto aluno:

setClass(
  "Aluno",
  representation(nome="character", idade="numeric"),
  prototype(nome=character(0), idade=numeric(0))
)

Esse objeto contêm, inicialmente, dois slots: a idade do aluno e seu nome. O nome é considerado nesse modelo uma variável caracter e a idade uma variável numérica. O comando prototype define quais devem ser os valores iniciais para esses dois atributos.

Entretanto, na universidade nós temos diversos tipos de aluno, eles podem ser: alunos da graduação, pós-graduação, extensão, aluno especial, ouvinte, etc. Suponha que tenhamos dois tipos de alunos apenas:

#Define as subclasses de alunos
setClass("Graduacao",
  representation(nota1="numeric",nota2="numeric",nota3="numeric"),
  contains="Aluno")

setClass("Posgraduacao",
  representation(artigo="character",nota1="numeric",nota2="numeric"),
  contains="Aluno")

Na primeira subclasse donominada Graduação o aluno faz três provas representadas pelos slots: nota1, nota2 e nota3. Já na segunda subclasse o aluno é um aluno de Pós-graduação e ele é avaliado de maneira diferente, tem que preparar um artigo o qual pode ser armazenado no slot artigo e faz duas provas cujas notas são armazenadas em nota1 e nota 2. Note que nesse caso, ambos são Alunos e por isso herdam os atributos de nome e idade.

Além do mais, na criação dessas subclasses não atribuímos os valores iniciais desses atributos por opção. A boa prática, no entanto, recomenda que todos esses atributos sejam definidos os valores de inicialização através do método prototype.

Cada tipo de aluno terá sua nota calculada de maneira diferente: para os alunos de graduação a nota final será a média das três notas nas provas, já os alunos de Pós-graduação a nota será 50% da média das duas provas e 50% se ele fez o artigo e zero se não fez. Nesse caso, usaremos a propriedade de Mutabilidade (Mutability) dos sistemas de programação orientada a objetos. Através dessa propriedade um cômputo diferente será realizado para cada tipo de objeto. O primeiro passo é a criação de um método genérico:

#Cria um método generico denominado calculaMencao
#o qual é aplicado para o tipo de objeto de entrada.
setGeneric(
  "calculaMencao",
  function(object) {
    standardGeneric("calculaMencao")
  }
)

Em seguida, vamos definir um cálculo de função para cada tipo de aluno:

#Para alunos de graduação
setMethod(
  "calculaMencao",
  signature("Graduacao"),
  function(object) {
    #Calcula a média
    media<-(object@nota1+object@nota2+object@nota3)/3
    #Retorna o resultado
    return(media)
  }
)

#Para alunos de pós-graduação
setMethod(
  "calculaMencao",
  signature("Posgraduacao"),
  function(object) {
    #Calcula a média
    media<-(object@nota1+object@nota2)/2
    #Fez o artigo ?
    artigo<-ifelse(object@artigo == "",0,10)
    #Calcula a nota final
    nota<-media*0.5+artigo*0.5
    #Retorna o resultado
    return(nota)
  }
)
Note que em ambas as definições dos métodos incluímos o método genérico calculaMencao. O que muda entre eles é a assinatura (signature) e, consequentemente, a maneira como a função é calculada. Vamos agora criar alguns alunos seguindo a classes S4 formada:
#Aluno de Graduação
joao <- new("Graduacao",
              nome="João Silva",
              idade=22,
              nota1=2.75,
              nota2=8.78,
              nota3=5.72)

#Aluno de Pós-Graduação
maria <- new("Posgraduacao",
            nome="Maria Silva",
            idade=27,
            nota1=8.54,
            nota2=3.21,
            artigo="Como programar em R")

#Aluno de Pós-Graduação
alex <- new("Posgraduacao",
             nome="Alexandre Silva",
             idade=33,
             nota1=2.15,
             nota2=5.48,
             artigo="")
Para não ter que repetir o mesmo código toda vez que desejarmos criar um objeto, podemos construir os constructors:
#Cria o constructor:
Graduacao <- function(nome, idade, nota1, nota2, nota3){
  new("Graduacao", nome=nome, idade=idade, nota1=nota1, nota2=nota2, nota3=nota3)
}

#Cria o constructor:
Posgraduacao <- function(nome, idade, nota1, nota2, artigo){
  new("Graduacao", nome=nome, idade=idade, nota1=nota1, nota2=nota2, artigo=artigo)
}
Podemos calcular a menção deles fazendo:
#Calcula a nota do João
resultadoJoao<-calculaMencao(joao)
resultadoJoao

#Calcula a nota da Maria
resultadoMaria<-calculaMencao(maria)
resultadoMaria

#Calcula a nota do Alex
resultadoAlex<-calculaMencao(alex)
resultadoAlex
Podemos ainda ter objetos que são oriundos de mais de uma classe, por exemplo, considere um aluno que faz tanto Graduação quanto Pós-graduação, podemos criar uma classe específica para ele:
setClass("Maluco",
         contains=c("Graduacao","Posgraduacao"))
A criação dos objetos e cálculo das notas é feito da seguinte forma:
#Registro do Sergio na graduação
sergioGraduacao <- new("Maluco",
                   nome="Sergio Silva",
                   idade=22,
                   nota1=8.36,
                   nota2=9.12,
                   nota3=5.64)

#Registro do Sergio na pós-graduação
sergioPosgraduacao <- new("Posgraduacao",
                      nome="Sergio Silva",
                      idade=22,
                      nota1=3.45,
                      nota2=8.34,
                      artigo="Não sei nada")


#Calcula da nota do Sergio
resultadoSergio1<-calculaMencao(sergioGraduacao)
resultadoSergio1

resultadoSergio2<-calculaMencao(sergioPosgraduacao)
resultadoSergio2
É comum ainda criar para os métodos alguns pontos de validação dos valores de entrada, por exemplo, nenhuma nota pode ser negativa ou maior do que 10. Podemos incluir essas restrições fazendo:
setClass("Graduacao",
         representation(
           nota1="numeric",
           nota2="numeric",
           nota3="numeric"),
         validity = function(object) { 
           if (object@nota1<0 | object@nota1>10 ) 
           { 
             stop("A nota na prova 1 tem que ser válida.") 
           } 
           if (object@nota2<0 | object@nota2>10 ) 
           { 
             stop("A nota na prova 2 tem que ser válida.") 
           } 
           if (object@nota3<0 | object@nota3>10) 
           { 
             stop("A nota na prova 3 tem que ser válida.") 
           } 
           return(TRUE) 
         },
         contains="Aluno"
         )
Nas classes S4 de objetos temos algumas funções especiais:
  • slotNames: fornece o nome dos slots.
  • getSlots: fornece o nome dos slots e seus tipos.
  • getClass: fornece o nome dos slots, seus tipos e das classes antecessoras.




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.

terça-feira, 15 de agosto de 2017

RcppEigen: Operações Matriciais e Rcpp - Parte 3

A biblioteca RcppEigen é uma das mais eficientes para o cômputo numérico de operações algébricas uma vez que utiliza algortimos optimizados para a memória cache e SIMD.

Por exemplo, considere o desafio de estimar uma Matriz de Variâncias e Covariâncias com base em um conjunto de dados. Essa estimação levaria em conta pelo menos $N\times(N-1)/2$ operações (calculando todas as variâncias e covariâncias na matriz triangular uma vez que a matriz é simétrica). Esse procedimento levaria em consideração o uso de loops aninhados (nested loops) os quais podem ser muito ineficientes e demorados quando implementados diretamente em R.

Nesse sentido, podemos usar as operações do pacote RcppEigen também documentadas na página do Eigen C++:

// [[Rcpp::export]]
#include 
// [[Rcpp::depends(RcppEigen)]]
#include 
using namespace Rcpp;

// [[Rcpp::export]]
Eigen::MatrixXd covCalcula(Eigen::MatrixXd mat) {
  //Centra as variáveis com respeito a média
  Eigen::MatrixXd centered = mat.rowwise() - mat.colwise().mean();
  //Calcula a covariância
  Eigen::MatrixXd cov = (centered.adjoint() * centered) / double(mat.rows() - 1);
  return cov;
}

Observe que as operações .rowwise() e .colwise() se referem, respectivamente, aos cômputos para cada linha e para cada coluna. Já o método .adjoint() calcula a Matriz Adjunta.

A invocação desse método no R é realizada da seguinte forma:

#Limpa o Working Directory
#Chama um conjunto de dados
data("ChickWeight")
#Converte todas as variáveis para numeric
dat2 <- data.frame(lapply(ChickWeight, function(x) as.numeric(as.character(x))))
#Transforma para matriz
mat<-as.matrix(dat2)
#Calcula a matriz de covariância
covCalcula(mat)


segunda-feira, 17 de julho de 2017

RcppEigen: Operações Matriciais e Rcpp - Parte 2


Uma das vantagens de se trabalhar com o RcppEigen e Rcpp, além da velocidade, é o fato de que não é necessário implementar todas as funções de interesse em C++. Podemos importar funções já existentes do R para que essas sejam executadas no ambiente Rcpp.

Por exemplo, um problema comum na formação de portfólios é o fato das matrizes de correlações estimadas entre os retornos dos ativos não ser Positivas Definidas. Nesses casos, problemas teóricos e numéricos podem surgir impedindo a obtenção de uma solução ótima para o problema de alocação de portfólios. Higham, N. (2002) propôs uma abordagem numérica para se encontrar a matriz Positiva Definida mais próxima da matriz de ineteresse.

Essa função está disponível no pacote Matrix e é denominada nearPD. Podemos importá-la no ambiente Rcpp da seguinte forma:

// [[Rcpp::export]]
Eigen::MatrixXd nearPDefinite(Eigen::MatrixXd X, int maxit=1e+6, double eigtol = 1e-06, double conv_tol = 1e-07, double posd_tol = 1e-08, bool keepDiagonal=false){
  //Cria o objeto com o Environment.
  Rcpp::Environment G = Rcpp::Environment::global_env();
  //Objeto com as informações do pacote Matrix
  Rcpp::Environment Matrix("package:Matrix");
  //Importa a função de ineteresse
  Function f = Matrix["nearPD"];
  //Obtêm o resultado da função nearPD
  Rcpp::List res = f(X, false, keepDiagonal, true, false, true, false, false,  eigtol, conv_tol, posd_tol, 100, "I", false);
  //Pega o primeiro elemento da lista e transforma para matrix
  Rcpp::NumericMatrix mat = internal::convert_using_rfunction(wrap(res[0]), "as.matrix");
  //Converte de Rcpp::NumericMatrix para Eigen::MatrixXd
  Eigen::MatrixXd Xpd=convertMatrix(mat);
  //Retorna o resultado
  return(Xpd);
}

A função anterior, após compilada, pode ser utilizada no ambiente R de maneira direta pelo nome nearPDefinite. Entretanto, é necessário que a biblioteca Matrix tenha sido habilitada e esteja disponível no environment do R para que seja utilizada. Essa estratégia pode ser para outras funções de outros pacotes também.