perdascarteira_p1.rmd

---
title: "Trabalho sobre Perdas em Carteira"
author: "André Camatta"
date: "20 de agosto de 2019"
output: 
  pdf_document:
    toc: true
    toc_depth: 2
    number_sections: true
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,tidy.opts=list(width.cutoff=60),tidy=TRUE)
```

# Objetivo

Esse trabalho visa modelar a distribuição de perdas de uma carteira. A primeira parte do trabalho consiste em:

* Realizar o ajuste de distribuiçõess para frequência e severidade.

* Considerar diferentes distribuições para severidade - Gama, LogNormal, Pareto, Inverse Gaussian, etc.
Considerar Poisson e Binomial Negativa para frequência.

* Utilizar o método de máxima verossimilhança para a estivamativa de parâmetros e escolher o melhor ajuste via AIC, teste Kolmogov-Smirnov para distribuição contínua ou outros.

# Modelagem da Frequência 

## Carga dos dados de frequência ##
```{r showfreqarq, eval=F}
frequencia.arquivo <- "C:\\Users\\andre\\OneDrive\\Documentos\\
                Atuaria\\TeoriaDoRisco\\perdascarteiras\\1frequencias.txt"
```

```{r runfreqarq, echo=F}
frequencia.arquivo <- "C:\\Users\\andre\\OneDrive\\Documentos\\Atuaria\\TeoriaDoRisco\\perdascarteiras\\1frequencias.txt"
```

```{r freqload}
frequencias.entrada <- read.table(frequencia.arquivo, header=T)
frequencias <- frequencias.entrada[,1]
```

## Estimativa da distribuição de Poisson para a frequência ##

### Maximização da função log de verossimilhança da distribuição de Poisson ###
```{r veropois}
vero.pois <- function(dados,param){
  -sum(log(dpois(dados,param)))
}
estima.num <- nlminb(100,vero.pois,dados=frequencias,lower=0)
lambda<-estima.num$par
lambda
```

### Uso da biblioteca fitdistrplus para máxima verossimilhança da distribuição de Poisson ###
```{r fitpois}
library(fitdistrplus)

summary(frequencias)

fpoisMLE <- fitdist(frequencias, "pois", method="mle")

fpoisMLE

summary(fpoisMLE)
```

É possível observar que o resultado para $\lambda$ é o mesmo.

## Estimativa da distribuição binomial negativa para a frequência ##

### Maximização da função log de verossimilhança da distribuição binomial negativa ###
```{r verobineg}
vero.nbinom <- function(dados,param){
  -sum(log(dnbinom(dados,size=param[1],prob=param[2])))
}
estima.num.nbinom <- nlminb(c(1000,0.5),vero.nbinom,dados=frequencias,lower=c(0,0),upper=c(Inf,1))
estima.num.nbinom$par
```

### Uso da biblioteca fitdistrplus para máxima verossimilhança da distribuição binomial negativa ###
```{r fitnbinom}
fnbinomMLE <- fitdist(frequencias, "nbinom", method="mle")

fnbinomMLE

summary(fnbinomMLE)
```

Nota-se um pequena diferença para o *size* e que foi estimado o parâmetro $\mu$ em vez de *prob*. Porém, é possível obter o valor de *prob* da maneira a seguir para a comparação.

```{r fitnbinom.prob}
prob.nbinomMLE <- fnbinomMLE$estimate["size"]/(fnbinomMLE$estimate["mu"]+fnbinomMLE$estimate["size"])
names(prob.nbinomMLE)<-"prob"
prob.nbinomMLE
```


Nota-se pequena diferença (a partir do quarto dígito significativo) para o valor de *prob*. 



## Análise gráfica da modelagem de frequência ##

A comparação de gráficos dos dados a função de densidade não traria informações relevantes, devido ao pequeno número de pontos de dados (12). Portanto, optou-se por apresentar os dados como um boxplot na horizontal, assim podemos ver a mediana, os quartis 25% e 75% (limites da caixa) e compará-los com as informações da distribuição de Poisson estimada (linhas em vermelho) e da distribuição binomial negativa (linhas em azul). 

```{r graphpois}
frequencias.df<-data.frame(frequencias)
frequencias.df["sinistros"]<-rep("Sinistros",length(frequencias))
library(ggplot2)
frequencias.boxplot<-ggplot(data = frequencias.df) +
  geom_boxplot(aes(x = sinistros, y = frequencias)) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  labs(y="Frequência por observação") +
  coord_flip() 

