Análise de Valores Extremos na cidade de Viçosa-Mg: Precipitação e Temperatura

Estudo de Caso: Mês de janeiro(1968-2022)

1 1. Introdução

1.1 Caracterização da Área de Estudo

O município de Viçosa está localizado na mesorregião da Zona da Mata do estado de Minas Gerais, a aproximadamente 230 km a sudeste da capital, Belo Horizonte. Com altitude média de 649 metros, o município ocupa uma área territorial de cerca de 300 km², dos quais 18 km² compreendem o perímetro urbano. Para o ano de 2025, a população foi estimada em 79.517 habitantes.

Figura 1: Vista da praça municipal de Viçosa-MG. FONTE: Google imagens

1.2 Base de Dados e Metodologia

Inicialmente, vamos analisar as precipitações máximas diárias históricas ocorridas exclusivamente no mêses de janeiro, posteriormente a análise será expandida para os demais mêses chuvosos do ano. A série temporal abrange o período de 1968 a 2022, totalizando 55 observações anuais, um tamanho amostral considerado adequado para garantir a robustez estatística das estimativas, como indicado pelos autores, (Cai; Hames, 2011); (Baddour; Kontongomde, 2007); (Roslan; Su Na; Gabda, 2020).

Os dados pluviométricos foram obtidos a partir das estações do Instituto Nacional de Meteorelogia - INMET (Instituto Nacional de Meteorologia, 2024) e do banco de dados do Boletim Meteorológico do Programa de Pós-Graduação em Meteorologia Aplicada da Universidade Federal de Viçosa (UFV) (Universidade Federal de Viçosa, 2023).

Para a estruturação dos dados, adotou-se a metodologia de Máximos de Blocos, onde, para cada mês de janeiro da série histórica (bloco de 31 dias), selecionou-se o evento de maior precipitação diária acumulada.

1.3 Abordagem Estatística

A modelagem probabilística foi realizada utilizando a Distribuição Generalizada de Valores Extremos (GEV). A consistência da série temporal e a adequação do modelo foram verificadas através das seguintes etapas:

1.4 1. Verificação de pressupostos e estrutura de dependência

A modelagem de valores extremos exige uma compreensão rigorosa do comportamento estocástico da série temporal. A aplicação dos testes a seguir visa diagnosticar a presença de tendência e dependência serial, orientando a escolha entre uma abordagem estacionária (parâmetros constantes) ou não-estacionária (parâmetros variáveis no tempo):

• Teste de Cox-Stuart: Proposto por Cox-Stuart Cox; Stuart (1955) o teste é aplicado para verificar a existência de tendências determinísticas. Um asteriscos importante de se destacar é que a detecção de tendência não invalida a análise, mas indica a necessidade de incorporar o tempo como covariável nos parâmetros da distribuição GEV (modelagem não-estacionária) ou proceder com a remoção da tendência para garantir a estacionariedade.

• Teste Augmented Dickey-Fuller (ADF): Proposto por Dickey-Fuller Dickey; Fuller (1979) é utilizado para testar a hipótese de raiz unitária e confirmar o grau de estacionariedade da série. Este diagnóstico é crucial para definir se a modelagem deve ser realizada na série em nível ou em suas diferenças.

• Análise de Autocorrelação (FAC e FACP): As Funções de Autocorrelação (FAC) e Autocorrelação Parcial (FACP) foram examinadas para avaliar a independência e autocorrelação presentes na série, o teste de Ljumg Box foi utilizado para validar ou não a análise visual.

2. Ajuste e Validação: Estimativa dos parâmetros da GEV via máxima verossimilhança e validação da aderência do modelo aos dados empíricos através do teste de Kolmogorov-Smirnov.
3. Projeção de Risco: Estimativa dos níveis de retorno (precipitação máxima esperada) para períodos de recorrência de 5, 10 e 100 anos.

2 Carregamento e Preparação

Nesta etapa, carregamos os pacotes estatísticos necessários e realizamos a leitura dos dados brutos.

library(readxl)
library(purrr)
library(dplyr)
library(tseries)   # ADF Test
library(trend)     # Testes não-paramétricos
library(randtests) # Cox-Stuart
library(goftest)   # Bondade de ajuste (AD test)
library(extRemes)  # Modelagem GEV
library(ismev)     # Diagnósticos GEV
library(urca)      # Teste de Raiz Unitária Avançado
library(DT)        # Tabelas Interativas
library(kableExtra)

# Leitura e Definições
EXTREMOS <- read_excel("EXTREMOS.xlsx")
CHUVA <- EXTREMOS$PRECIP
year <- seq(1968, 2022, 1)

# 1. Encontrando os índices (posições) dos extremos
idx_max <- which.max(CHUVA) # Posição da maior chuva
idx_min <- which.min(CHUVA) # Posição da menor chuva

# 2. Exibindo os Destaques com 'cat'
cat(paste0("---------------------------------------------------\n",
           "RESUMO DOS EXTREMOS (1968-2022):\n",
           "---------------------------------------------------\n",
           "⬆ MAIOR valor observado: ", CHUVA[idx_max], " mm (Ano: ", year[idx_max], ")\n",
           "⬇ MENOR valor observado: ", CHUVA[idx_min], " mm (Ano: ", year[idx_min], ")\n",
           "---------------------------------------------------\n\n"))

## ---------------------------------------------------
## RESUMO DOS EXTREMOS (1968-2022):
## ---------------------------------------------------
## ⬆ MAIOR valor observado: 184.8 mm (Ano: 1986)
## ⬇ MENOR valor observado: 15.1 mm (Ano: 1989)
## ---------------------------------------------------

# 3. Resumo Estatístico Padrão
summary(EXTREMOS$PRECIP)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15.10   34.70   49.00   58.58   80.60  184.80

# 4. Tabela Interativa
datatable(EXTREMOS, options = list(pageLength = 5, scrollX = TRUE), 
          caption = "Tabela 1: Série Histórica de Extremos (Viçosa)")

Historicamente, a série de precipitações máximas diárias para o mês de janeiro apresenta uma média de 58,58 mm. O evento extremo de maior magnitude foi registrado em 1986, com um acumulado de 184,80 mm. Em contrapartida, o ano de 1989 apresentou o menor valor para uma máxima mensal, com apenas 15,10 mm em 24 horas.

3 Análise Exploratória e Testes de Estacionariedade

A análise inicia-se com uma abordagem descritiva dos dados, compreendendo a visualização da série histórica e a análise da distribuição através do histograma de frequências.

plot(CHUVA ~ year, type="b", main="Série Histórica: Precipitação (Jan)", 
     ylab="Precipitação (mm)", xlab="Ano", pch=19, col="#2980b9", lwd=2)
grid()

hist(CHUVA, main="Distribuição de Frequências", 
     xlab="Precipitação (mm)", col="#3498db", border="white", prob=TRUE)
lines(density(CHUVA), col="#c0392b", lwd=2)

3.1 Análise da Folha 2 (Fevereiro)

Verificação breve dos dados referentes à segunda aba da planilha que serão abordados em outra análise.

folha2 <- read_excel("EXTREMOS.xlsx", sheet = 2)
CHUVA2 <- folha2$PRECIP   
year_f2 <- seq(1968, 2022, 1)

par(mfrow=c(1,2))
plot(CHUVA2 ~ year_f2, type="b", main="Série Histórica (Fev)", ylab="Precipitação")
hist(CHUVA2, main="Histograma (Fev)", col="lightgreen")

3.2 Testes Formais

Verificamos a existência de tendência, por meio do teste de Cox; Stuart (1955, cox_stuart_1955), bem como estacionariedade via teste de Dickey; Fuller (1979, 1979) o teste de sequência Wald; Wolfowitz (1940, 1940) foi utilizado para verificar a aleatoriedade da série. O teste de Ljung-Box em conjunto com os gráficos de autocorrelação e autocorrelação parcial foram utilizados para verificar a independência. Como a análise principal foca nos máximos sempre do mesmo mês, não houve necessidade de analisar sazonalidade. As definições formais dos testes, podem também ser encontradas com mais detalhes em Morettin; Toloi (2006, 2006).

Em complemento, a principal referência utilizada para embasamento teórico e dos códigos é referente ao livro Coles (2001).

3.2.1 Definição Teórica: Teste de Tendência de Cox-Stuart

O Teste de Cox-Stuart é um método estatístico não-paramétrico utilizado para verificar a existência de tendências monotônicas (crescentes ou decrescentes) em uma série temporal. Devido à sua natureza baseada em sinais, o teste apresenta baixa sensibilidade a outliers, tornando-o particularmente robusto para análises de valores extremos hidrológicos.

3.2.1.1 Formulação Algébrica

Seja \(X = \{x_1, x_2, \dots, x_n\}\) uma série temporal de observações independentes. O algoritmo do teste consiste nos seguintes passos:

1. Segmentação da Série: A série original é dividida em duas metades iguais. O número de pares (\(c\)) é definido como: \[c = \begin{cases} \frac{n}{2}, & \text{se } n \text{ é par} \\ \frac{n+1}{2}, & \text{se } n \text{ é ímpar} \end{cases}\] (Nota: No caso de \(n\) ímpar, a observação central da série é usualmente descartada para garantir simetria).

2. Cálculo das Diferenças Pares: Cada observação da primeira metade é comparada com sua correspondente na segunda metade da série temporal: \[D_i = x_{i+c} - x_i \quad \text{para } i = 1, \dots, c\]

3. Estatística do Teste (\(S\)): O teste avalia o sinal das diferenças \(D_i\). Define-se a variável indicadora: \[I_i = \begin{cases} 1, & \text{se } D_i > 0 \text{ (crescimento)} \\ 0, & \text{se } D_i < 0 \text{ (decrescimento)} \end{cases}\]

   A estatística \(S\) é dada pela soma dos sinais positivos (ou negativos, dependendo da hipótese lateral): \[S = \sum_{i=1}^{m} I_i\] Onde \(m\) é o número de diferenças não nulas (\(D_i \neq 0\)).

