library(tidyverse)
library(forecast)
library(strucchange)
library(vars)
library(rugarch)

Part I Introduction

Our analysis focuses primarily on rental vacancy rates in the United States provided by FRED. Simply, the rental vacancy rate is the proportion of rental properties that are currently vacant in the United States. We also utilize the unemployment rate in the United States in an attempt to improve our model. We believe that higher rates of unemployment could impact vacancy rates as some individuals begin to be unable to pay rent. We expect these dynamics to become even more apparent through shocks during recessions as both vacancy and unemployment rates increase.

Part II

# Reading the unemployment rate data and the rental vacancy rate data in the workspace
unrate <- read_csv('unrate.csv')
vacancy <- read_csv('rvr.csv')

# Converting the data into time series 
ts_unrate <- ts(unrate$UNRATE, start = c(2000,1), frequency = 4) # Quarterly Data
ts_vacancy <- ts(vacancy$RRVRUSQ156N, start = c(2000,1), frequency = 4) # Quarterly Data

unrate_train <- ts(unrate$UNRATE, start = c(2000,1), end = c(2016,4), frequency = 4)
unrate_test <- ts(unrate$UNRATE %>% tail(12), start = c(2017,1), frequency = 4)

vacancy_train <- ts(vacancy$RRVRUSQ156N, start = c(2000,1), end = c(2016,4), frequency = 4)
vacancy_test <- ts(vacancy$RRVRUSQ156N %>% tail(12), start = c(2017,1), frequency = 4)

II.1.(a)

# Producing time series plots & ACF, PACF plots of the data
# 1. U.S. Rental Vacancy Rate
plot(ts_vacancy, title="U.S. Rental Vacancy Rate", ylab="Rental Vacancy Rate") # Time Series Plot

acf(ts_vacancy, main="U.S. Rental Vacancy Rate") # ACF Plot

pacf(ts_vacancy, main="U.S. Rental Vacancy Rate") # PACF Plot

# 2. U.S. Unemployment Rate
plot(ts_unrate, title="U.S. Unemployment Rate", ylab="Unemployment Rate") # Time Series Plot

acf(ts_unrate, main="U.S. Unemployment Rate", lag.max=20) # ACF Plot

pacf(ts_unrate, main="U.S. Unemployment Rate", lag.max=20) # PACF Plot

II.1.(b)

# Fitting model with trend, seasonality, and cyclical components
# 1. U.S. Rental Vacancy Rate
vac_decomp <- decompose(vacancy_train) # Decompose time series components for vacancy rate.
dummies <- seasonaldummy(vac_decomp$seasonal) # Assign seasonal component to dummies
season_mod <- lm(vac_decomp$seasonal ~ dummies) 
arima_mod <- Arima(vacancy_train, order = c(1,2,1)) # Arima for trend and cycles. 

model1 <- tslm(vacancy_train ~ arima_mod$fitted + season_mod$fitted.values) # Combining values from multiple models.

plot(vacancy_train, main = "Rental Vacancy Rate & Model Fit", xlab = "Year", ylab = "Rental Vacancy Rate")
lines(model1$fitted, col = "red")
legend("topright",legend = c("Rental Vacancy Rate", "Fit"), text.col = 1:2, bty = "n")

# 2. U.S. Unemployment Rate
unrate_decomp <- decompose(unrate_train) # Decompose time series components for unemployment rate.
dummies2 <- seasonaldummy(unrate_decomp$seasonal) # Assign seasonal component to dummies
season_mod2 <- lm(unrate_decomp$seasonal ~ dummies2)
arima_mod2 <- Arima(unrate_train, order = c(1,1,0)) # Arima for trend and cycles.

model2 <- tslm(unrate_train ~ arima_mod2$fitted + season_mod2$fitted.values) # Combining values from multiple models. 

plot(unrate_train, main = "Unemployment Rate & Model Fit", xlab = "Year", ylab = "Unemployment Rate")
lines(model2$fitted, col = "red")
legend("topright",legend = c("Unemployment Rate", "Fit"), text.col = 1:2, bty = "n")

Model 1 does decently well capturing the dynamics of the rental vacancy rate time series. We captured seasonal components of the time series and incorporated it with an arima model based on acf, pacf, and other analytical diagnostics. We construced model 2 using the same tactics. The global financial crisis appears to cause some issues with both models, but the lagged aspect to our model captures the rapid increasee in the unemployment rate a little bit late.