poisson.iqr<-qpois(0.75,lambda)-qpois(0.25,lambda)
nbinom.size<-fnbinomMLE$estimate["size"]
nbinom.prob<-prob.nbinomMLE
nbinom.mu<-fnbinomMLE$estimate["mu"]
nbinom.p75<-qnbinom(0.75,nbinom.size,nbinom.prob)
nbinom.p25<-qnbinom(0.25,nbinom.size,nbinom.prob)
nbinom.iqr<-nbinom.p75-nbinom.p25

frequencias.boxplot + 
  geom_hline(yintercept = lambda, color="red", size=0.5) + 
  geom_hline(yintercept = qpois(0.75,lambda), linetype="dashed", color="red", size=0.5) +
  geom_hline(yintercept = qpois(0.25,lambda), linetype="dashed", color="red", size=0.5) +
  geom_hline(yintercept = qpois(0.75,lambda), linetype="dashed", color="red", size=0.5) +
  geom_hline(yintercept = qpois(0.25,lambda), linetype="dashed", color="red", size=0.5) +
  geom_hline(yintercept = qpois(0.75,lambda)+1.5*poisson.iqr, linetype="dotdash", color="red", size=0.5) +
  geom_hline(yintercept = qpois(0.25,lambda)-1.5*poisson.iqr, linetype="dotdash", color="red", size=0.5) +
  geom_hline(yintercept = nbinom.mu, color="blue", size=0.5) + 
  geom_hline(yintercept = nbinom.p75, linetype="dashed", color="blue", size=0.5) +
  geom_hline(yintercept = nbinom.p25, linetype="dashed", color="blue", size=0.5) +
  geom_hline(yintercept = nbinom.p75, linetype="dashed", color="blue", size=0.5) +
  geom_hline(yintercept = nbinom.p25, linetype="dashed", color="blue", size=0.5) +
  geom_hline(yintercept = nbinom.p75+1.5*nbinom.iqr, linetype="dotdash", color="blue", size=0.5) +
  geom_hline(yintercept = nbinom.p25-1.5*nbinom.iqr, linetype="dotdash", color="blue", size=0.5)
```

Apresentando, também, os dados de maneira tabular:

```{r freqcountr}
library(Countr)
breaks_ <- c(190, 200, 210, 220, 230, 240, 250, 260, 270)
count_table(count = frequencias, breaks=breaks_, formatChar=TRUE)
```

Observamos que o boxplot dos dados apresentou uma leva assimetria (com mediana um pouco menor que a média) e um *outlier* em seu próprio conjunto, mas que foi abarcado por ambas distribuições. Notamos que a binomial negativa é mais dispersa, aproximando seu desvio padrão dos dados, mas é possível observar que o primeiro (25%) e o terceiro quartil (75%) ficaram muito relaxados em relação aos dados. 

Embora os dados apresentem sobre-dispersão (*overdispersion*) em relação à Poisson, não necessariamente lançar mão da binomial negativa trará ganho significativo.

Abaixo a comparação do desvio padrão dos dados e distribuições estimadas.
```{r varpois}
#Desvio padrão dos dados
sd(frequencias) 
#Desvio padrão da distr. de Poisson
sqrt(lambda) 
#Desvio padrão da distr. binomial negativa
unname(sqrt(nbinom.size*(1-nbinom.prob)/nbinom.prob^2)) 

```

## Utilização de AIC para a escolha da distribuição  ##

A biblioteca *fitdistrplus* já calculou os AICs (Akaike Information Criterion) dessas distribuições. 

```{r AICs}
fpoisMLE$aic #AIC da distribuição de Poisson

fnbinomMLE$aic #AIC da distribuição binomial negativa
```

Da mesma maneira, o AIC pode ser obtido pelos cálculos já realizados pelas funções de log da verossimilhança.
```{r AIC_LOG_VERO}
#AIC da distribuição de Poisson
2*vero.pois(frequencias,estima.num$par) + 2*1 

#AIC da distribuição binomial negativa
2*vero.nbinom(frequencias,estima.num.nbinom$par) + 2*2 
```

Sendo a distribuição de Poisson com menor AIC, ela é a escolhida para representar a frequência de ocorrências de sinistro por período. Ou seja, ~Poisson($\lambda$), onde $\lambda$ é o valor abaixo:

```{r lambda}
lambda
```

# Modelagem da Severidade 

## Carga dos dados de severidade ##
```{r showsinarq, eval=F}
sinistros.arquivo <- "C:\\Users\\andre\\OneDrive\\Documentos\\
                Atuaria\\TeoriaDoRisco\\perdascarteiras\\1sinistros.txt"
