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 = FALSE,biasadj=FALSE) plot(ARMA_POP,main="Lincoln County Population Forecast",xlab="Year",ylab="Population") #####City level #See data http://eadiv.state.wy.us/pop/ KEM <- c(3273,3523,3688,3667,3626,3637,3611,3388,3156,3040,3020,3029,2989,2959,2976,2963,2910,2807,2729,2690,2657,2608,2575,2561,2574,2579,2603,2640,2679,2692,2642,2597,2551,2575,2578,2554,2544,2499,2457,2435,2413,2445,2445,2404,2378) DIAMOND <- c(1000,1070,1114,1101,1082,1078,1063,991,916,876,864,863,847,835,835,827,808,774,748,732,705,695,690,689,695,700,710,723,738,745,731,704,677,667,652,629,613,586,559,540,523,526,527,521,517) AREA_POP <- KEM+DIAMOND LN <- c(12177,13254,14031,14110,14111,14319,14384,13658,12875,12552,12625,12975,13124,13329,13759,14073,14206,14099,14114,14338,14621,14697,14858,15117,15539,15917,16429,17013,17629,18082,18083,17946,17822,18148,18346,18473,18766,18899,19042,19379,19658,20174,20690,20909,21000) NO_CITY <- c(10095,10392,10747,10944,11043) YEAR <- 1980:2024 ####Old data addtion:Period Ends in 1970 #See in part http://eadiv.state.wy.us/demog_data/cntycity_hist.htm LN_OLD <- c(12487,10894,10286,9023,9018,8640) KEM_OLD <- c(843,1517,1884,2026,1667,2028,2292) DIAMOND_OLD <- c(696,726,812,586,415,398,485) AREA_OLD <- KEM_OLD+DIAMOND_OLD AREA2 <- c(AREA_OLD,AREA_POP) LN2 <- c(NA,LN_OLD,LN) YEAR2 <- c(seq(1910,1980,by=10),1981:2024) A <- cbind(YEAR2,LN2) %>% as_tibble %>% rename(Population=LN2) %>% mutate(Region='Lincoln County') B <- cbind(YEAR2,AREA2) %>% as_tibble %>% rename(Population=AREA2) %>% mutate(Region='Kemmerer & Diamondvile') DATA <- rbind(A,B) %>% rename(Year=YEAR2) ggplot(aes(x=Year,y=Population,group=Region,color=Region),data=DATA2) +geom_line(linewidth=1.5) ###Kemmerer ARMA KEM_TS <- DATA %>% filter(Year>=1980,Region=='Kemmerer & Diamondvile') %>% pull(Population) %>% ts(start=c(1980),end=c(2024),frequency=1) BC <- BoxCox.lambda(KEM_TS) MODEL2 <- auto.arima(KEM_TS, lambda = BC) ARMA_KEM <- forecast(MODEL2,fan=TRUE,h=40,biasadj=FALSE) plot(ARMA_KEM,main="Kemmerer & Diamondvile Population",xlab="Year",ylab="Population",ylim=c(0,6000)) ###IMPLAN adjustment #TerraPower: Commuting SAM LOW_EMP <- 310.75 #Employees #TerraPower: No commuting HIGH_EMP <- 325 #Employees #Employment in Licolin LIN_EMP <- 625 #Regional (zip code level) ratios EMPLOYMENT <- 2425.84 HOUSEHOLDS <- 1850.20 POPULATION <- 4038 #Half of Lincolin Employment in some low industried can be added to Kemmerer effect. For example restruants may be scaled based on current output, when new ones will be built. ADJ_LOCAL_EMP <- 78.15 #85% commute from OUTSIDE Kemmer #15% local labor RATIO <- POPULATION/EMPLOYMENT #Note that number in household will be much larger. This acounts for average commuting rates by looking at population compared to employment LOW_GROWTH <- RATIO*LOW_EMP RATIO*HIGH_EMP MID_GROWTH <- RATIO*(LOW_EMP+ADJ_LOCAL_EMP) HIGH_GROWTH <- RATIO*(HIGH_EMP+ADJ_LOCAL_EMP) MID_GROWTH HIGH_GROWTH ############## 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