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(readxl)
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(ggplot2)
library(scales)

file_path <- "CEPEA.xlsx"

# 1) Ler a primeira aba
Preco_milho_raw <- read_excel(file_path, sheet = 1)
## New names:
## • `` -> `...2`
## • `` -> `...3`
# 2) Renomear primeiras colunas para padrão Data / Preço
names(Preco_milho_raw)[1:2] <- c("date", "price")

# 3) Converter datas (dd/mm/aaaa) e preço para numérico
Preco_milho <- Preco_milho_raw %>%
  mutate(
    date  = as.Date(date, format = "%d/%m/%Y"),
    price = as.numeric(price)
  )
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `price = as.numeric(price)`.
## Caused by warning:
## ! NAs introduzidos por coerção
# 4) Filtrar desde 2005 e ordenar corretamente
Preco_milho <- Preco_milho %>%
  filter(date >= as.Date("2005-01-01")) %>%
  arrange(date)

# 5) Criar vetores para usar no Modelo Cobweb
corn_prices <- Preco_milho$price
corn_dates  <- Preco_milho$date

# 6) Plotar série temporal do milho
ggplot(Preco_milho, aes(x = date, y = price)) +
  geom_line(color = "darkblue", size = 1) +
  geom_point(color = "blue", size = 0.4) +
   scale_x_date(
    date_breaks = "2 years",
    date_labels = "%Y"
  ) +
  scale_y_continuous(
    breaks = seq(floor(min(Preco_milho$price)/10)*10,
                 ceiling(max(Preco_milho$price)/10)*10,
                 by = 10)
  ) +
  labs(
    title = "Preço do Milho (CEPEA) – Série Mensal",
    x = "Ano",
    y = "Preço (R$/sc)"
  ) +
  theme_minimal()
## 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.

###############################################################
###  COBWEB COMPLETO – Oferta REAL (CONAB) + Demanda TEÓRICA
###  Um único chunk – pronto para rodar no RMarkdown
###############################################################

library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)

### =======================
### 1) PREÇO CEPEA (LIMPO)
### =======================

Preco_raw <- read_excel("CEPEA.xlsx", col_names = FALSE)
## New names:
## • `` -> `...1`
## • `` -> `...2`
## • `` -> `...3`
Preco_milho <- Preco_raw %>%
  slice(-1) %>%
  rename(date = ...1, price = ...2) %>%
  mutate(
    price = gsub("R\\$", "", price),
    price = gsub(",", ".", price),
    price = trimws(price),
    price = as.numeric(price),
    date  = as.Date(date, format = "%d/%m/%Y"),
    ano   = as.numeric(format(date, "%Y"))
  ) %>%
  filter(!is.na(price), ano >= 2015)
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `price = as.numeric(price)`.
## Caused by warning:
## ! NAs introduzidos por coerção
Preco_anual <- Preco_milho %>%
  group_by(ano) %>%
  summarise(preco_medio = mean(price))


### ===========================
### 2) PRODUÇÃO CONAB (mil t)
### ===========================

Prod_raw <- read_excel("MilhoTotalSerieHist (1).xlsx", sheet = "Produção")
Prod_BR <- Prod_raw[1, ]

Producao <- Prod_BR %>%
  pivot_longer(cols = -Anos,
               names_to = "ano",
               values_to = "producao") %>%
  mutate(
    ano = as.numeric(ano),
    producao = as.numeric(producao)   # mil toneladas
  )


### ===========================================
### 3) MERGE CEPEA + CONAB + NORMALIZAÇÃO
### ===========================================

Base <- Producao %>%
  inner_join(Preco_anual, by = "ano") %>%
  arrange(ano) %>%
  mutate(
    preco_n = preco_medio / 100,   # normalizado (0.3–1)
    prod_n  = producao / 1000,     # mil t → milhare
    preco_lag_n = lag(preco_n)
  )


### ===============================================
### 4) OFERTA REAL ESTIMADA (via regressão lm)
### ===============================================

modelo_oferta <- lm(prod_n ~ preco_lag_n, data = Base)
c <- coef(modelo_oferta)[1]
d <- coef(modelo_oferta)[2]
if (d < 0) d <- abs(d)   # força oferta crescente


### ===============================================
### 5) DEMANDA TEÓRICA NEGATIVA (para cruza oferta)
### ===============================================

A <- 120       # intercepto – ponto alto da produção normalizada
B <- 50        # inclinação negativa moderada

# Qd = A - B*p


### ==========================
### 6) Criar o Market()
### ==========================

Market <- function(A, B, c, d){
  list(
    demand = function(p) A - B*p,
    supply = function(p) c + d*p
  )
}

market_milho <- Market(A, B, c, d)


### ============================================
### 7) CURVAS DE OFERTA E DEMANDA (COBWEB BASE)
### ============================================

p_grid <- seq(min(Base$preco_n)*0.7,
              max(Base$preco_n)*1.3,
              length.out = 200)

df_market <- data.frame(
  p  = p_grid,
  Qd = market_milho$demand(p_grid),
  Qs = market_milho$supply(p_grid)
)

ggplot(df_market, aes(x=p)) +
  geom_line(aes(y=Qd, color="Demanda"), size=1.2) +
  geom_line(aes(y=Qs, color="Oferta"), size=1.2) +
  scale_color_manual(values=c("Demanda"="blue","Oferta"="red")) +
  labs(title="Oferta REAL (CONAB) + Demanda TEÓRICA (negativa)",
       x="Preço normalizado (p/100)",
       y="Quantidade normalizada (mil t/1000)") +
  theme_minimal()