```

```{r runsinistrosarq, echo=F}
sinistros.arquivo <- "C:\\Users\\andre\\OneDrive\\Documentos\\Atuaria\\TeoriaDoRisco\\perdascarteiras\\1sinistros.txt"
```

```{r severidadesload}
severidades.entrada <- read.table(sinistros.arquivo, header=T)
severidades <- severidades.entrada[,1]
```

## Estimativas de distribuições para a severidade ##

Abaixo, com o auxílio das bibliotecas fitdistrplus e actuar, estimamos parâmetros para as seguintes distribuições: Log Normal, Gamma, Pareto, Gaussiana Inversa, Weibull e Burr, todas por máxima verossimilhança.


### Estimativa da distribuição Log Normal ###

```{r severidades.lnorm}
slnormMLE <- fitdist(severidades, "lnorm", method="mle")
summary(slnormMLE)
```

### Estimativa da distribuição Gamma ###
```{r severidades.gamma}
sgammaMLE <- fitdist(severidades, "gamma", method="mle", lower=c(0,0))
summary(sgammaMLE)
```

### Estimativa da distribuição de Pareto ###

```{r severidades.pareto}
library(actuar)

sparetoMLE <- fitdist(severidades, "pareto", method="mle")
summary(sparetoMLE)
```

### Estimativa da distribuição Gaussiana Inversa ###

```{r severidades.invgauss}
sinvgaussMLE <- fitdist(severidades, "invgauss", method="mle", start = list(mean = mean(severidades), shape = 1))
summary(sinvgaussMLE)
```

### Estimativa da distribuição de Weibull ###

```{r severidades.weibull}
sweibullMLE <- fitdist(severidades, "weibull", method="mle")
summary(sweibullMLE)
```

### Estimativa da distribuição de Burr ###

```{r severidades.burr}
sburrMLE <- fitdist(severidades, "burr", method="mle", start = list(shape1=0.9,shape2=1,scale=1))
summary(sburrMLE)
```

## Comparação gráfica das distribuições estimadas para a severidade ##


### Comparação das distribuições de frequência ###

Abaixo, a comparação das distribuições de frequências para a severidade que foram estimadas por máxima verossimilhança em comparação com os dados de sinistros que foram fornecidos (representados por histograma). Observa-se que os ajustes da Log Normal, Gaussiana inversa e Burr foram as distribuições que não aparenteram superestimar as ocorrências de sinistros de baixo valor, enquanto o ajuste de distribuição de Burr capturou melhor o maior número de ocorrências entre os valores entre 1000 e 2000 para a severidade. A distribuição de Burr também aparenta o melhor ajuste para valores médios de severidade entre 5000 e 10000. Todas as outras superestimam os valores médios. Não é possível tirar conclusões para valores altos (cauda longa) com esses gráficos. 

```{r severidades.dist.graphs}
par(mfrow=c(2,3))
hist(severidades,nclass=1000,prob=T, xlim=c(0,20000), main="Ajuste Log-Normal")
curve(dlnorm(x,meanlog=slnormMLE$estimate["meanlog"],sdlog=slnormMLE$estimate["sdlog"]),add=T,col="blue",lwd=2)
hist(severidades,nclass=1000,prob=T, xlim=c(0,20000),main="Ajuste Gamma")
curve(dgamma(x,shape=sgammaMLE$estimate["shape"],rate=sgammaMLE$estimate["rate"]),add=T,col="red",lwd=2)
hist(severidades,nclass=1000,prob=T, xlim=c(0,20000), main="Ajuste Pareto")
curve(dpareto(x,shape=sparetoMLE$estimate["shape"],scale=sparetoMLE$estimate["scale"]),add=T,col="green",lwd=2)
hist(severidades,nclass=1000,prob=T, xlim=c(0,20000), main="Ajuste Gaussiana Inversa" )
curve(dinvgauss(x,mean=sinvgaussMLE$estimate["mean"],shape=sinvgaussMLE$estimate["shape"]),add=T,col="purple",lwd=2)
hist(severidades,nclass=1000,prob=T, xlim=c(0,20000), main="Ajuste Weibull" )
curve(dweibull(x,shape=sweibullMLE$estimate["shape"],scale=sweibullMLE$estimate["scale"]),add=T,col="orange",lwd=2)
hist(severidades,nclass=1000,prob=T, xlim=c(0,20000), main="Ajuste Burr" )
curve(dburr(x,shape1=sburrMLE$estimate["shape1"],shape2=sburrMLE$estimate["shape2"],scale=sburrMLE$estimate["scale"]),add=T,col="pink",lwd=2)
```

### Comparação das funções de probabilidade acumuladas ###

Abaixo, comparamos os dados com as três aparentes melhores distribuições (Log Normal, Gaussiana Inversa e Burr) do item anterior em um gráfico de distribuição de probabilidade acumulada e com o eixo x (valor da severidade) em log. Notamos que distribuição de Burr apresenta um ajuste em que não é possível ver diferenças com os dados, enquanto a Gaussiana Inversa e a Log Normal mostram descolamento em certas partes.

```{r severidades.distr.acumuladas}
par(mfrow = c(1, 1))
plot.legend <- c("invgauss", "lognormal", "burr")
cdfcomp(list(sinvgaussMLE, slnormMLE, sburrMLE), legendtext = plot.legend, xlogscale=TRUE)
```


## Escolha da distribuição por qualidade de ajuste ##

Uma função interessante do pacote *fitdistrplus* é a *gofstat*, onde fica mais fácil comparar o AIC, assim como estatísticas de ajuste conhecidas. Para todos os casos, quanto menor, melhor.

```{r severidades.gofstat}
g<-gofstat(list(slnormMLE, sgammaMLE, sparetoMLE, sinvgaussMLE, sweibullMLE, sburrMLE), fitnames = c("lnorm", "gamma", "pareto", "invgauss", "weibull", "burr"))
g
```

Podemos averiguar que a distribuição de Burr apresenta o melhor AIC com larga vantagem, mas com problemas na estatística de Anderson-Darling, aparesentando valor infinito. A distribuição de Burr também apresenta valores infinitos para quantis altos, conforme podemos ver a seguir. 

```{r severidades.quantis}
quantile(severidades,0.99)
quantile(slnormMLE,0.99)
quantile(sgammaMLE,0.99)
quantile(sparetoMLE,0.99)
quantile(sinvgaussMLE,0.99)
quantile(sweibullMLE,0.99)
quantile(sburrMLE,0.99)
```

Isso pode apresentar problemas no futuro, por exemplo, se formos realizar simulações utilizando essa distribuição. O problema pode estar relacionado à precisão computacional dos números. 

### Testes de Kolmogorov-Smirnov para as distribuições estudadas ##

Abaixo, os testes Kolmogov-Smirnov (KS) entre os dados e as distribuições estimadas. POdemos notar que somente a distribuição de Burr não é rejeitada pelos testes. 
```{r kstests}
ks.test(severidades,"plnorm",meanlog=slnormMLE$estimate["meanlog"],sdlog=slnormMLE$estimate["sdlog"])
ks.test(severidades,"pgamma",shape=sgammaMLE$estimate["shape"],rate=sgammaMLE$estimate["rate"])
ks.test(severidades,"ppareto",shape=sparetoMLE$estimate["shape"],scale=sparetoMLE$estimate["scale"])
ks.test(severidades,"pinvgauss",mean=sinvgaussMLE$estimate["mean"],shape=sinvgaussMLE$estimate["shape"])
ks.test(severidades,"pweibull",shape=sweibullMLE$estimate["shape"],scale=sweibullMLE$estimate["scale"])
ks.test(severidades,"pburr", shape1=sburrMLE$estimate["shape1"],shape2=sburrMLE$estimate["shape2"],scale=sburrMLE$estimate["scale"])
```

A função $goftstat$ da biblioteca $fitdistrplus$ também já fez testes KS, que podem ser exibidos como abaixo.
```{r gof_kstest}
g$kstest
```

### Consideração final

A distribuição escolhida para representar a severidade será a Gaussiana Inversa com os parâmetros de média `r sinvgaussMLE$estimate["mean"]` e $shape$ igual a `r sinvgaussMLE$estimate["shape"]`, pois foi a segunda distribuição que apresentou melhor valor para o AIC. 

A escolha pela segunda ocorreu visto que a distribuição Burr, embora tenha apresentado resultados melhores em todos os outros quesitos, apresenta resultados infinitos para determinados quantis, o que pode atrapalhar sobremaneira os resultados das simulações. A distribuição de Burr também poderá ser utilizada desde que seja limitada.

A disitribuição Gaussiana Inversa embora tenha sido rejeitada pelo teste KS, não necessariamente deve ser eliminada para a aproximar os resultados esperados para a severidade.