4. Distribuição e Decisão: Sob a Hipótese Nula (\(H_0\)) de ausência de tendência (aleatoriedade), a estatística \(S\) segue uma distribuição Binomial: \[S \sim \text{Binomial}(m, 0.5)\]

   Para amostras com tamanho suficiente (como no presente estudo, onde \(n=55\)), utiliza-se a aproximação pela Distribuição Normal Padronizada para obter o valor \(Z\): \[Z = \frac{|S - \frac{m}{2}| - 0.5}{\sqrt{\frac{m}{4}}}\]

   Rejeita-se \(H_0\) se o p-valor associado a \(Z\) for menor que o nível de significância \(\alpha\) (usualmente 0.05).

3.2.2 Definição Teórica: Teste de Tendência de Mann-Kendall

O Teste de Mann-Kendall é um método estatístico não-paramétrico amplamente utilizado para avaliar a significância de tendências monotônicas (crescimento ou decrescimento) em séries temporais climatológicas e ambientais. Assim como o teste de Cox-Stuart, ele baseia-se em postos (ranks) e diferenças de sinais, conferindo-lhe alta robustez contra outliers e dados faltantes, características essenciais na modelagem de extremos de temperatura.

3.2.2.1 Formulação Algébrica

Seja \(X = \{x_1, x_2, \dots, x_n\}\) uma série temporal de \(n\) observações temporais sequenciais. O algoritmo do teste consiste nos seguintes passos:

1. Comparação de Pares: Cada observação \(x_i\) da série é iterativamente comparada com todas as observações subsequentes \(x_j\), de modo que o tempo \(j > i\).

2. Função Sinal (Sign Function): A diferença entre as observações de cada par é avaliada e padronizada pela função sinal, extraindo-se apenas a direção da mudança: \[sgn(x_j - x_i) = \begin{cases} 1, & \text{se } x_j - x_i > 0 \\ 0, & \text{se } x_j - x_i = 0 \\ -1, & \text{se } x_j - x_i < 0 \end{cases}\]

3. Estatística do Teste (\(S\)): A estatística \(S\) de Mann-Kendall é dada pelo somatório de todos os sinais resultantes dessas comparações pares ao longo de toda a série: \[S = \sum_{i=1}^{n-1} \sum_{j=i+1}^{n} sgn(x_j - x_i)\] Valores muito positivos de \(S\) são indicativos de uma tendência de aquecimento (crescimento), enquanto valores muito negativos sugerem resfriamento (decrescimento).

4. Distribuição e Decisão: Sob a Hipótese Nula (\(H_0\)) de ausência de tendência (os dados são independentes e aleatoriamente ordenados), a estatística \(S\) aproxima-se de uma distribuição Normal com média zero (\(E[S] = 0\)) para amostras com \(n \ge 10\).

   A variância de \(S\), ajustada para a presença de valores repetidos (empates) na série, é dada por: \[Var(S) = \frac{n(n-1)(2n+5) - \sum_{p=1}^{g} t_p(t_p-1)(2t_p+5)}{18}\] Onde \(g\) é o número de grupos com valores empatados e \(t_p\) é a quantidade de dados contidos no \(p\)-ésimo grupo.

   A estatística de teste padronizada \(Z\) é então calculada aplicando-se uma correção de continuidade: \[Z = \begin{cases} \frac{S - 1}{\sqrt{Var(S)}}, & \text{se } S > 0 \\ 0, & \text{se } S = 0 \\ \frac{S + 1}{\sqrt{Var(S)}}, & \text{se } S < 0 \end{cases}\]

   Rejeita-se a hipótese \(H_0\) em favor da existência de uma tendência estatisticamente significativa se o p-valor associado a \(Z\) for menor que o nível de significância \(\alpha\) (usualmente 0.05), ou equivalentemente, se \(|Z| > Z_{1-\alpha/2}\).

3.2.3 Comparação entre os Testes de Mann-Kendall e Cox-Stuart

Embora os testes de Mann-Kendall e Cox-Stuart sejam métodos não-paramétricos destinados à identificação de tendências monotônicas em séries temporais, eles diferem substancialmente na forma como processam a informação temporal, o que impacta diretamente o poder estatístico de cada abordagem.

O teste de Cox-Stuart adota uma estratégia de segmentação, dividindo a série pela metade e comparando apenas os pares de observações correspondentes (\(n/2\) comparações). Esta simplificação descarta as variações que ocorrem entre os pontos adjacentes dentro da mesma metade da série, resultando em uma considerável perda de informação e, consequentemente, conferindo ao teste um baixo poder estatístico para detectar tendências mais sutis.

Em contrapartida, o teste de Mann-Kendall realiza uma avaliação exaustiva da série. O algoritmo compara cada observação com todas as observações subsequentes na linha do tempo, totalizando \(n(n-1)/2\) comparações pareadas. Ao utilizar a totalidade da estrutura topológica dos dados, o Mann-Kendall apresenta um poder estatístico significativamente superior. Por esta razão, é amplamente reconhecido na literatura hidrológica e climatológica como o método mais robusto e confiável para a detecção de tendências, sendo o critério principal para justificar a adoção de parâmetros variantes no tempo (não-estacionariedade) na modelagem de valores extremos.

3.2.4 Definição Teórica: Teste Augmented Dickey-Fuller (ADF)

Para a verificação de estacionariedade (condição em que a média e a variância da série são constantes ao longo do tempo), utilizou-se o Teste Augmented Dickey-Fuller (ADF). Este método testa a presença de uma raiz unitária na série temporal, o que indicaria um comportamento não estacionário.

3.2.4.1 Formulação Algébrica

O teste ADF baseia-se na estimação da seguinte equação de regressão para a série temporal \(y_t\):

\[\Delta y_t = \alpha + \beta t + \gamma y_{t-1} + \sum_{i=1}^{p} \delta_i \Delta y_{t-i} + \epsilon_t\]

Onde:

• \(\Delta y_t\): Diferença de primeira ordem (\(y_t - y_{t-1}\)).
• \(\alpha\): Constante (intercepto ou deriva).
• \(\beta t\): Coeficiente de tendência determinística temporal.
• \(\gamma\): O parâmetro de interesse (coeficiente da raiz unitária, onde \(\gamma = \rho - 1\)).
• \(p\): Número de defasagens (lags) incluídas para corrigir a correlação serial dos erros (parte “Augmented” do teste).
• \(\epsilon_t\): Termo de erro (ruído branco).

3.2.4.2 Hipóteses e Estatística do Teste

O foco do teste recai sobre o coeficiente \(\gamma\). As hipóteses são definidas como:

1. Hipótese Nula (\(H_0: \gamma = 0\)): Existe uma raiz unitária. A série não é estacionária (possui comportamento de passeio aleatório).
2. Hipótese Alternativa (\(H_1: \gamma < 0\)): Não existe raiz unitária. A série é estacionária.

A estatística de teste (\(DF_{\tau}\)) é calculada de forma análoga a um teste-t convencional:

\[DF_{\tau} = \frac{\hat{\gamma}}{SE(\hat{\gamma})}\]

No entanto, sob a hipótese nula de raiz unitária, essa estatística não segue uma distribuição t-Student padrão, mas sim a distribuição assintótica de Dickey-Fuller. Se a estatística calculada for menor (mais negativa) que o valor crítico da tabela de Dickey-Fuller (ou se o p-valor < 0.05), rejeita-se \(H_0\), concluindo-se pela estacionariedade da série.

3.3 Teste de Sinal (Sign Test)

O Teste de Sinal é um dos procedimentos não paramétricos mais antigos e fundamentais para a inferência estatística. Ele é utilizado principalmente para testar hipóteses sobre a mediana populacional (\(\theta\)) de uma distribuição contínua, servindo como uma alternativa robusta ao Teste-t quando a suposição de normalidade dos dados não é atendida ou quando a distribuição apresenta assimetria severa e outliers.

3.3.1 Formulação das Hipóteses

Seja \(X_1, X_2, \dots, X_n\) uma amostra aleatória de tamanho \(n\) proveniente de uma população contínua com mediana \(\theta\). O objetivo é testar se a mediana populacional é igual a um valor hipotético \(\theta_0\). As hipóteses podem ser formuladas da seguinte maneira:

• Bilateral: \(H_0: \theta = \theta_0\) contra \(H_1: \theta \neq \theta_0\)
• Unilateral à Direita: \(H_0: \theta \leq \theta_0\) contra \(H_1: \theta > \theta_0\)
• Unilateral à Esquerda: \(H_0: \theta \geq \theta_0\) contra \(H_1: \theta < \theta_0\)

3.3.2 Estatística do Teste

A construção da estatística do teste baseia-se nas diferenças entre as observações e a mediana hipotética: \(D_i = X_i - \theta_0\).

Sob a hipótese nula (\(H_0\)), a probabilidade de uma observação ser maior que a mediana é igual à probabilidade de ser menor, ou seja, \(P(X_i > \theta_0) = P(X_i < \theta_0) = 0,5\).

Definimos a variável indicadora \(\psi_i\) tal que:

\[ \psi_i = \begin{cases} 1, & \text{se } X_i > \theta_0 \quad (\text{sinal positivo}) \\ 0, & \text{se } X_i < \theta_0 \quad (\text{sinal negativo}) \end{cases} \]

Nota: Observações onde \(X_i = \theta_0\) são, por convenção, descartadas da análise, reduzindo o tamanho amostral \(n\).

A estatística do teste, denotada por \(S\), é o número total de diferenças positivas (sucessos):

\[ S = \sum_{i=1}^{n} \psi_i \]

3.3.3 Distribuição Amostral e Regra de Decisão

Como \(S\) é uma soma de variáveis de Bernoulli independentes, sob \(H_0\), a estatística \(S\) segue uma Distribuição Binomial com parâmetros \(n\) e \(p=0,5\):

\[ S \sim \text{Bin}(n; 0,5) \]

A decisão estatística depende do tamanho da amostra:

1. Para pequenas amostras (\(n \leq 20\)): Utiliza-se a probabilidade exata da distribuição Binomial. O \(p\)-valor é calculado somando-se as probabilidades das caudas correspondentes ao valor observado de \(S\).

