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

Figura 1: Vista da praça municipal de Viçosa-MG

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.

Os dados pluviométricos foram obtidos a partir das estações do Instituto Nacional de Meteorologia (INMET) 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).

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: Aplicado para verificar a existência de tendências determinísticas (monotonicidade). 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): 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.

  1. 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.
  2. 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
# 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 anual, 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.

# Importa especificamente a folha 2
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 e/ou estacionariedade nas séries bem como o teste de sequência para verificar a independência. Como a análise principal foca nos máximos semre do mesmo mês, não houve necessidade de analisar sazonalidade.

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

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

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

# Retorna a configuração gráfica ao padrão (1x1)
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 estrutura matemática do modelo GEV, apresenta-se a seguir com a implementação manual da função de Log-Verossimilhança e o cálculo dos parâmetros via otimização numérica (nlm), além do cálculo de incerteza via Método Delta.

A função de distribuição acumulada GEV é 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 utilozaçã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

O ajuste abaixo utiliza o método de máxima verossimilhança otimizado.

# 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 em Perfil de Verossimilhança (Profile Likelihood), que são mais precisos para extremos.

par(mfrow=c(1,3))

# Nota: O argumento 'xrange' foi ajustado visualmente para cobrir a região de verossimilhança
# 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 Temperatura máxima

6 Análise de Temperaturas Máximas

Nesta seção, replicamos a metodologia estatística para analisar o comportamento das temperaturas máximas extremas registradas no mês de janeiro ao longo da série histórica.

#Análise Exploratória

Visualização inicial da série temporal e distribuição dos dados para identificação de padrões visuais.

# Definição das variáveis
TEMP_MAX <- EXTREMOS$TEMP
year <- seq(1968, 2022, 1)

# Plotagem
plot(TEMP_MAX ~ year, type="b", main="Série Histórica: Temperatura (Jan)", 
     ylab="Temperatura (ºC)", xlab="Ano", col="#c0392b", pch=19, lwd=2)
grid()

hist(TEMP_MAX, main="Distribuição de Frequências", 
     xlab="Temperatura (ºC)", col="#e74c3c", border="white", prob=TRUE)
lines(density(TEMP_MAX), col="#2c3e50", lwd=2)

# 1. Teste Augmented Dickey-Fuller (Estacionariedade)
# H0: Série Não-Estacionária | H1: Série Estacionária
teste_completo2 <- ur.df(TEMP_MAX, type = "trend", selectlags = "AIC")
summary(teste_completo2)
## 
## ############################################### 
## # 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 
## -2.8117 -0.9867 -0.1703  1.1349  3.4481 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 27.31149    5.98963   4.560 3.44e-05 ***
## z.lag.1     -0.85765    0.18659  -4.597 3.04e-05 ***
## tt           0.02269    0.01456   1.558    0.126    
## z.diff.lag  -0.05061    0.13743  -0.368    0.714    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.548 on 49 degrees of freedom
## Multiple R-squared:  0.4754, Adjusted R-squared:  0.4433 
## F-statistic:  14.8 on 3 and 49 DF,  p-value: 5.458e-07
## 
## 
## Value of test-statistic is: -4.5965 7.0883 10.5841 
## 
## 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
# 2. Teste de Cox-Stuart (Tendência)
# H0: Sem tendência | H1: Com tendência
print(cox.stuart.test(TEMP_MAX))
## 
##  Cox Stuart test
## 
## data:  TEMP_MAX
## statistic = 18, n = 27, p-value = 0.1221
## alternative hypothesis: non randomness
fit_ismev2 <- gev.fit(TEMP_MAX)
## $conv
## [1] 0
## 
## $nllh
## [1] 101.8466
## 
## $mle
## [1] 32.0030452  1.4402150 -0.1722711
## 
## $se
## [1] 0.2209966 0.1589081 0.1110962
gev.diag(fit_ismev2)

6.1 Ajuste e Validação via extRemes

O ajuste abaixo utiliza o método de máxima verossimilhança otimizado.

