library(tidyverse)
library(fixest)

NAMES <- c("Trumbull Energy Center","Orange County Advanced Power Station","Shady Hills Combined Cycle Facility","Cumberland (TN)","Kingston","Lincoln Land Energy Center","Delta Blues Advanced Power Station","Homer","Vicksburg","Cheyenne Prairie","Legend and Lone Star","Wolf Summit Energy","Franklin Farms","Viola Generating Station","CPV Basin Ranch Energy Center","Smarr Combined Cycle Energy Facility","Waterford 5 & 6","Lake Charles","Jefferson Power Station")
CAP <- c(308.7*2+360,
	 453*2+400,
	 612,
	 454.8+323,
	323,
	 638.4*2,
	 477+341.7,
	 4678,
	 819.4,
	 132,
	754+453,
	600,
	1560,
	705,
	1490,
	1425,
	1640,
	994,
	754
)

DETAILED_REPORT <- c(1,1,1,1,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0)
EMP <- c(25,27,12.5,30,NA,34,21,250,21,12,73,30,44,30,45,30,30,30,22)
PEAK_CONST <- c(920,NA,600,NA,300,500,300,2500,560,400,NA,400,(1500+1800)/2,500,1000,1200,700,1100,NA)
DATE <- c(2017,2018,2018,2022,2015,2020,2024,2025,2025,2012,2024,2025,2024,2024,2026,2026,2025,2020,2025)
DATA <- cbind(CAP,EMP,PEAK_CONST,DETAILED_REPORT,DATE)  %>% as.matrix %>% as_tibble 
DATA$NAMES <- NAMES
DATA
colnames(DATA ) <- c("Capacity","Emp","Peak_Const","Report_Quality","Report_Year","Plant")
REG_DATA <-DATA %>% select(Emp,Capacity,Report_Quality,Report_Year)
REG_DATA_TRIMMED <-DATA %>% filter(Plant!='Homer') %>% select(Emp,Capacity,Report_Quality,Report_Year)
DATA
MOD1 <- feols(log(Emp)~Capacity+Report_Quality+Report_Year,REG_DATA,vcov='hetero')
feols(log(Emp)~Capacity+Report_Quality+Report_Year,REG_DATA,vcov='hetero')
PRED_DATA <- DATA[1:4,]
PRED_DATA$Capacity <- c(250,750,1000,1500)
PRED_DATA$Report_Year <- 2026
PRED_DATA$Plant <- 'Enyo'
PRED_DATA$Emp<- NA 
round(exp(predict(MOD1,newdata=PRED_DATA)))
#PRED_DATA$Capacity<- c(640,1280,1280,3680) #Larger number includes the single cyle turbines that look to be turned off most of the time
PRED_OPERATING <- round(exp(predict(MOD1,newdata=PRED_DATA)))
PRED_OPERATING 

ggplot(REG_DATA,aes(x=Capacity,y=log(Emp)))+geom_point()+geom_smooth(method='lm')
#feols(log(Emp)~Capacity+Report_Quality+Report_Year,REG_DATA_TRIMMED,vcov='hetero')


MOD2 <- feols(Peak_Const~Capacity+Report_Year+Report_Quality,DATA,vcov='hetero')
feols(Peak_Const~Capacity+Report_Year+Report_Quality,DATA,vcov='hetero')
ggplot(DATA,aes(x=Capacity,y=Peak_Const))+geom_point()+geom_smooth(method='lm')

###Construction Data
CONST_SCHEDULE <- c(45,75,90,115,175,235,285,350,400,350,255,175,95,65)
PEAK_RATIO <- (sum(CONST_SCHEDULE )/12)/400
PEAK_RATIO 
PRED_DATA_CONST <- PRED_DATA
PRED_DATA_CONST$Capacity <- c(640,640,0,2400)
CONST_EMPLOYMENT_PRED <- round(predict(MOD2,PRED_DATA_CONST)*PEAK_RATIO)
CONST_EMPLOYMENT_PRED 

MOD2
CONST_SCHEDULE 
95*2


PRED_DATA <- DATA[1:4,]
PRED_DATA