II.1.(c)

# Plotting the respective residuals vs fitted values
# 1. U.S. Rental Vacancy Rate
plot(model1$fitted, model1$residuals, pch=20, col="lightblue", lwd=1, xlab="Fitted Values", ylab="Residuals")
abline(h=0, lwd=2, col="saddlebrown")

# 2. U.S. Unemployment Rate
plot(model2$fitted, model2$residuals, pch=20, col="orange", lwd=1, xlab="Fitted Values", ylab="Residuals")
abline(h=0, lwd=2, col="saddlebrown")

When viewing the vacancy rate residuals and fitted values there appears to be a lot of bunch around certain values (7,8 and 10) of the vacancy rate. This could be the result of persistence around conditional means at various points throughout the time series. Residuals are small except for a few outliers particularly where the model overestimates. Lots of grouping around the lower unemployment rates which probably results from characteristics also seen in the vacancy rate. Residuals are small and random aside from a major outlier in 2009 which is likely due to dynamics from the global financial crisis.

II.1.(e)

# Plotting the ACF and PACF of the respective residuals
# 1. U.S. Rental Vacancy Rate
acf(model1$residuals, main="U.S. Rental Vacancy Rate - Respective residuals from the model", lag.max=500)

pacf(model1$residuals, main="U.S. Rental Vacancy Rate - Respective residuals from the model", lag.max=500)

# 2. U.S. Unemployment Rate
acf(model2$residuals, main="U.S. Unemployment Rate - Respective residuals from the model", lag.max=500)

pacf(model2$residuals, main="U.S. Unemployment Rate - Respective residuals from the model", lag.max=500)

ACF and PACF of residuals appear random with only a couple of statistically significant spikes at very high lags.

II.1.(f)

# Plotting the respective CUSUM for each model
# 1. U.S. Rental Vacancy Rate
plot(efp(model1$residuals ~ 1, type="Rec-CUSUM"))

# 2. U.S. Unemployment Rate
plot(efp(model2$residuals ~ 1, type="Rec-CUSUM"))

CUSUM plot for the vacancy rate stays within the bands throughout the time series and also stay close to zero. Overall, model does tend to overestimate, so CUSUM tends to stay above zero, but eventually returns to zero in most cases. Unemployment model does not perform as well, but trend back toward zero before errors are significant. Large spikes above zero indicate the model dealing with outliers during recessions.

II.1.(g)

# Plotting the respective Recursive Residuals for each model
# 1. U.S. Rental Vacancy Rate
y1 <- recresid(model1$residuals ~ 1)
plot(y1, pch=16, main="Recursive Residuals - U.S. Rental Vacancy Rate", ylab="Recursive Residuals", col="lightblue", type = "l")

# 2. U.S. Unemployment Rate
y2 <- recresid(model2$residuals ~ 1)
plot(y2, pch=16, main="Recursive Residuals - U.S. Unemployment Rate", ylab="Recursive Residuals", col="orange", type = "l")

Recursive residuals for the rental vacancy rate appear random but exhibit persistence in areas. The recursive residuals for the unemployment rate generally stay around zero, but a large positive spike appears during the global financial crisis. Overall, there appears to be a lot of persistence below zero.

II.1.(h)

# Generating the associated diagnostic statistics for each model
# 1. U.S. Rental Vacancy Rate
summary(model1) # Summary statisitcs of model.
## 
## Call:
## tslm(formula = vacancy_train ~ arima_mod$fitted + season_mod$fitted.values)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83843 -0.14341 -0.02557  0.20178  0.72576 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               0.56718    0.26935   2.106  0.03910 *  
## arima_mod$fitted          0.93359    0.02947  31.680  < 2e-16 ***
## season_mod$fitted.values  1.33238    0.40475   3.292  0.00161 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2916 on 65 degrees of freedom
## Multiple R-squared:  0.9395, Adjusted R-squared:  0.9376 
## F-statistic: 504.6 on 2 and 65 DF,  p-value: < 2.2e-16
accuracy(model1) # Accuracy measures of model.
##                         ME      RMSE     MAE        MPE     MAPE      MASE
## Training set -3.917816e-17 0.2850547 0.21935 -0.1021463 2.427132 0.4373334
##                   ACF1
## Training set 0.1281688
# 2. U.S. Unemployment Rate
summary(model2)
## 
## Call:
## tslm(formula = unrate_train ~ arima_mod2$fitted + season_mod2$fitted.values)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.70414 -0.22746 -0.04274  0.18141  1.14861 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                0.50550    0.19883   2.542   0.0134 *  
## arima_mod2$fitted          0.90126    0.03892  23.154   <2e-16 ***
## season_mod2$fitted.values  1.27717    3.54907   0.360   0.7201    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3575 on 65 degrees of freedom
## Multiple R-squared:  0.8919, Adjusted R-squared:  0.8885 
## F-statistic: 268.1 on 2 and 65 DF,  p-value: < 2.2e-16
accuracy(model2)
##                         ME      RMSE       MAE       MPE     MAPE      MASE
## Training set -6.544675e-18 0.3494868 0.2630632 -0.529065 5.297892 0.3380732
##                    ACF1
## Training set 0.09439724

