Amazon Web Service - AWS. Usando RStudio e CUDA na Amazon.

Um dos principais desafios do Análise de Dados atualmente é realizar operações rapidamente em muito conjunto de dados. Nesse sentido, uma vez que o custo de equipamentos para computação em paralelo (Compute Unified Device Architecture) são muito caros para os pesquisadores.

Uma solução é pegar "emprestado" alguns minutos ou horas de servidores que tenham os hardwares necessários. Para isso podemos usar o Amazon Web Service. Nesse post vou mostrar como configurar um servidor Amazon com GPU (Graphics processing unit).

Assumo aqui que você já tem uma conta na Amazon e na Amazon Web Service com cartão de crédito registrado.

O primeiro passo é entrar no Amazon Web Service e fazer o login. Em seguida visitar o AWS Management Console e escolher a opção Compute - EC2:

Escolha agora a opção Launch Instance:

Escolha AWS Market Place e na caixa de Busca digite RStudio:

Nesse exemplo como queremos um servidor com GPU escolhemos RStudio Server with Tensorflow-GPU for AWS:

Aparecerá uma tela com os valores. Se estiver de acordo, clique em Continue:

Em seguida precisamos escolher o tipo da instância. Quanto mais potente for o servidor e seus hardwares mais caro será o minuto de processamento. Os servidores que possuem GPU estão habilitados na página:

Nesse exemplo eu escolhi GPU compute: p2.xlarge e após clicar em Next: Configure Instance Details você deve manter tudo como está, exceto a opção Auto-assign Public IP que deve ser marcada como Enable:

Após essa alteração clique em Review and Launch:

E após revisar, se tudo estiver de acordo clique em Launch:

Por simplicidade vamos proceder sem chave:

Agora seu servidor está sendo inciado. Você deve então clicar no nome da instância (guarde esse nome pois será a senha do RStudio Server):

Quando ele estiver pronto, deve aparecer um sinal verde:

Para acessar o servidor, copie o endereço de IP na parte inferior nomeado como IPv4 Public IP e digite no navegador o IP seguido da porta 8787, algo como: XX.XXX.XX.XX:8787. Copie também o Instance ID na parte inferior pois será a senha do servidor AWS. Caso tudo tenha dado certo a tela de login surgirá no seu navegador:

Utilize como login: rstudio-user e a senha deve ser o Instance ID. Clique em Login e Pronto!! Pode começar a usar...

Agora, quando terminar o uso é preciso terminar a instância para evitar ser cobrado, para isso viste https://console.aws.amazon.com/ec2/. Lá deverá constar todas as suas instâncias, escolha a instância que você quer terminar e clique com o botão direito escolhendo a opção Instace State e Terminate:

Confirme que você realmente deseja terminar a instância:

.

A partir de agora você não será mais cobrado pelo uso do servidor. Caso deseje mais detalhes sobre o servidor visite: https://tensorflow.rstudio.com/tools/cloud_server_gpu.html#amazon-ec2

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:

?

1 2 3 4 5 6

#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:

?

1 2 3 4 5 6

#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:

?

1 2	#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?

?

1 2 3 4 5 6 7 8 9 10

#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:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14

#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:

?

1 2 3 4 5 6

#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.

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)$
onde $\nu$ são os parâmetros variacionais.
3. Encontre a Evidence Lower Bound (ELBO) dada por $L(\nu)=E_{Q}[\log(p(z,x))-\log(q(z|\nu))]$
4. Minimize a quantia $L(\nu)$
5. 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:

?

1 2 3 4 5 6 7 8

#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:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

#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:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

#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)

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:

?

1 2 3 4

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:

?

1 2 3 4 5 6 7 8

#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:

?

1 2 3 4 5 6

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:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

#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.

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:

?

1 2 3 4 5

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:

?

1 2 3 4 5 6 7 8

#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:

?

1 2 3 4 5 6 7 8

#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:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

#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:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

#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:

?

1 2 3 4 5 6 7 8 9

#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:

?

1 2 3 4 5 6 7 8 9 10 11

#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:

?

1 2	setClass("Maluco",          contains=c("Graduacao","Posgraduacao"))

A criação dos objetos e cálculo das notas é feito da seguinte forma:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

#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:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

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.

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:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14

#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:

?

1 2 3 4 5 6 7 8 9 10 11 12

#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:

?

1 2 3 4

#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:

?

1 2 3 4 5 6 7 8 9 10 11

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:

?

1 2 3 4 5 6

#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.

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++:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

// [[Rcpp::export]] #include <rcppeigen.h> // [[Rcpp::depends(RcppEigen)]] #include <cmath> 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; } </cmath></rcppeigen.h>

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:

?

1 2 3 4 5 6 7 8 9

#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)