\[ P(S = k) = \binom{n}{k} (0,5)^n \]

2. Aproximação para grandes amostras (\(n > 20\)): Pelo Teorema Central do Limite, a distribuição de \(S\) aproxima-se de uma Distribuição Normal. Para melhorar a precisão da aproximação de uma variável discreta por uma contínua, aplica-se a correção de continuidade.

A estatística padronizada \(Z\) é dada por:

\[ Z = \frac{S - E[S]}{\sqrt{Var(S)}} = \frac{(S \pm 0,5) - \frac{n}{2}}{\frac{\sqrt{n}}{2}} \]

Onde: * \(E[S] = np = \frac{n}{2}\) (Esperança sob \(H_0\)) * \(Var(S) = np(1-p) = \frac{n}{4}\) (Variância sob \(H_0\)) * O termo \(\pm 0,5\) refere-se à correção de continuidade (subtrai-se 0,5 se \(S > n/2\) e soma-se 0,5 se \(S < n/2\)).

Se \(|Z| > Z_{\alpha/2}\) (para teste bilateral), rejeita-se a hipótese nula ao nível de significância \(\alpha\).

---

Relação com o Teste de Cox-Stuart: O teste de Cox-Stuart, utilizado na verificação de tendência desta análise, é uma aplicação direta do Teste de Sinal. Nele, a série temporal é dividida em duas metades e o Teste de Sinal é aplicado às diferenças pareadas (\(X_{i+c} - X_i\)), testando se a “mediana das diferenças” difere de zero, o que evidenciaria uma tendência determinística.

3.3.3.1 O Teste de Ljung-Box

O teste de Ljung-Box (também conhecido como teste Q portmanteau) é utilizado para verificar se as autocorrelações de uma série temporal são, em conjunto, significativamente diferentes de zero. Em diagnósticos de modelos, ele é fundamental para confirmar se os resíduos se comportam como ruído branco.

3.3.3.2 Formulação Algébrica

A estatística de teste \(Q\) é calculada com base nas autocorrelações amostrais até uma defasagem \(h\):

\[Q = n(n+2) \sum_{k=1}^{h} \frac{\hat{\rho}_k^2}{n-k}\]

Onde:

• \(n\): Tamanho da amostra (número de observações).
• \(\hat{\rho}_k\): Autocorrelação amostral na defasagem \(k\).
• \(h\): Número de defasagens (lags) testadas em conjunto.
• \(k\): O índice de cada defasagem individual.

[Image of Chi-square distribution for different degrees of freedom]

3.3.3.3 Hipóteses e Estatística do Teste

O teste avalia se a estrutura de dependência serial foi completamente capturada ou se ainda existe informação sistemática na série. As hipóteses são:

1. Hipótese Nula (\(H_0\)): \(\rho_1 = \rho_2 = \dots = \rho_h = 0\). Os dados são independentes (não há autocorrelação; a série é ruído branco).
2. Hipótese Alternativa (\(H_1\)): Ao menos um \(\rho_k \neq 0\) para \(k \in \{1, \dots, h\}\). Os dados apresentam dependência serial (não são ruído branco).

Sob a hipótese nula, a estatística \(Q\) segue uma distribuição Qui-quadrado (\(\chi^2\)). A definição dos graus de liberdade (\(gl\)) depende da aplicação:

• Para séries originais: \(gl = h\).
• Para resíduos de modelos (ex: ARIMA): \(gl = h - m\), onde \(m\) é a soma dos parâmetros estimados (ex: \(p + q\)).

Se a estatística \(Q\) calculada for superior ao valor crítico da distribuição \(\chi^2\) para o nível de significância adotado (ou se o p-valor < 0,05), rejeita-se \(H_0\), indicando que o modelo não capturou toda a estrutura de correlação da série.

# 1. Análise Visual de Dependência (ACF e PACF)
# Configura a área de plotagem para 1 linha e 2 colunas
par(mfrow = c(1, 2)) 

acf(CHUVA, main = "Função de Autocorrelação (FAC)", ylab = "ACF")
pacf(CHUVA, main = "Autocorrelação Parcial (FACP)", ylab = "PACF")

Gráficos de Autocorrelação (ACF) e Autocorrelação Parcial (PACF).

