New Haven Economic Performance Index
Francesca Micallef
11/2/2021
This script downloads data from FRED for the Economic Performance Index. The index contains 5 series. Revised: Sys.Date() # “2020-10-27” Revised: Sys.Date() # “2020-10-10”
Download the Series by FRED code
#Use pdfetch to draw, assemble and prepare data
permits_CT = pdfetch_FRED(c("CTBPPRIVSA")) # Building Permits in CT
AWH_NH = pdfetch_FRED(c("SMU09757000500000002SA")) # Average Weekly Hours in New Haven
AWE_NH = pdfetch_FRED(c("SMU09757000500000011SA")) # Average Weekly Earnings New
EdsMeds_NH = pdfetch_FRED(c("NEWH709EDUH")) # Eds & Meds Employment New Haven
UR_NH = pdfetch_FRED(c("NEWH709UR")) # Unemployment Rate in New Haven
#adjust unemployment so increases are positive
UR_NH = max(UR_NH) - UR_NH
plot(UR_NH) 
plot(EdsMeds_NH)
plot(AWH_NH)
plot(AWE_NH)
plot(permits_CT)
#Shorten them cause UR lags 1 month extra
tail(AWH_NH); head(UR_NH)## SMU09757000500000002SA
## 2021-05-31 33.7
## 2021-06-30 33.4
## 2021-07-31 33.6
## 2021-08-31 33.8
## 2021-09-30 33.7
## 2021-10-31 33.9## NEWH709UR
## 1990-01-31 6.2
## 1990-02-28 6.2
## 1990-03-31 6.2
## 1990-04-30 6.1
## 1990-05-31 6.0
## 1990-06-30 5.9permits_CT <- permits_CT["2010-01-31/2021-08-31"] #reduce the length of all series
AWH_NH <- AWH_NH["2010-01-31/2021-08-31"] #reduce the length of all series
AWE_NH <- AWE_NH["2010-01-31/2021-08-31"] #reduce the length of all series
EdsMeds_NH <- EdsMeds_NH["2010-01-31/2021-08-31"] #reduce the length of all series
UR_NH <- UR_NH["2010-01-31/2021-08-31"] #reduce the length of all series
#Assemble into a data frame
ct_data = data.frame(permits_CT, AWH_NH, AWE_NH, EdsMeds_NH, UR_NH)
ct_data_m <- ts(ct_data,start=c(2010,1),frequency=12)
plot(ct_data_m)
lo = nrow(ct_data)Method: Indexing Rebase each series to 2010:1 = 100
index_rebased <- dplyr::mutate(ct_data,
permits_CT=100*permits_CT/as.numeric(permits_CT[1]),
AWH_NH=100*AWH_NH/as.numeric(AWH_NH[1]),
AWE_NH=100*AWE_NH/as.numeric(AWE_NH[1]),
EdsMeds_NH=100*EdsMeds_NH/as.numeric(EdsMeds_NH[1]),
UR_NH=100*UR_NH/as.numeric(UR_NH[1]))
head(index_rebased)## CTBPPRIVSA SMU09757000500000002SA SMU09757000500000011SA NEWH709EDUH
## 2010-01-31 281 33.4 836 73.5
## 2010-02-28 317 32.0 817 73.4
## 2010-03-31 292 32.5 840 73.5
## 2010-04-30 335 32.7 848 74.1
## 2010-05-31 254 32.7 860 74.2
## 2010-06-30 181 32.8 857 74.1
## NEWH709UR CTBPPRIVSA SMU09757000500000002SA SMU09757000500000011SA
## 2010-01-31 1.1 100.0 100.0 100.0
## 2010-02-28 0.9 112.7 96.0 97.7
## 2010-03-31 0.8 104.0 97.5 100.5
## 2010-04-30 0.8 119.2 98.1 101.4
## 2010-05-31 0.8 90.6 98.0 102.9
## 2010-06-30 0.8 64.5 98.4 102.5
## NEWH709EDUH NEWH709UR
## 2010-01-31 100.0 100.0
## 2010-02-28 99.9 81.8
## 2010-03-31 100.0 72.7
## 2010-04-30 100.9 72.7
## 2010-05-31 100.9 72.7
## 2010-06-30 100.9 72.7 tail(index_rebased) ## CTBPPRIVSA SMU09757000500000002SA SMU09757000500000011SA NEWH709EDUH
## 2021-03-31 450 34.3 1125 82.1
## 2021-04-30 734 33.9 1126 81.0
## 2021-05-31 311 33.7 1159 81.6
## 2021-06-30 373 33.4 1113 81.3
## 2021-07-31 284 33.6 1123 81.8
## 2021-08-31 294 33.8 1136 81.7
## NEWH709UR CTBPPRIVSA SMU09757000500000002SA SMU09757000500000011SA
## 2021-03-31 3.3 160 103 135
## 2021-04-30 3.2 261 102 135
## 2021-05-31 4.2 111 101 139
## 2021-06-30 4.2 133 100 133
## 2021-07-31 4.2 101 101 134
## 2021-08-31 4.9 105 101 136
## NEWH709EDUH NEWH709UR
## 2021-03-31 112 300
## 2021-04-30 110 291
## 2021-05-31 111 382
## 2021-06-30 111 382
## 2021-07-31 111 382
## 2021-08-31 111 445 names(index_rebased)## [1] "CTBPPRIVSA" "SMU09757000500000002SA" "SMU09757000500000011SA"
## [4] "NEWH709EDUH" "NEWH709UR" "permits_CT"
## [7] "AWH_NH" "AWE_NH" "EdsMeds_NH"
## [10] "UR_NH" #sum across columns
index = rowSums(index_rebased[6:10], dims = 1)
tail(index)## 2021-03-31 2021-04-30 2021-05-31 2021-06-30 2021-07-31 2021-08-31
## 809 899 843 859 829 899 plot.ts(index)
#rebase the index to 100
index=100*index/index[1]
tail(index)## 2021-03-31 2021-04-30 2021-05-31 2021-06-30 2021-07-31 2021-08-31
## 162 180 169 172 166 180 plot.ts(index)
Index <- ts(index, start=c(2010, 1), frequency=12) # UPDATED- October 2021!!!!!
autoplot(Index) +
theme_bw(base_size = 8) +
ylab("NHREPI = January 2010 = 100") +
labs(title="New Haven Region", subtitle = "Economic Performance") +
geom_hline(yintercept= 100) +
geom_forecast(h=6, fan = "True", level = c(90)) +
geom_line(size = 0.8, col = "darkred")
#Obtain a 6-month centered moving average
NH_MA = tsutils::cmav(Index) # this function interpolates the
# 6 NAs createdHOLTS WINTER FORECAST
hw(Index, h=6)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 2021 181 160 203 148 214
## Oct 2021 181 155 208 142 221
## Nov 2021 180 150 210 134 226
## Dec 2021 183 149 216 131 234
## Jan 2022 183 146 219 127 239
## Feb 2022 187 147 227 126 247autoplot(hw(Index, h=6),
main="New Haven Region Economic performance",
ylab = "NHEPI=100=January 2010", xlab = "Time") +
autolayer(NH_MA) +
theme(legend.position = "none")
accuracy(hw(Index)) ## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.132 15.8 10.7 -0.617 6.79 0.415 0.0312AUTO.ARIMA FORECAST
aa.model=auto.arima(Index)
autoplot(forecast(aa.model, h=6)) +
autolayer(NH_MA) +
theme(legend.position = "none")
forecast(aa.model,h=6)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 2021 173 153 194 142 205
## Oct 2021 174 148 199 135 212
## Nov 2021 174 145 202 130 217
## Dec 2021 174 143 204 126 221
## Jan 2022 173 140 206 123 223
## Feb 2022 173 138 207 120 226 accuracy(forecast(aa.model,h=6))## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.16 15.8 10 0.0986 6.36 0.39 -0.0212 autoplot(Index) +
autolayer(forecast(aa.model, h=6))+
autolayer(NH_MA) +
ylab("NHREPI = January 2010 = 100") +
labs(title="New Haven Region", subtitle = "Economic Performance") +
geom_hline(yintercept= 100) +
geom_line(size = 0.8, col = "blue")+
theme_bw(base_size = 12) +
theme(legend.position = "none")
Ensemble Forecast
fit1 <- hybridModel(Index, weights="equal")## Fitting the auto.arima model## Fitting the ets model## Fitting the thetam model## Fitting the nnetar model## Fitting the stlm model## Fitting the tbats model fc1 <- forecast(fit1, h=12) # forecast here is 12 months
autoplot(fc1) + autolayer(NH_MA)+
ggtitle("NH Region Economic Performance") +
xlab("Year") +
ylab("New Haven Economic Performance Index") +
geom_hline(yintercept= 100) +
geom_line(size = 0.8, col = "blue")+
theme_bw(base_size = 12) +
theme(legend.position = "none")
accuracy(fc1)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.06 15.6 10.1 0.0137 6.23 0.391 0.0583Different Version of Hybrid Graph
autoplot(fc1) +
autolayer(Index)+
autolayer(NH_MA)+
ylab("NHREPI = January 2010 = 100") +
xlab("Year")+
labs(title="New Haven Region", subtitle = "Economic Performance") +
geom_hline(yintercept= 100) +
geom_line(size = 0.8, col = "darkred")+
theme_bw(base_size = 12) +
theme(legend.position = "none")
Shorter Window - to highlight forecast
autoplot(fc1) +
autolayer(Index)+
autolayer(NH_MA)+
ylab("NHREPI = January 2010 = 100") +
xlab("Year")+
labs(title="New Haven Region",
subtitle = "Economic Performance and 12-M Forecast",
caption = "Created by: F. Micallef & J. Hickey, New Haven BANL") +
geom_hline(yintercept= 100) +
geom_line(size = 1.2, col = "darkred")+
theme_bw(base_size = 12) +
xlim(x = c(2019,2023)) +
theme(legend.position = "none")## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 108 row(s) containing missing values (geom_path).
## Warning: Removed 108 row(s) containing missing values (geom_path).
## Warning: Removed 108 row(s) containing missing values (geom_path).
forecast(fc1,h=6)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 2021 178 150 205 138 217
## Oct 2021 179 148 209 135 222
## Nov 2021 178 145 209 130 228
## Dec 2021 179 143 213 126 235
## Jan 2022 179 140 217 123 241
## Feb 2022 181 138 222 120 247 library(dygraphs) ## Warning: package 'dygraphs' was built under R version 4.0.5 NH_Dat = cbind(NH_MA, Index, fc1$mean) #
dygraph(NH_Dat,
main = "New Haven Economic Performance Index",
ylab = "NHREPI = January 2010 = 100" ) %>%
dySeries("fc1$mean", label = "12-MO Forecast, (%)") %>%
dySeries("NH_MA", label = "6-Mo Moving Average") %>%
dyShading(from = "2020-03-1", to = "2021-03-01",
color = "#FFE6E6") %>%
dyRangeSelector() %>%
dyLimit(min(Index[1], na.rm = TRUE), "2010:1 = 100 ",
strokePattern = "solid", color = "darkred")library(pdfetch)
suppressPackageStartupMessages(library(tidyverse))## Warning: package 'tidyverse' was built under R version 4.0.5## Warning: package 'tibble' was built under R version 4.0.5## Warning: package 'tidyr' was built under R version 4.0.5## Warning: package 'readr' was built under R version 4.0.5## Warning: package 'purrr' was built under R version 4.0.5## Warning: package 'stringr' was built under R version 4.0.5## Warning: package 'forcats' was built under R version 4.0.5 library(forecast)
library(knitr)## Warning: package 'knitr' was built under R version 4.0.5 # ================================= #
vix = rnorm(100, 98, 15)
vix.ts = ts(vix, start = c(1990,1), frequency = 12)
vix_fc = forecast(auto.arima(vix.ts))
vix_fc_table = vix_fc %>% as.data.frame() %>% tail(10)
knitr::kable(vix_fc_table, "simple",caption = "VIX ARIMA Forecast")| Point Forecast | Lo 80 | Hi 80 | Lo 95 | Hi 95 | |
|---|---|---|---|---|---|
| Jul 1999 | 91.9 | 71.0 | 113 | 59.9 | 124 |
| Aug 1999 | 92.0 | 71.1 | 113 | 60.0 | 124 |
| Sep 1999 | 92.9 | 72.0 | 114 | 60.9 | 125 |
| Oct 1999 | 92.0 | 71.1 | 113 | 60.0 | 124 |
| Nov 1999 | 92.9 | 71.9 | 114 | 60.8 | 125 |
| Dec 1999 | 92.8 | 71.8 | 114 | 60.7 | 125 |
| Jan 2000 | 92.2 | 71.2 | 113 | 60.0 | 124 |
| Feb 2000 | 92.1 | 71.0 | 113 | 59.9 | 124 |
| Mar 2000 | 93.9 | 72.8 | 115 | 61.7 | 126 |
| Apr 2000 | 92.4 | 71.3 | 114 | 60.1 | 125 |