Both models perform well in terms of r-squared and signifigance of coefficients. RMSE and MAPE also look decent for fit on training data, but better comparison will come with the introduction of more models later in the report.

II.1.(i)

# Plotting 12-step ahead forecasts of the models
# 1. U.S. Rental Vacancy Rate
plot(forecast(model1$fitted.values, h=12), main = "12-Step Forecast of U.S. Rental Vacancy Rate", ylab = "Vacancy Rate", xlab = "Year") # Overall looking at a glance

plot(forecast(model1$fitted.values, h=12), xlim=c(2015, 2023), main = "12-Step Forecast of U.S. Rental Vacancy Rate", ylab = "Vacancy Rate", xlab = "Year") # A closer look

# 1. U.S. Unemployment Rate
plot(forecast(model2$fitted.values, h=12), main = "12-Step Forecast of U.S. Unemployment Rate", ylab = "Unemployment Rate", xlab = "Year") # Overall looking at a glance

plot(forecast(model2$fitted.values, h=12), xlim=c(2015, 2021), main = "12-Step Forecast of U.S. Unemployment Rate", ylab = "Unemployment Rate", xlab = "Year") # A closer look

II.1.(j)

# US vacancy rates

# Generate ARIMA model automatically
vacancy_arima <- auto.arima(vacancy_train)

# Generate ETS model
vacancy_ets <- ets(vacancy_train)

# Generate Holt-Winters model
vacancy_hw <- hw(vacancy_train, seasonal = "multiplicative", h = 12)

# Forecast every model 12 steps ahead
vacancy_arima_fcast <- vacancy_arima %>% forecast(12)
vacancy_ets_fcast <- vacancy_ets %>% forecast(12)
vacancy_hw_fcast <- vacancy_hw %>% forecast(12)
vacancy_custom_fcast <- model1$fitted.values %>% forecast(12)

# Output a table of MAPE values and the respective model
c("ARIMA" = accuracy(vacancy_arima_fcast,ts_vacancy)["Test set","MAPE"],
"ETS" =  accuracy(vacancy_ets_fcast, ts_vacancy)["Test set","MAPE"],
"Holtz-Winters" = accuracy(vacancy_hw_fcast,ts_vacancy)["Test set","MAPE"],
"Custom" = accuracy(vacancy_custom_fcast, ts_vacancy)["Test set","MAPE"])
##         ARIMA           ETS Holtz-Winters        Custom 
##      8.984562      8.176662      8.768658     12.472394
# US unemployment rate

# Generate ARIMA model automatically
unrate_arima <- auto.arima(unrate_train)

# Generate ets model
unrate_ets <- ets(unrate_train)

# Generate Holt-Winters model
unrate_hw <- hw(unrate_train, seasonal = "multiplicative", h = 12)

# Forecast every model 12 steps ahead
unrate_arima_fcast <- unrate_arima %>% forecast(12)
unrate_ets_fcast <- unrate_ets %>% forecast(12)
unrate_hw_fcast <- unrate_hw %>% forecast(12)
unrate_custom_fcast <- model2$fitted.values %>% forecast(12)

# Output a table of MAPE values and the respective model
c("ARIMA" = accuracy(unrate_arima_fcast,ts_unrate)["Test set","MAPE"],
"ETS" =  accuracy(unrate_ets_fcast, ts_unrate)["Test set","MAPE"],
"Holtz-Winters" = accuracy(unrate_hw_fcast,ts_unrate)["Test set","MAPE"],
"Custom" = accuracy(vacancy_custom_fcast, ts_vacancy)["Test set","MAPE"])
##         ARIMA           ETS Holtz-Winters        Custom 
##      18.43015      19.01909      19.11167      12.47239

