Upload New File

LabCode/lab_5.R

+
+# Lab 5
+
+# install packages
+install.packages('tseries')       # for time series tests
+install.packages('seastests')     # for seasonality test
+install.packages('car')           # for calculating VIF
+
+# load packages 
+library(tidyverse)
+library(ggplot2)
+library(lmtest)
+library(sandwich)
+library(stargazer)
+library(tseries)
+library(seastests)
+library(car)
+
+
+# clear 
+rm(list = ls())
+
+# import data file
+data0 <- read_csv("https://code.osu.edu/liu.6200/sp21-aede6120/-/raw/master/LabData/depos0.csv")
+
+
+
+# Time series models ----
+
+# Y ~ X
+
+# dif_Y ~ dif_X
+# Y_t - Y_(t-1)
+
+# %_Y ~ %_X
+# Y_t/Y_(t-1) - 1
+
+
+
+# Variable transformation ----
+
+# transform the dependent variable
+data0 <- data0 %>% mutate(Deposits_a1 = lag(Deposits))              # lag
+data0 <- data0 %>% mutate(Deposits_d1 = Deposits - Deposits_a1)     # difference
+data0 <- data0 %>% mutate(Deposits_p1 = Deposits/Deposits_a1 - 1)   # % change 
+
+# transform independent variables
+
+# set up my_lag
+lag_names <- paste0("a", 1:4)
+
+my_lag <- setNames(paste0("lag(.,", 1:4, ")"),
+                   lag_names)
+
+# create lags
+var_names <- c('rdi', 'hpi')
+data0 <- data0 %>% mutate_at(var_names, funs_(my_lag))
+
+
+# create differences
+data0 <- data0 %>% mutate(rdi_d1 = rdi - rdi_a1)
+data0 <- data0 %>% mutate(hpi_d1 = hpi - hpi_a1)
+
+
+# create lags of the differences
+var_names_diff <- c('rdi_d1', 'hpi_d1')
+data0 <- data0 %>% mutate_at(var_names_diff, funs_(my_lag))
+
+
+# create % changes
+data0 <- data0 %>% mutate(rdi_p1 = rdi/rdi_a1 - 1)
+data0 <- data0 %>% mutate(hpi_p1 = hpi/hpi_a1 - 1)
+
+
+# create lags of the % changes
+var_names_p <- c('rdi_p1', 'hpi_p1')
+data0 <- data0 %>% mutate_at(var_names_p, funs_(my_lag))
+
+
+
+# Create seasonal dummy variables ----
+
+# S1.extract months
+data0$date <- data0$date %>% as.Date('%m/%d/%Y')
+data0$m <- data0$date %>% months
+
+# S2.create dummy vars 
+dumv <- model.matrix(~m, data0)
+dumv <- dumv %>% as.data.frame %>% select(-'(Intercept)')
+
+# S3.put dummy vars into data0
+data0 <- cbind(data0, dumv)
+
+
+
+# Test for seasonality ----
+actual.na <- data0 %>% filter(scenario=='actuals')       # put 'actuals' into a new dataset
+wo(actual.na$Deposits, freq = 4) %>% summary
+
+
+
+# Detecting trends in time-series data ----
+actual.na %>% select(Deposits) %>% ts.plot
+actual.na %>% select(Deposits_d1) %>% ts.plot
+actual.na %>% select(Deposits_p1) %>% ts.plot
+
+
+
+# Augmented Dickey-Fuller test for unit root ----
+actual <- actual.na %>% na.omit                          # delete missing values
+
+# ADF test
+adf.test(actual$Deposits, k = 4)
+adf.test(actual$Deposits_p1, k = 4)
+
+
+
+# Model selection review ----
+
+# create the starting model
+reg0 <- lm(Deposits_p1 ~ mMarch + mJune + mSeptember, data = actual)
+summary(reg0)
+
+# run model selection
+reg1 <- step(reg0, 
+             direction = "forward",
+             scope = ( ~ mMarch + mJune + mSeptember + 
+                         rdi_d1 + rdi_d1_a1 + rdi_d1_a2 + rdi_d1_a3 + rdi_d1_a4 + 
+                         hpi_p1 + hpi_p1_a1 + hpi_p1_a2 + hpi_p1_a3 + hpi_p1_a4))
+summary(reg1)
+
+
+# see VIFs of the variables
+vif(reg1)
+
+# DW test for autocorrelation
+dwtest(reg1)
+
+# Breusch-Godfrey test for autocorrelation
+bgtest(reg1)
+
+# histogram of residuals for detecting violation of normality assumption
+hist(residuals(reg1))
+
+# plot residuals vs fitted values for detecting heteroscedasticity 
+plot(reg1, 1)
+plot(reg1)
+
+
+
+# Present results with robust s.e. ----
+
+# calculate standard s.e. and robust s.e.
+ols_se <- reg1 %>% vcov %>% diag %>% sqrt
+hc_se <- reg1 %>% vcovHC(type = "HC2") %>% diag %>% sqrt
+
+# present results
+stargazer(reg1, reg1, 
+          se = list(ols_se, hc_se),
+          type = 'text')
+
+
+
+# Predict deposits based on the selected model ----
+
+# merge scenarios with actuals.na
+fbase <- bind_rows(actual.na, data0[data0$scenario=="fbase", ])
+fadvr <- bind_rows(actual.na, data0[data0$scenario=="fadvr", ])
+
+# get out-of-sample predictions
+fbase$p_fbase_p1 <- predict(reg1, fbase)
+fadvr$p_fadvr_p1 <- predict(reg1, fadvr)
+
+# delete un-wanted vars 
+fbase <- fbase %>% select(p_fbase_p1, date, scenario, Deposits)
+fadvr <- fadvr %>% select(p_fadvr_p1, date, scenario, Deposits)
+
+# export result
+write_csv(fbase, 'fbase.csv')
+write_csv(fadvr, 'fadvr.csv')
+
+
+
+