# H0: Os dados são independentes, não há autocorrelação.
Box.test(CHUVA, lag = 10, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  CHUVA
## X-squared = 6.6114, df = 10, p-value = 0.7616

par(mfrow = c(1, 1)) 

# 2. Teste de Estacionariedade (Augmented Dickey-Fuller)
# type = "trend" considera intercepto e tendência na equação do teste
teste_adf <- ur.df(CHUVA, type = "trend", selectlags = "AIC")
summary(teste_adf)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -46.535 -23.250  -7.984  19.973 125.164 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  61.7103    14.7502   4.184 0.000118 ***
## z.lag.1      -1.1769     0.1999  -5.887 3.49e-07 ***
## tt            0.2421     0.3209   0.755 0.454088    
## z.diff.lag    0.1654     0.1413   1.171 0.247277    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.54 on 49 degrees of freedom
## Multiple R-squared:  0.5175, Adjusted R-squared:  0.488 
## F-statistic: 17.52 on 3 and 49 DF,  p-value: 7.303e-08
## 
## 
## Value of test-statistic is: -5.8871 11.5592 17.3385 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2  6.50  4.88  4.16
## phi3  8.73  6.49  5.47

# NOTA DE INTERPRETAÇÃO (ADF):
# Compare a estatística de teste (tau3) com os valores críticos (Critical Values).
# H0: Existe Raiz Unitária (Não estacionária).
# Se t-statistic < Critical Value (5%), rejeita-se H0 (Série é Estacionária).

# 3. Teste de Tendência (Cox-Stuart)
# H0: Não existe tendência monotônica.
cox.stuart.test(CHUVA)

## 
##  Cox Stuart test
## 
## data:  CHUVA
## statistic = 16, n = 27, p-value = 0.4421
## alternative hypothesis: non randomness

# 4. Teste de Aleatoriedade (Runs Test)
# H0: A série é aleatória.
runs.test(CHUVA)

## 
##  Runs Test
## 
## data:  CHUVA
## statistic = -0.27477, runs = 27, n1 = 27, n2 = 27, n = 54, p-value =
## 0.7835
## alternative hypothesis: nonrandomness

---

4 Fundamentação Teórica e Implementação Manual

A Distribuição Generalizada de Valores Extremos (GEV), introduzida por Jenkinson (1955), fornece uma estrutura analítica robusta para modelar máximos de blocos. A função de distribuição acumulada é dada por:

\[G(z) = \exp\left\{ -\left[ 1 + \xi \left( \frac{z - \mu}{\sigma} \right) \right]^{-1/\xi} \right\}\]

4.1 Definição da Log-Likelihood Manual

# Variável global para as funções manuais
dataset <- CHUVA

gev.loglik <- function(theta){
  mu <- theta[1]
  sigma <- theta[2]
  xi <- theta[3]
  
  m <- min((1+(xi*(dataset-mu)/sigma)))
  if(m < 0.00001) return(as.double(1e6))
  if(sigma < 0.00001) return(as.double(1e6))
  
  if(xi == 0){ # Caso Gumbel
    loglik <- (-length(dataset)*log(sigma)-sum((dataset-mu)/sigma)
               -sum(exp(-((dataset-mu)/sigma))))
  } else { # Caso GEV
    loglik <- (-length(dataset)*log(sigma)
               -(1/xi+1)*sum(log(1+(xi*(dataset-mu)/sigma)))
               -sum((1+(xi*(dataset-mu)/sigma))^(-1/xi)))
  }
  return(-loglik)
}

4.2 Otimização e Estimativa de Parâmetros

# Chute inicial
theta_start <- c(mean(dataset), sd(dataset), 0.1)

# Otimização via Newton-type algorithm
A <- nlm(gev.loglik, theta_start, hessian=TRUE)
print(A$estimate) # Parâmetros Estimados (Mu, Sigma, Xi)

## [1] 40.7879940 21.2843617  0.2322431

# Matriz de Variância-Covariância
varcovar <- solve(A$hessian)
cat("Erros Padrão dos Parâmetros:\n")

## Erros Padrão dos Parâmetros:

print(sqrt(diag(varcovar)))

## [1] 3.383208 2.775215 0.137969

4.3 Diagnóstico Gráfico Manual

# Funções auxiliares CDF e Inversa
GEV.DF <- function(data, mu, sigma, xi){
  if(xi==0) { GEV=exp(-exp(-((data-mu)/sigma))) }
  else { GEV=exp(-(1+xi*((data-mu)/sigma))^(-1/xi)) }
  return(GEV)
}

GEV.INV <- function(data, mu, sigma, xi){
  if(xi==0) { INV=mu-sigma*log(-log(1-data)) } # Ajuste teórico para p
  else { INV=mu+(sigma/xi)*(((-log(data))^(-xi))-1) }
  return(INV)
}

# Preparação dos dados ordenados
ordered <- sort(CHUVA)
n <- length(ordered)
empirical <- (1:n)/(n+1)

# Geração dos valores do modelo
model_probs <- numeric(n)
model_quantile <- numeric(n)

for(i in 1:n){
  model_probs[i] <- GEV.DF(ordered[i], A$est[1], A$est[2], A$est[3])
  model_quantile[i] <- GEV.INV(empirical[i], A$est[1], A$est[2], A$est[3])
}

# Plots
par(mfrow=c(1,2))
plot(model_probs ~ empirical, main='PP-Plot (Manual)', xlab="Empírico", ylab="Modelo")
abline(0,1, col="red")
plot(model_quantile ~ ordered, main='QQ-Plot (Manual)', xlab="Dados", ylab="Modelo")
abline(0,1, col="red")

4.4 Cálculo do Erro Padrão do Nível de Retorno (Método Delta)

Cálculo manual da variância para o nível de retorno de 10 anos.

y10 <- -log(1-(1/10))
del <- matrix(ncol=1, nrow=3)

# Derivadas parciais
del[1,1] <- 1
del[2,1] <- -((A$est[3])^(-1))*(1-(y10^(-A$est[3])))
del[3,1] <- (((A$est[2])*((A$est[3])^(-2))*(1-((y10)^(-A$est[3]))))
             -((A$est[2])*((A$est[3])^(-1))*((y10)^(-(A$est[3])))*log(y10)))

del.transpose <- t(del)
se_return <- sqrt(del.transpose %*% varcovar %*% del)

cat("Erro Padrão (10 anos):", se_return)

## Erro Padrão (10 anos): 12.76559

---

5 Modelagem Otimizada e Validação

Para fins de robustez e otimização, diagnóstico e cálculo de intervalos de confiança via Profile Likelihood, é útil a utilização de alguns pacotes já existentes, como: extRemes e ismev.

5.1 Diagnóstico via ismev

O pacote ismev fornece um painel completo de diagnósticos visuais que permite uma análise visual da adequação do modelo, bem como intervalos de confiança dos níveis de retorno.

fit_ismev <- gev.fit(CHUVA)

## $conv
## [1] 0
## 
## $nllh
## [1] 262.355
## 
## $mle
## [1] 40.7846235 21.2816111  0.2321556
## 
## $se
## [1] 3.3823308 2.7733506 0.1378968

gev.diag(fit_ismev)

5.2 Ajuste e Validação via extRemes

# Ajuste do Modelo
fit1 <- fevd(CHUVA, units = "mm")
print(fit1)

## 
## fevd(x = CHUVA, units = "mm")
## 
## [1] "Estimation Method used: MLE"
## 
## 
##  Negative Log-Likelihood Value:  262.355 
## 
## 
##  Estimated parameters:
##   location      scale      shape 
## 40.7964849 21.2901341  0.2320783 
## 
##  Standard Error Estimates:
##  location     scale     shape 
## 3.3837941 2.7749166 0.1379423 
## 
##  Estimated parameter covariance matrix.
##            location       scale       shape
## location 11.4500626  5.84873820 -0.18164368
## scale     5.8487382  7.70016205 -0.07692724
## shape    -0.1816437 -0.07692724  0.01902808
## 
##  AIC = 530.7099 
## 
##  BIC = 536.7319

5.2.1 Testes de Adequação

Foram aplicados os testes de Kolmogorov-Smirnov e Anderson-Darling para confirmar estatisticamente a aderência do modelo.

• Hipótese Nula (\(H_0\)): O modelo ajustado descreve bem os dados.
• Critério: P-valor > 0.05 indica bom ajuste.

# Recuperando parâmetros estimados
params <- fit1$results$par
loc_est <- params["location"]
scale_est <- params["scale"]
shape_est <- params["shape"]

# Teste KS
ks_result <- ks.test(CHUVA, "pevd", loc=loc_est, scale=scale_est, shape=shape_est, type="GEV")
print(ks_result)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  CHUVA
## D = 0.095798, p-value = 0.6938
## alternative hypothesis: two-sided

# Teste AD (Anderson-Darling)
ad_result <- ad.test(CHUVA, "pevd", loc=loc_est, scale=scale_est, shape=shape_est, type="GEV")
print(ad_result)

## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: distribution 'pevd'
##  with parameters loc = 40.7964848566098, scale = 21.2901340682425, shape
##  = 0.232078322647953, type = GEV
##  Parameters assumed to be fixed
## 
## data:  CHUVA
## An = 0.41043, p-value = 0.8377

5.3 Estimativa de Níveis de Retorno

Por fim, foram calculados os níveis de retorno para 5, 10 e 20 anos com intervalos de confiança baseados na Verossimilhança perfilhada (Profile Likelihood), que são mais precisos para extremos.

par(mfrow=c(1,3))
# 5 Anos
ci(fit1, method = "proflik", xrange = c(65, 100), verbose = TRUE, return.period = 5)

## 
##  Preparing to calculate  95 % CI for  5-year return level 
## 
##  Model is   fixed 
## 
##  Using Profile Likelihood Method.
## 
##  Calculating profile likelihood.  This may take a few moments.

## 
##  Profile likelihood has been calculated.  Now, trying to find where it crosses the critical value =  -264.2757

## fevd(x = CHUVA, units = "mm")
## 
## [1] "Profile Likelihood"
## 
## [1] "5-year return level: 78.994"
## 
## [1] "95% Confidence Interval: (66.5829, 98.2412)"

# 10 Anos
ci(fit1, method = "proflik", xrange = c(80, 160), verbose = TRUE, return.period = 10)

## 
##  Preparing to calculate  95 % CI for  10-year return level 
## 
##  Model is   fixed 
## 
##  Using Profile Likelihood Method.
## 
##  Calculating profile likelihood.  This may take a few moments.

## 
##  Profile likelihood has been calculated.  Now, trying to find where it crosses the critical value =  -264.2757

## fevd(x = CHUVA, units = "mm")
## 
## [1] "Profile Likelihood"
## 
## [1] "10-year return level: 103.713"
## 
## [1] "95% Confidence Interval: (84.9722, 143.2425)"

# 100 Anos
ci(fit1, method = "proflik", xrange = c(120, 470), verbose = TRUE, return.period = 100)

## 
##  Preparing to calculate  95 % CI for  100-year return level 
## 
##  Model is   fixed 
## 
##  Using Profile Likelihood Method.
## 
##  Calculating profile likelihood.  This may take a few moments.

## 
##  Profile likelihood has been calculated.  Now, trying to find where it crosses the critical value =  -264.2757

## fevd(x = CHUVA, units = "mm")
## 
## [1] "Profile Likelihood"
## 
## [1] "100-year return level: 215.865"
## 
## [1] "95% Confidence Interval: (144.7157, 451.5789)"

As estimativas indicam que um evento de precipitação de 78,94 mm possui um tempo de recorrência de 5 anos para o mês de janeiro. Para períodos de retorno mais longos, de 10 e 100 anos, os volumes estimados são de 103 mm e 215,865 mm, respectivamente. Estatisticamente, isso implica que tais magnitudes devem ser igualadas ou superadas, em média, uma vez a cada intervalo de tempo considerado.

5.3.0.1 Generalizando para os demais mêses

# ==============================================================================
# FUNÇÃO GERAL DE ANÁLISE CLIMÁTICA 
# ==============================================================================
analisar_clima <- function(dados, ano_vec, nome_mes, cor_tema) {

  cat(paste0("\n### Análise Descritiva: ", nome_mes, "\n"))
  
  # --- 1. Descritiva e Gráficos ---
  idx_max <- which.max(dados)
  cat(paste0("**Evento Máximo:** ", round(dados[idx_max], 2), " mm em ", ano_vec[idx_max], "\n\n"))
  
  par(mfrow=c(2,2))
  df_plot <- data.frame(Ano = ano_vec, Chuva = dados); df_plot <- na.omit(df_plot)
  
  plot(Chuva ~ Ano, data = df_plot, type="b", main=paste("Série -", nome_mes), 
       ylab="Precipitação (mm)", xlab="Ano", pch=19, col=cor_tema, lwd=1.5)
  grid()
  hist(df_plot$Chuva, main=paste("Histograma -", nome_mes), 
       xlab="mm", col=cor_tema, border="white", prob=TRUE)
  lines(density(df_plot$Chuva), col="black", lwd=2)
  acf(df_plot$Chuva, main="ACF", col=cor_tema, lag.max = 15)
  pacf(df_plot$Chuva, main="PACF", col=cor_tema, lag.max = 15)
  par(mfrow=c(1,1)) 
  
  # --- Função de Badge ---
  criar_badge <- function(texto, cor) {
    paste0("<span style='display: block; padding: 2px 8px; border-radius: 4px; ",
           "color: white; background-color: ", cor, "; font-weight: bold; text-align: center;'>", 
           texto, "</span>")
  }
  
  # --- 2. Bateria de Testes de Diagnóstico ---
  cat("\n#### Diagnóstico da Série Temporal\n")
  
  # A) Cox-Stuart (Tendência Sinais)
  cs_test <- cox.stuart.test(dados)
  p_cs <- cs_test$p.value
  
  # B) Mann-Kendall (Tendência Rankeada - Mais Potente)
  mk_test <- trend::mk.test(dados)
  p_mk <- mk_test$p.value
  status_mk <- if(p_mk >= 0.05) criar_badge("Sem Tendência", "#27ae60") else criar_badge("Tendência (MK)", "#e67e22")
  
  # C) Teste de Pettitt (Quebra Estrutural / Salto)
  pt_test <- trend::pettitt.test(dados)
  p_pt <- pt_test$p.value
  # Identifica o ano da provável quebra se houver
  ano_quebra <- if(p_pt < 0.05) paste0(" (Em ", ano_vec[pt_test$estimate], ")") else ""
  status_pt <- if(p_pt >= 0.05) criar_badge("Homogênea", "#27ae60") else criar_badge(paste0("Quebra", ano_quebra), "#c0392b")

  # D) ADF (Estacionariedade)
  adf_check <- tseries::adf.test(dados, alternative = "stationary")
  p_adf <- adf_check$p.value
  status_adf <- if(p_adf < 0.05) criar_badge("Estacionária", "#27ae60") else criar_badge("Não Estacionária", "#c0392b")
  
  # E) Ljung-Box (Independência)
  lb_check <- Box.test(dados, type = "Ljung-Box", lag = 10)
  p_lb <- lb_check$p.value
  status_lb <- if(p_lb >= 0.05) criar_badge("Independente", "#27ae60") else criar_badge("Dependência", "#c0392b")
  
  # Tabela Diagnóstico Atualizada
  df_pre <- data.frame(
    Teste = c("Cox-Stuart (Tendência)", "Mann-Kendall (Tendência)", "Pettitt (Quebra)", "ADF (Estacionariedade)", "Ljung-Box (Memória)"),
    Estatistica = c(round(cs_test$statistic, 4), round(mk_test$statistic, 4), round(pt_test$statistic, 4), round(adf_check$statistic, 4), round(lb_check$statistic, 4)),
    P_Valor = c(format.pval(p_cs, digits=4, eps=0.001), 
                format.pval(p_mk, digits=4, eps=0.001),
                format.pval(p_pt, digits=4, eps=0.001),
                format.pval(p_adf, digits=4, eps=0.001),
                format.pval(p_lb, digits=4, eps=0.001)),
    Diagnostico = c(
      if(p_cs >= 0.05) criar_badge("Sem Tendência", "#27ae60") else criar_badge("Tendência", "#e67e22"),
      status_mk,
      status_pt,
      status_adf,
      status_lb
    )
  )
  
  tabela_html_pre <- df_pre %>%
    kbl(caption = paste("Tabela A: Diagnóstico Avançado -", nome_mes), align = "c", escape = FALSE, row.names = FALSE) %>%
    kable_styling(bootstrap_options = c("hover", "condensed"), full_width = F, position = "left") %>%
    column_spec(1, bold = T) %>%
    column_spec(4, width = "180px")
    
  cat(tabela_html_pre)
  cat("\n\n")
  
  # --- 3. Ajuste GEV e Validação ---
  cat("\n#### Modelagem GEV e Validação\n")
  fit <- fevd(dados, type = "GEV", units = "mm")
  params <- fit$results$par
  loc_e <- params["location"]; sc_e <- params["scale"]; sh_e <- params["shape"]
  
  # Testes de Bondade de Ajuste
  ks_res <- ks.test(dados, "pevd", loc=loc_e, scale=sc_e, shape=sh_e, type="GEV")
  ad_res <- goftest::ad.test(dados, "pevd", loc=loc_e, scale=sc_e, shape=sh_e, type="GEV")
  
  df_gof <- data.frame(
    Teste = c("Kolmogorov-Smirnov", "Anderson-Darling"),
    Estatistica = c(round(ks_res$statistic, 4), round(ad_res$statistic, 4)),
    P_Valor = c(format.pval(ks_res$p.value, digits=4, eps=0.001), format.pval(ad_res$p.value, digits=4, eps=0.001)),
    Conclusao = c(
      if(ks_res$p.value > 0.05) criar_badge("Bom Ajuste", "#27ae60") else criar_badge("Mau Ajuste", "#c0392b"),
      if(ad_res$p.value > 0.05) criar_badge("Bom Ajuste", "#27ae60") else criar_badge("Mau Ajuste", "#c0392b")
    )
  )
  
  tabela_html_gof <- df_gof %>%
    kbl(caption = paste("Tabela B: Validação do Modelo GEV -", nome_mes), align = "c", escape = FALSE, row.names = FALSE) %>%
    kable_styling(bootstrap_options = c("hover", "condensed"), full_width = F, position = "left") %>%
    column_spec(1, bold = T) %>%
    column_spec(4, width = "140px")
    
  cat(tabela_html_gof)
  cat("\n\n")
  
  # --- 4. Gráficos e Retorno ---
  par(mfrow=c(1,2))
  plot(fit, type = "density", main = "Densidade Ajustada")
  plot(fit, type = "qq", main = "QQ-Plot")
  par(mfrow=c(1,1))
  
  cat("\n#### Estimativas de Risco (Nível de Retorno)\n")
  par(mfrow=c(1,3))
  plot_ci_seguro <- function(modelo, periodo, titulo) {
    tryCatch({
      ci(modelo, return.period = periodo, method = "proflik", verbose = FALSE, main = titulo)
    }, error = function(e) {
      ci(modelo, return.period = periodo, method = "normal", verbose = FALSE, main = paste(titulo, "(Normal Appx)"))
    })
  }
  plot_ci_seguro(fit, 5, "5 Anos"); plot_ci_seguro(fit, 10, "10 Anos"); plot_ci_seguro(fit, 50, "50 Anos")
  par(mfrow=c(1,1))
  
  return(fit)
}

# ==============================================================================
# LEITURA DOS DADOS 
# ==============================================================================

# Definição do vetor de anos
anos <- seq(1968, 2022, 1)

# Lendo Jan (Sheet 1)
df_jan <- read_excel("EXTREMOS.xlsx", sheet = 1)
dados_jan <- df_jan$PRECIP

# Lendo Fev (Sheet 2)
df_fev <- read_excel("EXTREMOS.xlsx", sheet = 2)
dados_fev <- df_fev$PRECIP

# Lendi Março (Sheet 3)
df_mar<- read_excel("EXTREMOS.xlsx", sheet = 3)
dados_mar<- df_mar$PRECIP

# Lendo Dez (Sheet 12)
df_dez <- read_excel("EXTREMOS.xlsx", sheet = 12)
dados_dez <- df_dez$PRECIP

#Lendo Nov (Sheet 11)
df_nov <- read_excel("EXTREMOS.xlsx", sheet = 11)
dados_nov<- df_nov$PRECIP

# ==============================================================================
# INVESTIGAÇÃO DETALHADA: ESTACIONARIEDADE DE FEVEREIRO
# ==============================================================================
library(readxl)
library(urca)
library(ggplot2)

# 1. Leitura isolada dos dados de Fevereiro (para garantir que estamos testando o certo)
df_fev_investigacao <- read_excel("EXTREMOS.xlsx", sheet = 2)
dados_fev_check <- na.omit(df_fev_investigacao$PRECIP)

# 2. Visualização Rápida

plot(dados_fev_check, type="l", col="#c0392b", lwd=2, 
     main="Investigação Visual: Série de Fevereiro", ylab="Precipitação (mm)")
abline(reg=lm(dados_fev_check~time(dados_fev_check)), col="blue", lty=2) # Linha de tendência
legend("topleft", legend=c("Dados", "Tendência Linear"), col=c("#c0392b", "blue"), lty=c(1,2))

# 3. Aplicação do Teste ADF Completo (Pacote urca)
# type = "trend": Testa considerando que pode haver uma tendência e uma constante.
# selectlags = "AIC": O R escolhe automaticamente quantos anos para trás olhar (lags).
teste_urca <- ur.df(dados_fev_check, type = "trend", selectlags = "AIC")

# 4. Exibição do Sumário Detalhado
summary(teste_urca)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.257 -13.270  -6.980   7.999  77.642 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 49.74471   11.35325   4.382 6.20e-05 ***
## z.lag.1     -1.26488    0.23017  -5.496 1.39e-06 ***
## tt           0.02244    0.21131   0.106    0.916    
## z.diff.lag  -0.02596    0.14355  -0.181    0.857    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.53 on 49 degrees of freedom
## Multiple R-squared:  0.6477, Adjusted R-squared:  0.6261 
## F-statistic: 30.03 on 3 and 49 DF,  p-value: 3.656e-11
## 
## 
## Value of test-statistic is: -5.4955 10.0686 15.1017 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2  6.50  4.88  4.16
## phi3  8.73  6.49  5.47

Resultados por Mês {.tabset .tabset-fade} A seguir, apresentamos a análise detalhada para cada um dos meses críticos da estação chuvosa.

#Janeiro

fit_jan <- analisar_clima(dados_jan, anos, "Janeiro", "#2980b9") # Azul

5.3.1 Análise Descritiva: Janeiro

Evento Máximo: 184.8 mm em 1986

#### Diagnóstico da Série Temporal

Tabela A: Diagnóstico Avançado - Janeiro

Teste	Estatistica	P_Valor	Diagnostico
Cox-Stuart (Tendência)	16.0000	0.4421	Sem Tendência
Mann-Kendall (Tendência)	0.1742	0.8617	Sem Tendência
Pettitt (Quebra)	186.0000	0.5873	Homogênea
ADF (Estacionariedade)	-3.5385	0.04644	Estacionária
Ljung-Box (Memória)	6.6114	0.7616	Independente

5.3.1.1 Modelagem GEV e Validação

Tabela B: Validação do Modelo GEV - Janeiro

Teste	Estatistica	P_Valor	Conclusao
Kolmogorov-Smirnov	0.0958	0.6938	Bom Ajuste
Anderson-Darling	0.4104	0.8377	Bom Ajuste

#### Estimativas de Risco (Nível de Retorno)

# O modelo parece superestimar os dados, uma investigação maior deve ser realizada.

#Fevereiro

fit_fev <- analisar_clima(dados_fev, anos, "Fevereiro", "#27ae60") # Verde

5.3.2 Análise Descritiva: Fevereiro

Evento Máximo: 119 mm em 2020

#### Diagnóstico da Série Temporal

Tabela A: Diagnóstico Avançado - Fevereiro

Teste	Estatistica	P_Valor	Diagnostico
Cox-Stuart (Tendência)	17.0000	0.2478	Sem Tendência
Mann-Kendall (Tendência)	-0.1452	0.8845	Sem Tendência
Pettitt (Quebra)	133.0000	1	Homogênea
ADF (Estacionariedade)	-2.8185	0.2447	Não Estacionária
Ljung-Box (Memória)	14.7417	0.1418	Independente

5.3.2.1 Modelagem GEV e Validação

Tabela B: Validação do Modelo GEV - Fevereiro

Teste	Estatistica	P_Valor	Conclusao
Kolmogorov-Smirnov	0.0718	0.9395	Bom Ajuste
Anderson-Darling	0.2799	0.9523	Bom Ajuste

#### Estimativas de Risco (Nível de Retorno)

# O ajuste para o mês de fevereiro no geral indica um bom ajuste de acordo com os gráficos e testes realizados.

6 Março

fit_mar <-analisar_clima(dados_mar, anos, "Março", "gold")

6.0.1 Análise Descritiva: Março

Evento Máximo: 110.4 mm em 2015

#### Diagnóstico da Série Temporal

Tabela A: Diagnóstico Avançado - Março

Teste	Estatistica	P_Valor	Diagnostico
Cox-Stuart (Tendência)	15.0000	0.7011	Sem Tendência
Mann-Kendall (Tendência)	1.2923	0.1962	Sem Tendência
Pettitt (Quebra)	248.0000	0.2264	Homogênea
ADF (Estacionariedade)	-3.6404	0.0379	Estacionária
Ljung-Box (Memória)	7.8612	0.6424	Independente

6.0.1.1 Modelagem GEV e Validação

Tabela B: Validação do Modelo GEV - Março

Teste	Estatistica	P_Valor	Conclusao
Kolmogorov-Smirnov	0.0758	0.91	Bom Ajuste
Anderson-Darling	0.2467	0.9721	Bom Ajuste

#### Estimativas de Risco (Nível de Retorno)

# O mês de março apresentou um bom ajuste, bem como indicou que o modelo se adequou bem aos dados de acordo com os testes realizados.

#Dezembro

fit_dez <- analisar_clima(dados_dez, anos, "Dezembro", "#e67e22") # Laranja

6.0.2 Análise Descritiva: Dezembro

Evento Máximo: 163.9 mm em 1979

#### Diagnóstico da Série Temporal

Tabela A: Diagnóstico Avançado - Dezembro

Teste	Estatistica	P_Valor	Diagnostico
Cox-Stuart (Tendência)	14.0000	0.845	Sem Tendência
Mann-Kendall (Tendência)	0.4647	0.6422	Sem Tendência
Pettitt (Quebra)	162.0000	0.7895	Homogênea
ADF (Estacionariedade)	-3.8931	0.02073	Estacionária
Ljung-Box (Memória)	13.5805	0.193	Independente

6.0.2.1 Modelagem GEV e Validação

Tabela B: Validação do Modelo GEV - Dezembro

Teste	Estatistica	P_Valor	Conclusao
Kolmogorov-Smirnov	0.1129	0.4845	Bom Ajuste
Anderson-Darling	0.9195	0.4019	Bom Ajuste

#### Estimativas de Risco (Nível de Retorno)

# Para o mês de dezembro o modelo aparenta subestimar os valores, uma investigação mais profunda também deve ser realizada. 

#Novembro

fit_nov <- analisar_clima(dados_nov, anos, "Novembro", "navy") #

6.0.3 Análise Descritiva: Novembro

Evento Máximo: 93 mm em 1975

#### Diagnóstico da Série Temporal

Tabela A: Diagnóstico Avançado - Novembro

Teste	Estatistica	P_Valor	Diagnostico
Cox-Stuart (Tendência)	10.0000	0.2478	Sem Tendência
Mann-Kendall (Tendência)	-2.2943	0.02178	Tendência (MK)
Pettitt (Quebra)	300.0000	0.08253	Homogênea
ADF (Estacionariedade)	-4.2874	0.01	Estacionária
Ljung-Box (Memória)	7.0596	0.7198	Independente

6.0.3.1 Modelagem GEV e Validação

Tabela B: Validação do Modelo GEV - Novembro

Teste	Estatistica	P_Valor	Conclusao
Kolmogorov-Smirnov	0.0748	0.9184	Bom Ajuste
Anderson-Darling	0.3567	0.8898	Bom Ajuste

#### Estimativas de Risco (Nível de Retorno)

Analisando a tendencia de Novembro

# ==============================================================================
# 1. PREPARAÇÃO E DIAGNÓSTICO
# ==============================================================================
# Criamos o dataframe primeiro
df_nov_trend <- data.frame(Ano = anos, Prec = dados_nov)
# Removemos NAs do DATAFRAME INTEIRO (preservando o alinhamento Ano-Prec)
df_nov_trend <- na.omit(df_nov_trend)

# Teste de Tendência e Estacionariedade
mk_res <- mk.test(df_nov_trend$Prec)
cs_res <- cox.stuart.test(df_nov_trend$Prec)
adf_res <- ur.df(df_nov_trend$Prec, type = "trend", selectlags = "AIC")

# Função para criar "Badges" coloridos
criar_badge <- function(p_valor) {
  cor <- if(p_valor < 0.05) "#c0392b" else "#27ae60" 
  texto <- if(p_valor < 0.05) "Significativo" else "Não Signif."
  paste0("<span style='color: white; background-color: ", cor, "; padding: 2px 6px; border-radius: 4px;'>", texto, "</span>")
}

df_diag <- data.frame(
  Teste = c("Mann-Kendall (Tendência)", "Cox-Stuart (Tendência)", "ADF (Estacionariedade)"),
  Estatistica = c(round(mk_res$statistic, 4), round(cs_res$statistic, 4), round(adf_res@teststat[1], 4)),
  P_Valor = c(format.pval(mk_res$p.value, digits=4), format.pval(cs_res$p.value, digits=4), "Ver Tab. Crítica"),
  Conclusao = c(criar_badge(mk_res$p.value), criar_badge(cs_res$p.value), "-")
)

print(kbl(df_diag, caption = "Tabela A: Diagnóstico Preliminar da Série (Novembro)", escape = FALSE) %>%
      kable_styling(bootstrap_options = c("hover", "condensed"), full_width = F, position = "center"))

Tabela A: Diagnóstico Preliminar da Série (Novembro)

	Teste	Estatistica	P_Valor	Conclusao
z	Mann-Kendall (Tendência)	-2.2943	0.02178	Significativo
statistic	Cox-Stuart (Tendência)	10.0000	0.2478	Não Signif.
	ADF (Estacionariedade)	-5.3930	Ver Tab. Crítica

cat("\n\n")

# ==============================================================================
# 1. MODELAGEM GEV (M0, M1 e M3) E SELEÇÃO - NOVEMBRO
# ==============================================================================
# Criamos o dataframe primeiro
df_nov_trend <- data.frame(Ano = anos, Prec = dados_nov)
# Removemos NAs mantendo alinhamento
df_nov_trend <- na.omit(df_nov_trend)

# --- Ajuste dos Modelos ---
# M0: Estacionário
fit0 <- fevd(Prec, data = df_nov_trend, type = "GEV", units = "mm")

# M1: Não-Estacionário em Locação (Mu ~ Ano)
fit1 <- fevd(Prec, data = df_nov_trend, location.fun = ~Ano, type = "GEV", units = "mm")

# M3: Não-Estacionário em Locação e Escala (Mu ~ Ano, Sigma ~ Ano)
fit3 <- fevd(Prec, data = df_nov_trend, location.fun = ~Ano, scale.fun = ~Ano, type = "GEV", units = "mm")

# --- Testes de Razão de Verossimilhança (LRT) Sequenciais ---
# 1. M0 vs M1 (Existe tendência na média?)
lrt_0_1 <- lr.test(fit0, fit1)
p_0_1   <- lrt_0_1$p.value

# 2. M1 vs M3 (Existe tendência na variabilidade, dado que ajustamos a média?)
lrt_1_3 <- lr.test(fit1, fit3)
p_1_3   <- lrt_1_3$p.value

# --- Lógica de Decisão Automática ---
# Define qual objeto será o 'fit_nov_prec' usado nas próximas etapas
modelo_selecionado <- "M0"
fit_nov_prec <- fit0 # Default
cor_decisao <- "#95a5a6" # Cinza (neutro)

if(p_0_1 < 0.05) {
  # Se M1 venceu M0, testamos se M3 vence M1
  if(p_1_3 < 0.05) {
    modelo_selecionado <- "M3"
    fit_nov_prec <- fit3
    msg_conclusao <- "O modelo **M3 (Locação e Escala variáveis)** apresentou o melhor ajuste (menor AIC/BIC e P-valor < 0.05)."
  } else {
    modelo_selecionado <- "M1"
    fit_nov_prec <- fit1
    msg_conclusao <- "O modelo **M1 (Apenas Locação variável)** foi o mais adequado. A variabilidade manteve-se constante."
  }
} else {
  msg_conclusao <- "Não há evidências estatísticas suficientes para rejeitar a estacionariedade. O modelo **M0 (Estacionário)** foi mantido."
}

# --- Tabela de Comparação (Agora com AIC) ---
df_comp_modelos <- data.frame(
  Modelo = c("M0: Estacionário", "M1: Locação Var.", "M3: Loc. e Escala Var."),
  LogLik = c(round(fit0$results$value, 2), round(fit1$results$value, 2), round(fit3$results$value, 2)),
  # Extraindo AIC e BIC dos sumários
  AIC = c(round(summary(fit0)$AIC, 2), round(summary(fit1)$AIC, 2), round(summary(fit3)$AIC, 2)),
  BIC = c(round(summary(fit0)$BIC, 2), round(summary(fit1)$BIC, 2), round(summary(fit3)$BIC, 2)),
  P_Valor_LRT = c("-", format.pval(p_0_1, digits=3), format.pval(p_1_3, digits=3)),
  Status = c(
    if(modelo_selecionado == "M0") "SELECIONADO" else "Rejeitado",
    if(modelo_selecionado == "M1") "SELECIONADO" else "Rejeitado",
    if(modelo_selecionado == "M3") "SELECIONADO" else "Rejeitado"
  )
)

fevd(x = Prec, data = df_nov_trend, type = “GEV”, units = “mm”)

[1] “Estimation Method used: MLE”

Negative Log-Likelihood Value: 241.4625

Estimated parameters: location scale shape 42.9418708 18.4396956 -0.1940002

Standard Error Estimates: location scale shape 2.8381210 2.0473277 0.1144137

Estimated parameter covariance matrix. location scale shape location 8.0549311 1.1627683 -0.14089164 scale 1.1627683 4.1915509 -0.12805345 shape -0.1408916 -0.1280534 0.01309049

AIC = 488.925

BIC = 494.947

fevd(x = Prec, data = df_nov_trend, location.fun = ~Ano, type = “GEV”, units = “mm”)

[1] “Estimation Method used: MLE”

Negative Log-Likelihood Value: 241.0682

Estimated parameters: mu0 mu1 scale shape 94.9357275 -0.0261335 18.4247327 -0.1967721

Standard Error Estimates: mu0 mu1 scale shape 436.2967845 0.2187167 2.8419444 0.1326587

Estimated parameter covariance matrix. mu0 mu1 scale shape mu0 190354.88418 -95.42337657 -836.1573820 27.17840973 mu1 -95.42338 0.04783698 0.4197267 -0.01369492 scale -836.15738 0.41972672 8.0766481 -0.25925620 shape 27.17841 -0.01369492 -0.2592562 0.01759832

AIC = 490.1364

BIC = 498.1657

fevd(x = Prec, data = df_nov_trend, location.fun = ~Ano, scale.fun = ~Ano, type = “GEV”, units = “mm”)

[1] “Estimation Method used: MLE”

Negative Log-Likelihood Value: 240.6845

Estimated parameters: mu0 mu1 sigma0 sigma1 shape 125.77095611 -0.04142993 63.82452280 -0.02268187 -0.20003887

AIC = 491.3689

BIC = 501.4056

fevd(x = Prec, data = df_nov_trend, type = “GEV”, units = “mm”)

[1] “Estimation Method used: MLE”

Negative Log-Likelihood Value: 241.4625

Estimated parameters: location scale shape 42.9418708 18.4396956 -0.1940002

Standard Error Estimates: location scale shape 2.8381210 2.0473277 0.1144137

Estimated parameter covariance matrix. location scale shape location 8.0549311 1.1627683 -0.14089164 scale 1.1627683 4.1915509 -0.12805345 shape -0.1408916 -0.1280534 0.01309049

AIC = 488.925

BIC = 494.947

fevd(x = Prec, data = df_nov_trend, location.fun = ~Ano, type = “GEV”, units = “mm”)

[1] “Estimation Method used: MLE”

Negative Log-Likelihood Value: 241.0682

Estimated parameters: mu0 mu1 scale shape 94.9357275 -0.0261335 18.4247327 -0.1967721

Standard Error Estimates: mu0 mu1 scale shape 436.2967845 0.2187167 2.8419444 0.1326587

Estimated parameter covariance matrix. mu0 mu1 scale shape mu0 190354.88418 -95.42337657 -836.1573820 27.17840973 mu1 -95.42338 0.04783698 0.4197267 -0.01369492 scale -836.15738 0.41972672 8.0766481 -0.25925620 shape 27.17841 -0.01369492 -0.2592562 0.01759832

AIC = 490.1364

BIC = 498.1657

fevd(x = Prec, data = df_nov_trend, location.fun = ~Ano, scale.fun = ~Ano, type = “GEV”, units = “mm”)

[1] “Estimation Method used: MLE”

Negative Log-Likelihood Value: 240.6845

Estimated parameters: mu0 mu1 sigma0 sigma1 shape 125.77095611 -0.04142993 63.82452280 -0.02268187 -0.20003887

AIC = 491.3689

BIC = 501.4056

print(kbl(df_comp_modelos, caption = "Tabela B: Seleção de Modelo GEV (Critérios de Informação)", align = "c") %>%
      kable_styling(bootstrap_options = "striped", full_width = F, position = "center") %>%
      row_spec(which(df_comp_modelos$Status == "SELECIONADO"), bold = T, color = "white", background = "#27ae60"))

Tabela B: Seleção de Modelo GEV (Critérios de Informação)

Modelo	LogLik	AIC	BIC	P_Valor_LRT	Status
M0: Estacionário	241.46	488.93	494.95		SELECIONADO
M1: Locação Var.	241.07	490.14	498.17	0.375	Rejeitado
M3: Loc. e Escala Var.	240.68	491.37	501.41	0.381	Rejeitado

cat(paste0("\n<div class='alert-info'>**CONCLUSÃO:** ", msg_conclusao, "</div>\n"))

CONCLUSÃO: Não há evidências estatísticas suficientes para rejeitar a estacionariedade. O modelo M0 (Estacionário) foi mantido.

# ==============================================================================
# DIAGNÓSTICO E PLOTAGEM DO VENCEDOR
# ==============================================================================
par(mfrow=c(1,2))
plot(fit_nov_prec, type = "density", main = paste("Densidade:", modelo_selecionado))
plot(fit_nov_prec, type = "qq", main = paste("QQ-Plot:", modelo_selecionado))

par(mfrow=c(1,1))

# Gráfico de Retorno (Evolução Temporal)
# Plota apenas se o modelo selecionado NÃO for o estacionário
if(modelo_selecionado != "M0") {
  cat("\n#### Evolução do Risco Climático (1968-2022)\n")
  
  # Calcula os níveis de retorno para o período dos dados
  # Nota: Com scale.fun, o return.level retorna uma matriz variando ano a ano
  rl_vals <- return.level(fit_nov_prec, return.period = c(2, 20, 100))
  
  matplot(df_nov_trend$Ano, rl_vals, type = "l", lty = 1, lwd = 2,
          col = c("#27ae60", "#e67e22", "#c0392b"),
          main = paste0("Níveis de Retorno Efetivos (", modelo_selecionado, ")"),
          ylab = "Precipitação Máxima (mm)", xlab = "Ano")
  
  legend("topleft", legend = c("Retorno de 2 Anos", "Retorno de 20 Anos", "Retorno de 100 Anos"), 
         col = c("#27ae60", "#e67e22", "#c0392b"), lty = 1, lwd = 2, bty = "n", cex = 0.9)
  grid(col = "lightgray", lty = "dotted")
}
cat("\n\n")

7 Identificando as distribuições

# ==============================================================================
# TABELA DE CARACTERIZAÇÃO DAS DISTRIBUIÇÕES
# ==============================================================================
# Lista dos modelos (Certifique-se que fit_nov_prec é o vencedor do bloco anterior)
lista_modelos <- list(
  "Novembro"  = fit_nov_prec,
  "Dezembro"  = fit_dez,
  "Janeiro"   = fit_jan,
  "Fevereiro" = fit_fev,
  "Março"     = fit_mar
)

df_tipos <- data.frame(
  Mes = character(), Xi = character(), IC_95 = character(), 
  Tipo = character(), stringsAsFactors = FALSE
)

for (mes in names(lista_modelos)) {
  modelo <- lista_modelos[[mes]]
  
  # Extrai IC dos parâmetros (Robustez para encontrar o 'shape' ou 'xi')
  ci_params <- ci(modelo, type = "parameter", method = "normal", verbose = FALSE)
  
  # Busca inteligente pelo índice do parâmetro de forma
  idx_shape <- grep("shape|xi", rownames(ci_params), ignore.case = TRUE)
  
  if(length(idx_shape) > 0) {
    est <- ci_params[idx_shape, 2]
    low <- ci_params[idx_shape, 1]
    upp <- ci_params[idx_shape, 3]
    
    # Classificação
    if (low > 0) {
      tipo_dist <- "Fréchet (Cauda Pesada)"
      cor_tipo <- "#c0392b" # Vermelho
    } else if (upp < 0) {
      tipo_dist <- "Weibull (Cauda Limitada)"
      cor_tipo <- "#27ae60" # Verde
    } else {
      tipo_dist <- "Gumbel (Cauda Leve)"
      cor_tipo <- "#2980b9" # Azul
    }
    
    df_tipos <- rbind(df_tipos, data.frame(
      Mes = mes,
      Xi = round(est, 3),
      IC_95 = paste0("[", round(low, 3), "; ", round(upp, 3), "]"),
      Tipo = paste0("<span style='color:", cor_tipo, "; font-weight:bold'>", tipo_dist, "</span>")
    ))
  }
}

df_tipos %>%
  kbl(caption = "<div style='text-align:center; color:#2c3e50; font-size:18px; font-weight:bold'>Caracterização das Distribuições Ajustadas</div>", 
      align = "c", escape = FALSE) %>%
  kable_styling(bootstrap_options = c("hover", "striped"), full_width = F, position = "center") %>%
  column_spec(1, bold = T) %>%
  footnote(general = "Classificação baseada no IC (95%) do parâmetro de forma. Se inclui zero, assume-se Gumbel.", general_title = "Nota: ")

Caracterização das Distribuições Ajustadas

Mes	Xi	IC_95	Tipo
Novembro	-0.194	[-0.418; 0.03]	Gumbel (Cauda Leve)
Dezembro	0.000	[-0.132; 0.131]	Gumbel (Cauda Leve)
Janeiro	0.232	[-0.038; 0.502]	Gumbel (Cauda Leve)
Fevereiro	0.167	[-0.071; 0.404]	Gumbel (Cauda Leve)
Março	-0.032	[-0.208; 0.144]	Gumbel (Cauda Leve)
Nota:
Classificação baseada no IC (95%) do parâmetro de forma. Se inclui zero, assume-se Gumbel.

#Comparativo e Conclusão

Comparação das curvas de Nível de Retorno para os três meses analisados. Este gráfico permite identificar qual mês apresenta maior risco de eventos extremos para um mesmo período de retorno.

# Carrega bibliotecas de estética
if(!require(kableExtra)) install.packages("kableExtra")
if(!require(extRemes)) install.packages("extRemes")
library(kableExtra)
library(extRemes)

# 1. Definição dos Períodos
periodos <- c(5, 10, 20, 50, 100)

# ==============================================================================
# FUNÇÕES BLINDADAS (Para suportar modelos com e sem tendência)
# ==============================================================================

# Função para formatar Texto com IC (HTML)
get_est_ci_smart <- function(modelo, periods, ano_ref = 2022) {
  
  resultado <- tryCatch({
    if (!is.null(modelo$cov.data)) {
      # Se tem tendência (ex: Novembro), cria a matriz qcov para o ano 2022
      dados_modelo <- modelo$cov.data
      nome_col <- colnames(dados_modelo)[1]
      novos_dados <- data.frame(ano_ref)
      colnames(novos_dados) <- nome_col
      
      matriz_covariavel <- make.qcov(modelo, vals = novos_dados)
      ci(modelo, return.period = periods, method = "normal", verbose = FALSE, qcov = matriz_covariavel)
    } else {
      # Se é estacionário (ex: Jan/Fev/Mar)
      ci(modelo, return.period = periods, method = "normal", verbose = FALSE)
    }
  }, error = function(e) {
    # Fallback de segurança
    return(ci(modelo, return.period = periods, method = "normal", verbose = FALSE))
  })
  
  est <- round(resultado[,2], 1)
  low <- round(resultado[,1], 1)
  upp <- round(resultado[,3], 1)
  
  paste0("<strong>", est, "</strong><br><span style='font-size:11px; color:#7f8c8d'>[", low, "; ", upp, "]</span>")
}

# Função para extrair valores numéricos (Para o gráfico)
get_values_smart <- function(modelo, periods, ano_ref = 2022) {
  res <- tryCatch({
    if (!is.null(modelo$cov.data)) {
      dados_modelo <- modelo$cov.data
      nome_col <- colnames(dados_modelo)[1]
      novos_dados <- data.frame(ano_ref)
      colnames(novos_dados) <- nome_col
      
      matriz_covariavel <- make.qcov(modelo, vals = novos_dados)
      return.level(modelo, return.period = periods, qcov = matriz_covariavel)
    } else {
      return.level(modelo, return.period = periods)
    }
  }, error = function(e) return.level(modelo, return.period = periods))
  return(as.numeric(res))
}

# 3. Gerando as colunas formatadas (Texto)
col_nov <- get_est_ci_smart(fit_nov_prec, periodos, ano_ref = 2022)
col_dez <- get_est_ci_smart(fit_dez, periodos, ano_ref = 2022) 
col_jan <- get_est_ci_smart(fit_jan, periodos)
col_fev <- get_est_ci_smart(fit_fev, periodos)
col_mar <- get_est_ci_smart(fit_mar, periodos) # Inclusão de Março

# 4. Construção do DataFrame Final (Ordem Cronológica)
df_final <- data.frame(
  "Retorno"   = paste(periodos, "Anos"),
  "Novembro"  = col_nov,
  "Dezembro"  = col_dez,
  "Janeiro"   = col_jan,
  "Fevereiro" = col_fev,
  "Março"     = col_mar
)

# 5. Lógica do Vencedor (Numérica)
vals_nov <- get_values_smart(fit_nov_prec, periodos, ano_ref = 2022)
vals_dez <- get_values_smart(fit_dez, periodos, ano_ref = 2022)
vals_jan <- get_values_smart(fit_jan, periodos)
vals_fev <- get_values_smart(fit_fev, periodos)
vals_mar <- get_values_smart(fit_mar, periodos)

vencedores <- c()
for(i in 1:length(periodos)){
  linha <- c(Novembro=vals_nov[i], Dezembro=vals_dez[i], Janeiro=vals_jan[i], Fevereiro=vals_fev[i], Março=vals_mar[i])
  vencedores[i] <- names(linha)[which.max(linha)]
}
df_final$Risco_Maior <- vencedores

# 6. Tabela Estilizada
df_final %>%
  kbl(caption = "<div style='text-align:center; color:#2c3e50; font-size:18px; font-weight:bold'>Precipitação Máxima (mm): Comparativo Completo (Nov-Mar)</div>", 
      align = "c", escape = FALSE, row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("hover", "condensed"), full_width = F, position = "center", font_size = 14) %>%
  add_header_above(c(" " = 1, "Estação Chuvosa" = 5, "Conclusão" = 1)) %>% # Agora são 5 meses
  column_spec(1, bold = T, background = "#f8f9fa", width = "90px") %>%
  column_spec(7, bold = T, color = "white", background = "#e74c3c", width = "110px") %>% # Coluna Vencedor
  row_spec(0, bold = T, color = "white", background = "#34495e") %>%
  footnote(general = "Para os meses com tendência (Nov/Dez), os valores referem-se ao risco estimado para o ano de 2022.", general_title = "Nota: ")

Precipitação Máxima (mm): Comparativo Completo (Nov-Mar)

	Estação Chuvosa					Conclusão
Retorno	Novembro	Dezembro	Janeiro	Fevereiro	Março	Risco_Maior
5 Anos	66.9 [60.3; 73.6]	72.4 [64; 80.8]	79 [64.4; 93.6]	54.5 [45.1; 63.8]	61.5 [53.2; 69.8]	Janeiro
10 Anos	76.6 [68.7; 84.4]	85.8 [74.5; 97.1]	103.7 [78.7; 128.7]	70 [55; 85.1]	74.5 [63.3; 85.7]	Janeiro
20 Anos	84.6 [74.3; 94.9]	98.7 [83.5; 113.9]	131.8 [88.9; 174.8]	86.9 [62.6; 111.2]	86.7 [71.3; 102.1]	Janeiro
50 Anos	93.4 [78.4; 108.4]	115.4 [93.3; 137.5]	176 [94.3; 257.6]	111.9 [68.5; 155.4]	102.1 [79; 125.2]	Janeiro
100 Anos	99.1 [79.9; 118.2]	127.9 [99.3; 156.5]	215.9 [90.1; 341.6]	133.4 [69.3; 197.6]	113.3 [82.9; 143.8]	Janeiro
Nota:
Para os meses com tendência (Nov/Dez), os valores referem-se ao risco estimado para o ano de 2022.

# 7. Gráfico Comparativo (5 Meses)
# Define limites do eixo Y considerando TODOS os meses
todos_valores <- c(vals_nov, vals_dez, vals_jan, vals_fev, vals_mar)
y_min <- min(todos_valores) * 0.9
y_max <- max(todos_valores) * 1.1

plot(periodos, vals_jan, type="n", 
     ylim=c(y_min, y_max), log="x",
     main="Curvas de Risco: Comparativo Sazonal",
     xlab="Período de Retorno (Anos)", ylab="Precipitação (mm)", axes=FALSE)

grid(nx=NULL, ny=NULL, col="lightgray", lty="dotted")
axis(1, at=periodos, labels=paste(periodos, "anos")); axis(2, las=2)

# --- Plotagem das Linhas ---
# Novembro (Roxo)
lines(periodos, vals_nov, col="#8e44ad", lwd=3, lty=1)
points(periodos, vals_nov, pch=18, col="#8e44ad", cex=1.5)

# Dezembro (Laranja Escuro)
lines(periodos, vals_dez, col="#d35400", lwd=3, lty=5) 
points(periodos, vals_dez, pch=15, col="#d35400", cex=1.3)

# Janeiro (Azul)
lines(periodos, vals_jan, col="#2980b9", lwd=3, lty=1)
points(periodos, vals_jan, pch=19, col="#2980b9", cex=1.3)

# Fevereiro (Verde)
lines(periodos, vals_fev, col="#27ae60", lwd=3, lty=2)
points(periodos, vals_fev, pch=17, col="#27ae60", cex=1.3)

# Março (Dourado/Amarelo - Novo!)
lines(periodos, vals_mar, col="#f39c12", lwd=3, lty=4)
points(periodos, vals_mar, pch=8, col="#f39c12", cex=1.5) # pch=8 é um asterisco/estrela

# Legenda Atualizada
legend("topleft", 
       legend=c("Novembro", "Dezembro", "Janeiro", "Fevereiro", "Março"),
       col=c("#8e44ad", "#d35400", "#2980b9", "#27ae60", "#f39c12"), 
       lwd=3, pch=c(18, 15, 19, 17, 8), lty=c(1, 5, 1, 2, 4), bty="n", cex=0.85, ncol=2) # ncol=2 para economizar espaço vertical

8 Referências Bibliográficas

8.1 Softwares e Pacotes Utilizados

A análise estatística foi conduzida no ambiente R. Abaixo estão as citações oficiais dos softwares e bibliotecas aplicadas neste estudo:

R Core Team (2025). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.

Gilleland E, Katz RW (2016). “extRemes 2.0: An Extreme Value Analysis Package in R.” Journal of Statistical Software, 72(8), 1-39. doi:10.18637/jss.v072.i08 https://doi.org/10.18637/jss.v072.i08.

Pfaff B (2008). Analysis of Integrated and Cointegrated Time Series with R, Second edition. Springer, New York. ISBN 0-387-27960-1, https://www.pfaffikus.de.

Pohlert T (2023). trend: Non-Parametric Trend Tests and Change-Point Detection. doi:10.32614/CRAN.package.trend https://doi.org/10.32614/CRAN.package.trend, R package version 1.1.6, https://CRAN.R-project.org/package=trend.

Trapletti A, Hornik K (2024). tseries: Time Series Analysis and Computational Finance. doi:10.32614/CRAN.package.tseries https://doi.org/10.32614/CRAN.package.tseries, R package version 0.10-58, https://CRAN.R-project.org/package=tseries.

Referências

BADDOUR, Omar; KONTONGOMDE, H. The role of climatological normals in a changing climate. Geneva: World Meteorological Organization (WMO), 2007. Disponível em: <https://library.wmo.int/idurl/4/4166>.

CAI, Yuzhi; HAMES, Dominic. Minimum sample size determination for generalized extreme value distribution. Communications in Statistics-Simulation and Computation, v. 40, n. 1, p. 87–98, 2011.

COLES, Stuart. An Introduction to Statistical Modeling of Extreme Values. London: Springer, 2001.

COX, David R.; STUART, Alan. Some quick sign tests for trend in location and dispersion. Biometrika, v. 42, n. 1/2, p. 80–95, 1955.

DICKEY, David A.; FULLER, Wayne A. Distribution of the estimators for autoregressive time series with a unit root. Journal of the American Statistical Association, v. 74, n. 366a, p. 427–431, 1979.

INSTITUTO NACIONAL DE METEOROLOGIA. Banco de Dados Meteorológicos. Disponível em: <https://portal.inmet.gov.br/>, 2024.

JENKINSON, Arthur F. The frequency distribution of the annual maximum (or minimum) values of meteorological elements. Quarterly Journal of the Royal Meteorological Society, v. 81, n. 348, p. 158–171, 1955.

MORETTIN, Pedro A.; TOLOI, Clélia MC. Análise de Séries Temporais. 2. ed. São Paulo: Edgard Blücher, 2006.

ROSLAN, Razira Aniza; SU NA, Chin; GABDA, Darmesah. Parameter Estimations of the Generalized Extreme Value Distributions for Small Sample Size. Mathematics and Statistics, v. 8, n. 2A, p. 47–51, 2020.

UNIVERSIDADE FEDERAL DE VIÇOSA. Boletim Meteorológico 2023. Viçosa-MG: Departamento de Engenharia Agrícola. Estação Climatológica Principal de Viçosa, 2023.

WALD, Abraham; WOLFOWITZ, Jacob. On a test whether two samples are from the same population. The Annals of Mathematical Statistics, v. 11, n. 2, p. 147–162, 1940.