library(tidyverse)
library(forecast)
library(lmtest)

source("Scripts/Functions.r")
#source("Scripts/Load_Wyoming_Web_Data.r")

DF <- FRED_GET('WYLINC3POP','LN_POP') %>% select(-YEAR)

TS <- 1000*ts(DF,start=c(1970),end=c(2024),frequency=1)
BC <- BoxCox.lambda(TS)
MODEL <- auto.arima(TS, lambda = BC)
forecast(MODEL,h=20)
ARMA_POP <- forecast(MODEL,h=35,fan=TRUE,bootstrap = TRUE,npaths=100000,biasadj=FALSE)
ARMA_POP
plot(ARMA_POP,main="Lincoln County Population Forecast",xlab="Year",ylab="Population")
help("forecast")
####Employment to pop ratio
EMP <- FRED_GET('LAUCN560230000000005','EMP') %>% inner_join(FRED_GET('WYLINC3POP','LN_POP')) %>% mutate(LN_POP=1000*LN_POP)

EMP <- EMP %>% mutate(RATIO=LN_POP/EMP)
ggplot(aes(x=YEAR,y=RATIO),data=EMP)+geom_line()
AVG_POP_RATIO <- mean(EMP$RATIO)
SD_POP_RATIO <- sd(EMP$RATIO)
##############
TS <- EMP %>% select(EMP) %>% ts(start=c(1990),end=c(2024),frequency=1)
BC <- BoxCox.lambda(TS)
MODEL2 <- auto.arima(TS, lambda = BC)
MODEL2
ARMA_EMP <- forecast(MODEL2,h=35,fan=TRUE,bootstrap = TRUE,npaths=100000,biasadj=FALSE)
plot(ARMA_EMP,main="Lincoln County Employment",xlab="Year",ylab="Employed")

##########
DATA <- FRED_GET('BPPRIV056023','PRIV_HOUSING') %>% inner_join(FRED_GET('WYLINC3POP','LN_POP')) %>% inner_join(FRED_GET('ATNHPIUS56023A','HOUSE_PRICE_INDEX')) %>% left_join(FRED_GET('LAUCN560230000000005','EMP')) %>%left_join(FRED_GET('DCOILWTICO','WTI'))  %>%left_join(FRED_GET('PCOALAUUSDM','COAL'))  %>% mutate(LN_POP=1000*LN_POP) %>% select(-YEAR) %>% ts() %>% log() %>% diff()
grangertest(PRIV_HOUSING~HOUSE_PRICE_INDEX,data=DATA,order=1)
grangertest(HOUSE_PRICE_INDEX~PRIV_HOUSING,data=DATA,order=2)
grangertest(LN_POP~PRIV_HOUSING,data=DATA,order=1)
grangertest(HOUSE_PRICE_INDEX~COAL,data=DATA,order=2) 
grangertest(LN_POP~WTI,data=DATA,order=1)
grangertest(COAL~WTI,data=DATA,order=2)
grangertest(WTI~COAL,data=DATA,order=2)
grangertest(LN_POP~COAL,data=DATA,order=1)
grangertest(EMP~LN_POP,data=DATA,order=2)


grangertest(PRIV_HOUSING~LN_POP,data=DATA,order=2) ## Signficant
grangertest(EMP~LN_POP,data=DATA,order=2) ## Signficant
grangertest(EMP~PRIV_HOUSING,data=DATA,order=2) ##Signficant
grangertest(EMP~LN_POP,data=DATA,order=2) ##Signficant
grangertest(PRIV_HOUSING~COAL,data=DATA,order=1) ##Signficant
grangertest(PRIV_HOUSING~COAL,data=DATA,order=2) ##Signficant
grangertest(PRIV_HOUSING~COAL,data=DATA,order=3) ##Signficant
grangertest(EMP~COAL,data=DATA,order=1) ##Signficant






grangertest(LN_POP~EMP,data=DATA,order=2)
grangertest(PRIV_HOUSING~EMP,data=DATA,order=2)






plot(DATA)
library(BVAR)
MOD <- bvar(DATA,lags=1, n_draw=40000)
opt_irf <- bv_irf(horizon = 25, identification = TRUE)
plot(irf(MOD,opt_irf,conf_bands = c(0.05, 0.1,0.15)),area=TRUE)



head(DATA)




DATA
# ATNHPIUS56023A Housing price index
# Median income MHIWY56023A052NCEN
# Employed Persons in Lincoln LAUCN560230000000005

#####Plan and ideas
#1) Review IMPLAN for industry multipliers
#2) Review IMPLAN for employment to population multipliers (imparted)
#3) Find a list of all planned new projects
#4) Use the IMPLAN multipliers for each sector to estimate total change
#5) Develop survey to estimate likelihood of new projects
#6) Compare to the ARMA percentile
#7) Adjust the ARMA up assuming some of these outputs are known.
#8) Occupancy rate from IMPLAN as a housing cap when projecting
#9) Model housing construciton rate (Maybe) 
#10) Employment rate by age in IMPLAN
####Other ideas, develop larger plan? Maybe look at decline in other industries as a proportion of employment

###Seperate out Kemmer and Diamondville? http://eadiv.state.wy.us/pop/wyc&sc40.htm