### ===============================
### 8) Função g(p) — dinâmica
### ===============================

g <- function(model, p){
  -(model$supply(p) - model$demand(p)) / B
}


### ===============================
### 9) Diagrama 45 graus
### ===============================

plot45 <- function(model, pmin, pmax){
  p_grid <- seq(pmin, pmax, length.out = 200)
  g_vals <- sapply(p_grid, function(p) g(model, p))

  df <- data.frame(p = p_grid, g = g_vals)

  ggplot(df, aes(x=p)) +
    geom_line(aes(y=g), color="darkgreen", size=1.2) +
    geom_abline(slope=1, intercept=0, linetype="dashed") +
    labs(title="Diagrama 45° – Milho (Cobweb)",
         x=expression(p[t]), y=expression(p[t+1])) +
    theme_minimal()
}

plot45(market_milho,
       min(Base$preco_n)*0.7,
       max(Base$preco_n)*1.2)

### ====================================
### 10) Simulação temporal do cobweb
### ====================================

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 <- data.frame(t=1:ts_length, p=p)

  ggplot(df, aes(x=t, y=p)) +
    geom_line(color="blue") +
    geom_point(color="darkblue") +
    labs(title="Dinâmica do Preço – Cobweb Milho",
         x="t", y=expression(p[t])) +
    theme_minimal()
}

ts_plot_price(market_milho,
              p0 = Base$preco_n[nrow(Base)],
              ts_length = 15)

### ==========================================
### 11) Simulação Adaptativa (igual Sargent)
### ==========================================

ts_price_plot_adaptive <- function(model, p0,
                                   ts_length = 20,
                                   alpha = c(1.0, 0.8, 0.6)){

  par(mfrow=c(1, length(alpha)))

  for(a in alpha){
    pe <- p0
    p_vals <- numeric(ts_length)
    p_vals[1] <- p0

    for(i in 2:ts_length){
      p_vals[i] <- g(model, pe)
      pe <- a*p_vals[i] + (1-a)*pe
    }

    plot(p_vals, type="o",
         main=bquote(alpha == .(a)),
         xlab="t", ylab="preço")
  }
}

ts_price_plot_adaptive(market_milho,
                       p0 = Base$preco_n[nrow(Base)],
                       ts_length = 20)

Parte 2 — Produção do Milho (CONAB) + Normalização

 # Ler produção total do milho (CONAB)

Prod_raw <- read_excel("MilhoTotalSerieHist (1).xlsx", sheet = "Produção")

# A primeira linha contém a série histórica do Brasil

Prod_BR <- Prod_raw[1, ]

# Transformar em formato longo

Producao_milho <- Prod_BR %>%
pivot_longer(cols = -Anos,
names_to = "ano",
values_to = "producao") %>%
mutate(
ano      = as.numeric(ano),
producao = as.numeric(producao)       # mil toneladas
)

### MERGE de CEPEA anual com CONAB anual

Preco_anual <- Preco_milho %>%
mutate(ano = as.numeric(format(date, "%Y"))) %>%
group_by(ano) %>%
summarise(preco_medio = mean(price))

Base <- Producao_milho %>%
inner_join(Preco_anual, by = "ano") %>%
arrange(ano) %>%
mutate(
preco_n = preco_medio / 100,   # Preço normalizado (30–100 → 0.3–1)
prod_n  = producao / 1000,     # Produção normalizada (mil t → 1000)
preco_lag_n = lag(preco_n)     # Lag para oferta
)

Parte 3 — Oferta REAL do milho (via CONAB)

modelo_oferta <- lm(prod_n ~ preco_lag_n, data = Base)

c <- coef(modelo_oferta)[1]
d <- coef(modelo_oferta)[2]

# Garante inclinação positiva da oferta

if (d < 0) d <- abs(d)
c; d
## (Intercept) 
##    75.53868
## preco_lag_n 
##    64.25015

Parte 4 — Demanda TEÓRICA (negativa e cruzando a oferta)

pmax <- max(Base$preco_n, na.rm = TRUE)

a <- max(Base$prod_n) * 1.10
b <- (a - min(Base$prod_n)*0.85) / pmax  # garante declive negativo

a; b
## [1] 155.2046
## [1] 107.4253

Parte 5 — Construção do mercado (Market)

Market <- function(a, b, c, d){
list(
a = a, b = b, c = c, d = d,
demand = function(p) a - b*p,
supply = function(p) c + d*p
)
}

market_milho <- Market(a, b, c, d)

# Grid de preços normalizados

p_grid <- seq(min(Base$preco_n)*0.8,
max(Base$preco_n)*1.2,
length.out = 300)

df_curves <- data.frame(
p = p_grid,
demanda = market_milho$demand(p_grid),
oferta  = market_milho$supply(p_grid)
)

ggplot(df_curves, aes(x=p)) +
geom_line(aes(y=demanda, color="Demanda"), size=1.2) +
geom_line(aes(y=oferta,  color="Oferta"),  size=1.2) +
scale_color_manual(values=c("Demanda"="blue","Oferta"="red")) +
labs(title="Curvas de Oferta e Demanda – Milho",
x="Preço",
y="Quantidade") +
theme_minimal()