ETS performs best in terms of MAPE for the vacancy rate while ETS captures the unemployment rate data the best.In both cases ETS outperforms the other models but the difference is minimal.

II.1.(k)

# I dont know how we're supposed to weight
# Obtain weights using OLS weighting
vacancy_weight <- lm(vacancy_train ~ vacancy_arima$fitted + vacancy_ets$fitted + vacancy_hw$fitted + model1$fitted.values)[["coefficients"]]

vacancy_weight
##          (Intercept) vacancy_arima$fitted   vacancy_ets$fitted 
##           0.57795661           0.04990273           0.94633081 
##    vacancy_hw$fitted model1$fitted.values 
##           0.06838757          -0.13317020
# Get weighted forecast values
vacancy_combo <- (vacancy_arima_fcast[["mean"]] * vacancy_weight[2] 
                 + vacancy_ets_fcast[["mean"]] * vacancy_weight[3]
                 + vacancy_hw_fcast[["mean"]] * vacancy_weight[4]
                 + vacancy_custom_fcast[["mean"]] * vacancy_weight[5]
                 + vacancy_weight[1])

# Obtain accuracy statistics on the combined forecast model
accuracy(vacancy_combo, ts_vacancy)
##                 ME      RMSE       MAE      MPE     MAPE      ACF1 Theil's U
## Test set 0.3928728 0.4394901 0.3953152 5.630235 5.665126 -0.103688  1.391698
# Obtain weights using OLS weighting
unrate_weight <- lm(unrate_train ~ unrate_arima$fitted + unrate_ets$fitted + unrate_hw$fitted + model2$fitted.values)[["coefficients"]]

unrate_weight
##          (Intercept)  unrate_arima$fitted    unrate_ets$fitted 
##            0.1962787            0.6322738           -0.2748003 
##     unrate_hw$fitted model2$fitted.values 
##            0.1468749            0.4571788
# Get weighted forecast values
unrate_combo <- (unrate_arima_fcast[["mean"]] * unrate_weight[2] 
                 + unrate_ets_fcast[["mean"]] * unrate_weight[3]
                 + unrate_hw_fcast[["mean"]] * unrate_weight[4]
                 + unrate_custom_fcast[["mean"]] * unrate_weight[5]
                 + unrate_weight[1])

# Obtain accuracy statistics on the combined forecast model
accuracy(unrate_combo, ts_unrate)
##                 ME     RMSE      MAE      MPE     MAPE      ACF1 Theil's U
## Test set 0.7728285 1.788697 1.318974 6.931207 18.18428 0.8065233  2.052805

When combining the forecasts MAPE improves a significant amount, a reduction of over 3%, for the rental vacancy rate model. The combined model for the unemployment rate does not perform as well as the best performing individual model, ETS, and in fact does not perform as well as any of the unemployment rate models.

II.1.(l)

#Creating VARS model
y=cbind(unrate_train, vacancy_train) %>%
  data.frame()
