trabalhoFinalMat3

Trabalho final da disciplina de Economia Matemática 3 Estimação dos preços históricos do milho para o “Cobweb Model”

Instalando pacotes para preço e quantidade de milho

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(ggplot2)
library(readxl)
 
# Fator de conversão: 1000 toneladas * 1000 (kg/ton) / 60 (kg/saca)
FATOR_CONVERSAO <- 1000000 / 60 

# Carregar Preços (CEPEA)
preco <- read_xlsx("CEPEA.xlsx", skip = 3) 

# Carregar Produção (MilhoTotal)
producao <- read_xlsx("MilhoTotalSerieHist.xlsx") 

LIMPEZA E CÁLCULO DA MÉDIA ANUAL DO PREÇO (CEPEA)

# 1. Leitura do arquivo (Mantida)
preco <- read_excel("CEPEA.xlsx", skip = 3) 

# 2. TRATAMENTO E CONVERSÃO DE DATA (Forçado fora do pipe)
# O read_excel pode ter lido a data como número serial (comum no Excel) ou como texto.

# Tentativa 1: Converter a coluna 'Data' assumindo que ela é um número serial do Excel.
# (Se for texto, resultará em NA, que será corrigido na próxima linha.)
preco$Data_Convertida <- as.Date(as.numeric(preco$Data), origin = "1899-12-30")

## Warning in as.Date(as.numeric(preco$Data), origin = "1899-12-30"): NAs
## introduzidos por coerção

# Tentativa 2: Para os valores que FALHARAM (NA) na Tentativa 1, tenta converter como texto Dia/Mês/Ano (dmy).
preco$Data_Convertida[is.na(preco$Data_Convertida)] <- dmy(preco$Data[is.na(preco$Data_Convertida)])

# 3. Cria a coluna Ano a partir da coluna de Data convertida.
preco$Ano <- year(preco$Data_Convertida) 

# 4. Filtrar a partir de 2003 e calcular a Média Anual (USANDO O NOVO DATAFRAME)
preco_anual <- preco %>%
  # Garante que filtra apenas onde a conversão para Ano foi bem-sucedida e o ano é >= 2003
  filter(!is.na(Ano), Ano >= 2003) %>% 
  group_by(Ano) %>%
  summarise(
    Preco_Medio_Anual = mean(`À vista R$`, na.rm = TRUE), 
    .groups = 'drop'
  )

PREPARAÇÃO, PIVOTAGEM E CONVERSÃO DA PRODUÇÃO (MILHOTOTAL)

# Fator de conversão: 1000 toneladas * 1000 (kg/ton) / 60 (kg/saca)
FATOR_CONVERSAO <- 1000000 / 60 

# Leitura do arquivo (Mantida: Renomeada de df_producao para producao)
producao <- read_excel("MilhoTotalSerieHist.xlsx") 
library(tidyr)

# Pivotar os dados (transformar de colunas de anos para uma coluna de Ano)
producao_long <- producao %>%
  # A coluna 'Anos' é o identificador, todas as outras são medidas (produção)
  pivot_longer(
    cols = -Anos, # Pivotar todas as colunas exceto a primeira ('Anos')
    names_to = "Ano_str",
    values_to = "Producao_1000Ton"
  ) %>%
  # Limpar e converter o ano para numérico
  mutate(
    Ano = as.integer(Ano_str), 
    Producao_1000Ton = as.numeric(Producao_1000Ton)
  ) %>%
  # Filtra anos a partir de 2003
  filter(!is.na(Ano), Ano >= 2003) 

# 3.2 Conversão de 1000 ton para sacas de 60kg
producao_limpa <- producao_long %>%
  mutate(
    Producao_Sacas_60kg = Producao_1000Ton * FATOR_CONVERSAO
  ) %>%
  select(Ano, Producao_Sacas_60kg)

COMBINAR AS SÉRIES E VISUALIZAR SÉRIE TEMPORAL

dados_milho_anual <- inner_join(preco_anual, producao_limpa, by = "Ano")
dados_milho_anual <- na.omit(dados_milho_anual)