# Ajuste do Modelo
fit2 <- fevd(TEMP_MAX, units = "deg C")
print(fit2)
## 
## fevd(x = TEMP_MAX, units = "deg C")
## 
## [1] "Estimation Method used: MLE"
## 
## 
##  Negative Log-Likelihood Value:  101.8466 
## 
## 
##  Estimated parameters:
##   location      scale      shape 
## 32.0031274  1.4403017 -0.1722879 
## 
##  Standard Error Estimates:
##  location     scale     shape 
## 0.2210030 0.1589272 0.1111127 
## 
##  Estimated parameter covariance matrix.
##              location        scale        shape
## location  0.048842330  0.007655014 -0.010397961
## scale     0.007655014  0.025257849 -0.009215418
## shape    -0.010397961 -0.009215418  0.012346041
## 
##  AIC = 209.6931 
## 
##  BIC = 215.7151

6.1.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 do objeto fit2
params <- fit2$results$par
loc_est <- params["location"]
scale_est <- params["scale"]
shape_est <- params["shape"]

# Teste KS
ks_result2 <- ks.test(TEMP_MAX, "pevd", loc=loc_est, scale=scale_est, shape=shape_est, type="GEV")
print(ks_result2)
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  TEMP_MAX
## D = 0.098161, p-value = 0.6643
## alternative hypothesis: two-sided
# Teste AD (Anderson-Darling)
ad_result2 <- ad.test(TEMP_MAX, "pevd", loc=loc_est, scale=scale_est, shape=shape_est, type="GEV")
print(ad_result2)
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: distribution 'pevd'
##  with parameters loc = 32.0031274408662, scale = 1.44030173589631, shape
##  = -0.172287886670032, type = GEV
##  Parameters assumed to be fixed
## 
## data:  TEMP_MAX
## An = 0.50548, p-value = 0.7405

6.2 Estimativa de Níveis de Retorno

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

par(mfrow=c(1,3))

# Nível de Retorno: 4 Anos
ci(fit2, method = "proflik", xrange = c(33, 36), verbose = FALSE, return.period = 4, main = "Retorno 4 Anos")
## fevd(x = TEMP_MAX, units = "deg C")
## 
## [1] "Profile Likelihood"
## 
## [1] "4-year return level: 33.618"
## 
## [1] "95% Confidence Interval: (33.127, 34.1679)"
# Nível de Retorno: 10 Anos
ci(fit2, method = "proflik", xrange = c(34, 38), verbose = FALSE, return.period = 10, main = "Retorno 10 Anos")
## fevd(x = TEMP_MAX, units = "deg C")
## 
## [1] "Profile Likelihood"
## 
## [1] "10-year return level: 34.69"
## 
## [1] "95% Confidence Interval: (34.1333, 35.5805)"
# Nível de Retorno: 20 Anos
ci(fit2, method = "proflik", xrange = c(35, 41), verbose = FALSE, return.period = 20, main = "Retorno 20 Anos")
## fevd(x = TEMP_MAX, units = "deg C")
## 
## [1] "Profile Likelihood"
## 
## [1] "20-year return level: 35.352"
## 
## [1] "95% Confidence Interval: (35, 36.6996)"

7 Referências Bibliográficas

7.1 Literatura e Fontes de Dados

  • UNIVERSIDADE FEDERAL DE VIÇOSA – UFV. Departamento de Engenharia Agrícola. Estação Climatológica Principal de Viçosa. Boletim meteorológico 2023. Viçosa-MG: UFV, 2023.
  • MOOD, A. M.; GRAYBILL, F. A.; BOES, D. C. Introduction to the Theory of Statistics. 3rd. ed. [S.l.]: McGraw-Hill,1974.
  • INSTITUTO NACIONAL DE METEOROLOGIA (INMET). Banco de Dados Meteorológicos. Brasília. Disponível em: https://portal.inmet.gov.br/. Acesso em: 27 jan. 2026.

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