VARselect(data.frame(unrate_train, vacancy_train))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      2      2 
## 
## $criteria
##                  1            2           3          4           5           6
## AIC(n) -4.42651684 -4.624638618 -4.53597782 -4.4691771 -4.46229158 -4.37813494
## HQ(n)  -4.34349089 -4.486262037 -4.34225060 -4.2200993 -4.15786311 -4.01835583
## SC(n)  -4.21336756 -4.269389823 -4.03862950 -3.8297293 -3.68074423 -3.45448808
## FPE(n)  0.01195827  0.009815617  0.01074179  0.0115148  0.01164365  0.01274455
##                 7           8           9          10
## AIC(n) -4.2839012 -4.26330889 -4.27872581 -4.24777076
## HQ(n)  -3.8687714 -3.79282852 -3.75289481 -3.66658912
## SC(n)  -3.2181548 -3.05546299 -2.92878039 -2.75572582
## FPE(n)  0.0141242  0.01458312  0.01457302  0.01531402
y_model=VAR(y,p=6)
summary(y_model)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: unrate_train, vacancy_train 
## Deterministic variables: const 
## Sample size: 62 
## Log Likelihood: -29.532 
## Roots of the characteristic polynomial:
## 0.9732 0.9732 0.7888 0.7888 0.6681 0.6681 0.5778 0.5778 0.5745 0.5745 0.5243 0.2884
## Call:
## VAR(y = y, p = 6)
## 
## 
## Estimation results for equation unrate_train: 
## ============================================= 
## unrate_train = unrate_train.l1 + vacancy_train.l1 + unrate_train.l2 + vacancy_train.l2 + unrate_train.l3 + vacancy_train.l3 + unrate_train.l4 + vacancy_train.l4 + unrate_train.l5 + vacancy_train.l5 + unrate_train.l6 + vacancy_train.l6 + const 
## 
##                   Estimate Std. Error t value Pr(>|t|)    
## unrate_train.l1   1.432829   0.138054  10.379 5.78e-14 ***
## vacancy_train.l1  0.055350   0.146835   0.377  0.70784    
## unrate_train.l2  -0.612767   0.234984  -2.608  0.01205 *  
## vacancy_train.l2  0.005163   0.183700   0.028  0.97769    
## unrate_train.l3   0.027000   0.236021   0.114  0.90939    
## vacancy_train.l3 -0.056350   0.184375  -0.306  0.76118    
## unrate_train.l4  -0.097288   0.224681  -0.433  0.66691    
## vacancy_train.l4  0.132770   0.183309   0.724  0.47233    
## unrate_train.l5   0.180294   0.218148   0.826  0.41254    
## vacancy_train.l5 -0.139483   0.179939  -0.775  0.44197    
## unrate_train.l6  -0.056537   0.137875  -0.410  0.68355    
## vacancy_train.l6 -0.082123   0.144500  -0.568  0.57241    
## const             1.444707   0.531989   2.716  0.00911 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.3534 on 49 degrees of freedom
## Multiple R-Squared: 0.9137,  Adjusted R-squared: 0.8925 
## F-statistic: 43.23 on 12 and 49 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation vacancy_train: 
## ============================================== 
## vacancy_train = unrate_train.l1 + vacancy_train.l1 + unrate_train.l2 + vacancy_train.l2 + unrate_train.l3 + vacancy_train.l3 + unrate_train.l4 + vacancy_train.l4 + unrate_train.l5 + vacancy_train.l5 + unrate_train.l6 + vacancy_train.l6 + const 
## 
##                   Estimate Std. Error t value Pr(>|t|)    
## unrate_train.l1  -0.085882   0.132223  -0.650   0.5190    
## vacancy_train.l1  0.761171   0.140633   5.412 1.85e-06 ***
## unrate_train.l2   0.241204   0.225058   1.072   0.2891    
## vacancy_train.l2  0.317754   0.175941   1.806   0.0771 .  
## unrate_train.l3  -0.183748   0.226051  -0.813   0.4202    
## vacancy_train.l3 -0.179892   0.176587  -1.019   0.3133    
## unrate_train.l4   0.008544   0.215191   0.040   0.9685    
## vacancy_train.l4  0.186866   0.175566   1.064   0.2924    
## unrate_train.l5   0.215965   0.208934   1.034   0.3064    
## vacancy_train.l5 -0.048427   0.172338  -0.281   0.7799    
## unrate_train.l6  -0.102452   0.132051  -0.776   0.4416    
## vacancy_train.l6 -0.071268   0.138396  -0.515   0.6089    
## const            -0.183818   0.509518  -0.361   0.7198    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.3385 on 49 degrees of freedom
## Multiple R-Squared: 0.934,   Adjusted R-squared: 0.9178 
## F-statistic: 57.79 on 12 and 49 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##               unrate_train vacancy_train
## unrate_train      0.124892     -0.008914
## vacancy_train    -0.008914      0.114564
## 
## Correlation matrix of residuals:
##               unrate_train vacancy_train
## unrate_train       1.00000      -0.07452
## vacancy_train     -0.07452       1.00000
acf(residuals(y_model)[,1])

pacf(residuals(y_model)[,1])

acf(residuals(y_model)[,2])

pacf(residuals(y_model)[,2])

The VAR model with an order of 6 shows that neither variable variable have a affect on each other.

II.1.(m)

#Plotting impulse response functions
plot(irf(y_model))