# Visualização da Série Temporal do Preço
ggplot(data = dados_milho_anual, aes(x = Ano, y = Preco_Medio_Anual)) +
  geom_line(color = "blue") + geom_point(color = "darkblue") +
  labs(x = "Ano", y = "Preço Médio do Milho (R$/sc)") +
  theme_minimal()

ADAPTAÇÃO DO MODELO DA TEIA DE ARANHA (COBWEB MODEL)

#Calculando Pontos de Equilíbrio (P* e Q*)
preco_equilibrio_sugerido <- mean(dados_milho_anual$Preco_Medio_Anual, na.rm = TRUE)
Q_media_original <- mean(dados_milho_anual$Producao_Sacas_60kg, na.rm = TRUE)

#NORMALIZAÇÃO: Convertendo para Milhões de Sacas
Q_media_normalizada <- Q_media_original / FATOR_CONVERSAO

cat(paste("P* (Preço de Equilíbrio): R$", round(preco_equilibrio_sugerido, 2), "/sc\n"))

## P* (Preço de Equilíbrio): R$ 40.09 /sc

cat(paste("Q* (Quantidade de Equilíbrio):", round(Q_media_normalizada, 2), " Milhões de Sacas\n"))

## Q* (Quantidade de Equilíbrio): 83929.35  Milhões de Sacas

#Parâmetros de calibração na escala normalizada:
# b_valor: Declive da Demanda (quanto maior, mais horizontal/inelástica)
b_valor <- 5.0
# beta_valor: Controla a inclinação da Oferta (1/beta). Quanto MENOR beta, MAIS vertical/elástica a Oferta.
beta_valor <- 2.5 

lambda_valor <- 0.05 # Mantido no list, mas não mais usado na supply

#Recalculando 'a' para garantir Q_d = Q* em P=P*
a_recalibrado <- Q_media_normalizada + b_valor * preco_equilibrio_sugerido

#Definindo o modelo com Oferta Linear
Market_Milho <- function(a = a_recalibrado, 
                         b = b_valor, 
                         c = preco_equilibrio_sugerido, 
                         d = Q_media_normalizada, 
                         beta = beta_valor,
                         lambda = lambda_valor) { 
  
  list(
    a = a, b = b, c = c, d = d, beta = beta, lambda = lambda,
    demand = function(p) a - b * p, 
    # OFERTA LINEAR E ELÁSTICA
    supply = function(p) d + (1 / beta) * (p - c) 
  )
}

market_milho <- Market_Milho()

#Plot Oferta e Demanda (Gráfico com inclinações visíveis)
p_grid_milho <- seq(20, 50, length.out = 10) 
demand_milho <- market_milho$demand(p_grid_milho)
supply_milho <- market_milho$supply(p_grid_milho)
df_market_milho <- data.frame(p_grid = p_grid_milho, demand = demand_milho, supply = supply_milho)

ggplot(df_market_milho, aes(x = p_grid)) +
  geom_line(aes(y = demand, color = "Demanda")) +
  geom_line(aes(y = supply, color = "Oferta")) +
  labs(x = "Quantidade (Milhões de Sacas de 60kg)", y = "Preço (R$/sc)") +  
  theme_minimal()

FUNÇÕES E GRÁFICOS DINÂMICOS

# Função g para calcular o próximo preço (Cobweb Dynamics)
# P_t = (a - Q_s(P_{t-1})) / b
g <- function(model, current_price) {
  # Usando a função de Oferta LINEAR: Q_s(P_{t-1}) = d + (1 / beta) * (P_{t-1} - c)
  current_supply <- model$d + (1 / model$beta) * (current_price - model$c)
  # Preço de equilíbrio na Demanda: P_t = (a - Q_s) / b
  (model$a - current_supply) / model$b
}

