---
title: "01_Analysis_BufferAggregate"
output: html_document
date: "2024-02-14"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Code and output for the aggregate analysis in the buffer zone and in the buffer zone/PA combined

```{r Load Packages and Data}
library(tidyverse);library(lfe)

master_agg_analysis <- read.csv("DataKeep/master_agg_analysis_buffer.csv")
master_agg_combined <- read.csv("DataKeep/master_agg_combined.csv")
master_agg_analysis$precipitation <- master_agg_analysis$precipitation*12
master_agg_combined$Precipitation <- master_agg_combined$Precipitation*12

master_agg_analysis$Tourism
master_agg_combined$Tourism

```

```{r Buffer aggregate analysis}
#Model
fem_full_formula <- as.formula("PercentDeforest ~ Tourism + PopDensity + precipitation + OilPriceAdj | Year + Name | 0 | Regions")
dd_reg_full_agg_buffer <- felm(fem_full_formula, data = master_agg_analysis)

#Output
summary_dd_reg_full_agg_buffer <- summary(dd_reg_full_agg_buffer)
agg_3km <- as.data.frame(summary_dd_reg_full_agg_buffer$coefficients)
agg_3km$Analysis <- "PA"
agg_3km$Variable <- rownames(agg_3km)
agg_3km$Location <- "Buffer 3km"
agg_3km$`Pr(>|t|)` <- round(agg_3km$`Pr(>|t|)`, digits = 3)
agg_3km$Estimate <- round(agg_3km$Estimate, digits = 3)
agg_3km$`Cluster s.e.` <- round(agg_3km$`Cluster s.e.`, digits = 3)
agg_3km$`t value` <- round(agg_3km$`t value`, digits = 3)
agg_3km$Variable[agg_3km$Variable=="PopDensity"] <- "Pop. Density"
agg_3km$Variable[agg_3km$Variable=="OilPriceAdj"] <- "Oil Price"
agg_3km$Variable[agg_3km$Variable=="Tourism1"] <- "Tourism"
agg_3km$Variable[agg_3km$Variable=="PopDensityCommune"] <- "Pop. Density"

100*agg_3km[1,1]/mean(master_agg_analysis$PercentDeforest)

write.csv(agg_3km, "Outputs/regression_table_3km_agg.csv")


```

```{r Combined (Buffer and PA) aggregate analysis}
#Model

fem_full_formula <- as.formula("PercentDeforest ~ Tourism + PopDensity + Precipitation + OilPriceAdj | Year + Name | 0 | Regions")
dd_reg_full_agg_combined <- felm(fem_full_formula, data = master_agg_combined)

#Output
summary_dd_reg_full_agg_buffer <- summary(dd_reg_full_agg_combined)
agg_combined <- as.data.frame(summary_dd_reg_full_agg_buffer$coefficients)
agg_combined$Analysis <- "PA"
agg_combined$Variable <- rownames(agg_combined)
agg_combined$Location <- "Inside PA + Buffer 3km"
agg_combined$`Pr(>|t|)` <- round(agg_combined$`Pr(>|t|)`, digits = 3)
agg_combined$Estimate <- round(agg_combined$Estimate, digits = 3)
agg_combined$`Cluster s.e.` <- round(agg_combined$`Cluster s.e.`, digits = 3)
agg_combined$`t value` <- round(agg_combined$`t value`, digits = 3)
agg_combined$Variable[agg_combined$Variable=="PopDensity"] <- "Pop. Density"
agg_combined$Variable[agg_combined$Variable=="OilPriceAdj"] <- "Oil Price"
agg_combined$Variable[agg_combined$Variable=="Tourism1"] <- "Tourism"
agg_combined$Variable[agg_combined$Variable=="PopDensityCommune"] <- "Pop. Density"

100*agg_combined[1,1]/mean(master_agg_combined$PercentDeforest)

write.csv(agg_combined, "Outputs/regression_table_combined_agg.csv")



```