There is a large effect on unemployment rate from a shock of the unemployment rate and a small effect on the vacancy rate from a shock of the unemployment rate. The increase in vacancy rate is lagged, but appears statistically significant over several periods before trending back to zero. A shock to the vacancy rate exhibits a very small decrease then persistence above zero for the duration of the graph. The effect on the unemployment rate is negligible.

II.1.(n)

#Generating Granger Test
grangertest(unrate_train ~ vacancy_train, order = 6)
## Granger causality test
## 
## Model 1: unrate_train ~ Lags(unrate_train, 1:6) + Lags(vacancy_train, 1:6)
## Model 2: unrate_train ~ Lags(unrate_train, 1:6)
##   Res.Df Df      F Pr(>F)
## 1     49                 
## 2     55 -6 0.8199 0.5599
grangertest(vacancy_train ~ unrate_train, order = 6)
## Granger causality test
## 
## Model 1: vacancy_train ~ Lags(vacancy_train, 1:6) + Lags(unrate_train, 1:6)
## Model 2: vacancy_train ~ Lags(vacancy_train, 1:6)
##   Res.Df Df      F Pr(>F)
## 1     49                 
## 2     55 -6 1.1333 0.3573

It appears that shocks to the unemployment rate can Granger-cause fluctuations in the vacancy rate at the 99% significance level.

II.1.(o)

#Predicting with VARS model
var.predict = predict(object=y_model, n.ahead=12)
plot(var.predict)

accuracy(var.predict$fcst$vacancy_train[,1], vacancy_test)
##                  ME      RMSE       MAE       MPE     MAPE      ACF1 Theil's U
## Test set -0.6804978 0.9789689 0.7852581 -10.19459 11.60932 0.6578505  3.131454
accuracy(var.predict$fcst$unrate_train[,1], unrate_test)
##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set -2.502523 2.593456 2.502523 -64.11053 64.11053 0.6781094  14.65776

Understandably, the VAR model for unemployment does not perform very well as we have seen with the reults from the IRF and granger causality. The model for the vacancy rate does much better and performs almost nearly as well the combination model in terms of MAPE.

II.1.(p)

#load favorite model residuals
vacancy_resid <- vacancy_hw$residuals

# load specification for garch model
vacancy_spec <- ugarchspec(
  variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sstd")

# Fit garch model for unemployment rate data
vacancy_garch <- ugarchfit(spec = vacancy_spec, data = vacancy_resid)

# Forecast garch model
vacancy_garch_fcast <- ugarchforecast(vacancy_garch, n.ahead = 12, n.roll = 0, out.sample = 0)

plot(vacancy_garch_fcast, which = 1)

# Obtain accuracy metrics for garch model
#accuracy(vacancy_garch_fcast,ts_vacancy)

#load favorite model residuals
unrate_resid <- unrate_hw$residuals

# load specification for garch model
unrate_spec <- ugarchspec(
  variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 1), include.mean = TRUE),
  distribution.model = "sstd")

# Fit garch model for unemployment rate data
unrate_garch <- ugarchfit(spec = unrate_spec, data = unrate_resid)

# Forecast garch model
unrate_garch_fcast <- ugarchforecast(unrate_garch, n.ahead = 12)
plot(unrate_garch_fcast, which = 1)

# Obtain accuracy metrics for garch model
#accuracy(unrate_garch_fcast,ts_unrate)

III Conclusion

Through the application of our own model, several automatic models, VAR, and GARCH we found that the best performing model in terms of MAPE for our data was the combination of our own model, ETS, and Holtz-Winters. This combined model performed several percentage points better than the best performing solo model. We were also pleasantly surprised by the effectiveness of the IRF and VAR model. When choosing our data we thought maybe the unemployment rate would be too disconnected from the rental vacancy rate, but the delayed response found in the IRF for unemployment rate on vacancy rate mimics real world dynamics. For future work we would like to examine other variables and the possible Granger-causality on rental vacancy rates in the United States and abroad. It may also be beneficial to observe and analyze a stationary form of the rental vacancy rate time series through differences or logs of the data.

IV Sources

U.S. Bureau of Labor Statistics, Unemployment Rate [UNRATE], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/UNRATE, March 5, 2020.

U.S. Census Bureau, Rental Vacancy Rate for the United States [RRVRUSQ156N], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/RRVRUSQ156N, March 4, 2020.