Função para plotar a linha de 45 graus (Diagrama de Fase)

  plot45 <- function(model, pmin, pmax, p0, num_arrows = 5) {
    p_grid <- seq(pmin, pmax, length.out = 150)
  g_values <- sapply(p_grid, function(p) g(model, p))
  
  df_45 <- data.frame(p_grid, g_values)
  
  ggplot(df_45, aes(x = p_grid)) +
    geom_line(aes(y = g_values, color = "g"), size = 1) +
    geom_abline(slope = 1, intercept = 0, color = "darkgreen", linetype = "dashed") +
    # Ajustando os limites para focar no intervalo de preço relevante
    xlim(c(pmin, pmax)) + ylim(c(pmin, pmax)) +
    labs(x = expression(p[t]), y = expression(p[t+1])) +
    theme_minimal() +
    theme(legend.position = "none")
}

EXIBINDO O DIAGRAMA DE FASE

plot45(market_milho, pmin = 20, pmax = 100, p0 = 30) 

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Função para simular e plotar a série temporal do preço

ts_plot_price <- function(model, p0, ts_length = 20) {
  p <- numeric(ts_length)
  p[1] <- p0
  for (t in 2:ts_length) {
    p[t] <- g(model, p[t - 1])
  }
  df_ts <- data.frame(t = 1:ts_length, p = p)
  
# Adicionando uma linha horizontal no preço de equilíbrio (c) para referência
ggplot(df_ts, aes(x = t, y = p)) +
      geom_line(color = "blue") +
      geom_point(color = "darkblue") +
      geom_hline(yintercept = model$c, linetype="dashed", color="gray") +
      labs(x = "Período (t)", y = expression(p[t])) +
      theme_minimal()
}

# Simulação 1: Preço inicial (P0 = 35) ABAIXO do equilíbrio
ts_plot_price(market_milho, p0 = 35, ts_length = 50)

# Simulação 2: Preço inicial (P0 = 70) ACIMA do equilíbrio
ts_plot_price(market_milho, p0 = 70, ts_length = 20)

Função adaptada para calcular o preço atual (P_t) com base no preço esperado (P_e, t)

find_next_price_adaptive <- function(model, pe_t) {
  # 1. Oferta (Q_s) é baseada na Expectativa (P_e, t)
  # Q_s = d + (1 / beta) * (P_e, t - c)
  current_supply <- model$d + (1 / model$beta) * (pe_t - model$c)
  
  # 2. Preço de mercado (P_t) é definido pela Demanda (a - b * P_t = Q_s).
  # P_t = (a - Q_s) / b
  next_price <- (model$a - current_supply) / model$b
  return(next_price)
}

Função para simular e plotar a série temporal ADAPTATIVA

ts_plot_price_adaptive_ggplot <- function(model, p0, ts_length = 30, alpha = c(1.0, 0.9, 0.75)) 
  {all_results <- data.frame()
  
   for (a in alpha) {
    pe_last <- p0 
    p_values <- numeric(ts_length)
    p_values[1] <- find_next_price_adaptive(model, pe_last) 
    pe_last <- a * p_values[1] + (1 - a) * pe_last
    
    for (i in 2:ts_length) {
      p_values[i] <- find_next_price_adaptive(model, pe_last)
      pe_last <- a * p_values[i] + (1 - a) * pe_last
    }
  
    df_temp <- data.frame(
      t = 1:ts_length,
      p = p_values,
      alpha_label = paste("alpha =", a)
    )
    all_results <- rbind(all_results, df_temp)
  }
  
  # Plotar com ggplot2 e facet_wrap
  ggplot(all_results, aes(x = t, y = p)) +
    geom_line(color = "blue") +
    geom_point(color = "darkblue") +
    # Adicionando linha de equilíbrio para referência
    geom_hline(yintercept = model$c, linetype = "dashed", color = "gray") +
    facet_wrap(~ alpha_label) +
    labs(
      title = paste("Série Temporal Adaptativa para P0 =", p0),
      x = "Período (t)", 
      y = expression(p[t])
    ) +
    theme_minimal()
}

ts_plot_price_adaptive_ggplot(market_milho, p0 = 100, ts_length = 30)