Problem Set 1.

  1. Set working directory and load approval.csv. Explore the data set and select only month, year and approve columns. Describe the range of the time series. Do you think the data set contains enough observations?
# 1. Загрузка данных и выбор колонок
setwd("/Users/andreii/Documents/R_files/ds_hw1")
data <- read.csv("approval.csv")
data <- data %>% select(month, year, approve)
summary(data)
##      month             year         approve     
##  Min.   : 1.000   Min.   :2001   Min.   :35.67  
##  1st Qu.: 3.000   1st Qu.:2002   1st Qu.:48.50  
##  Median : 6.000   Median :2003   Median :54.67  
##  Mean   : 6.308   Mean   :2003   Mean   :57.17  
##  3rd Qu.: 9.000   3rd Qu.:2005   3rd Qu.:65.20  
##  Max.   :12.000   Max.   :2006   Max.   :88.00
  1. (0.5 points) Create a date column with paste0 function or any way you come up with. Apply the as.Date function and assign %Y-%m-%d date format.
data$date <- as.Date(paste0(data$year, "-", data$month, "-01"), "%Y-%m-%d")
  1. (0.5 points) Visualise the presidential approval time series with ggplot2. Add a mean line. Visually inspect the series and guess whether it is stationary or not. Justify your answer with the appropriate statistics.
data %>% 
  ggplot(aes(x = date, y = approve)) + 
  geom_line() + 
  geom_hline(yintercept = mean(data$approve, na.rm = TRUE), col = "red")

  1. (1 point) Conduct unit-root tests. Briefly describe the idea behind them. Report all the results in one table, e.g. which unit root tests you use, what the respective null hypotheses are and what the tests’ verdict is. Compare them to what you initially guessed.
adf_result <- adf.test(na.omit(data$approve))
pp_result <- pp.test(na.omit(data$approve))
kpss_result <- kpss.test(na.omit(data$approve))
## Warning in kpss.test(na.omit(data$approve)): p-value smaller than printed
## p-value
results_table <- data.frame(
  Test = c("ADF", "PP", "KPSS"),
  H0 = c("Non-stationary", "Non-stationary", "Stationary"),
  p_value = c(adf_result$p.value, pp_result$p.value, kpss_result$p.value),
  Decision = c(ifelse(adf_result$p.value < 0.05, "Stationary", "Non-stationary"),
               ifelse(pp_result$p.value < 0.05, "Stationary", "Non-stationary"),
               ifelse(kpss_result$p.value < 0.05, "Non-stationary", "Stationary"))
)

results_table
  1. (1 point) Use ACF and PACF to visually inspect the presence of MA and AR components. What would you conclude from the ACF graph? From PACF? Write it down in a couple of sentences. Guess \(p, i, q\) parameters of the would-be ARIMA model. Fit the model with the parameters chosen. Visualise observed time series and predicted values from your ARIMA(\(p,i,q\)) model. Describe how well your model fits the observed data.
par(mfrow = c(1, 2))
acf(na.omit(data$approve))
pacf(na.omit(data$approve))

arima_model <- arima(data$approve, c(1,0,1))
data$fitted <- fitted(arima_model)

ggplot(data, aes(x = date)) +
  geom_line(aes(y = approve), col = "black") +
  geom_line(aes(y = fitted), col = "red")

  1. (0.5 points) Decide on how many times you should difference your data. Plot the differenced data. Describe irregularities in behaviour of your time series.
data$approve_diff <- c(NA, diff(data$approve))
plot.ts(na.omit(data$approve_diff))

  1. (0.5 points) Inspect the differenced data on presence of AR/MA components. What would the appropriate model look like? Decide on \(p,i,q\). Do you find something suspicious? What’s the process your differenced data follows?
par(mfrow = c(1, 2))
acf(na.omit(data$approve_diff))
pacf(na.omit(data$approve_diff))

  1. Read approval gallup.dta. Examine the data set. Create an appropriate date column. Rename the column with an approval share to approve.
gallup <- read_dta("/Users/andreii/Documents/R_files/ds_hw1/approval_gallup.dta")
skim(gallup)
Data summary
Name gallup
Number of rows 541
Number of columns 13
_______________________
Column type frequency:
character 6
numeric 7
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
pres 0 1 4 10 0 9 0
poll 0 1 6 6 0 1 0
polltype 0 1 0 8 533 2 0
month2 0 1 1 2 0 12 0
year2 0 1 4 4 0 48 0
time 0 1 6 7 0 541 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
month 0 1.00 6.42 3.45 1 3 6 9 12 ▇▆▅▅▇
year 0 1.00 1976.79 13.95 1953 1965 1977 1989 2000 ▇▇▇▇▇
papp 0 1.00 56.19 11.94 24 48 57 65 89 ▁▅▇▅▁
disapprove 0 1.00 31.81 12.77 2 22 32 41 65 ▂▆▇▅▁
none 0 1.00 12.03 4.78 2 9 12 14 43 ▆▇▁▁▁
pollsize 123 0.77 1400.14 318.08 615 1039 1518 1558 3077 ▃▇▁▁▁
date 0 1.00 206.94 167.51 -83 63 212 355 491 ▇▇▇▇▇ 
gallup$date <- as.Date(paste0(gallup$year, "-", gallup$month, "-01"), "%Y-%m-%d")
names(gallup)[names(gallup) == "papp"] <- "approve"

  1. (1.5 points) Do all the steps from the previous part on the approval_gallup data set. Inspect time series properties both visually and technically (employing formal unit root tests). Describe all the steps by commenting the code.

    1. Visualise plot (judge on the stationarity with the plot)
    2. Run unit-root tests. Compare them to each other and your visual guess from the step 1.
    3. Run ACF/PACF to pick up the plausible components for AR, I, MA. Run ARIMA with the paramaters chosen.
    4. Run auto-arima to find the best model. Compare fits of the model of your choice and the model proposed by the auto-arima, both with quality metrics and visualisation of fitted values’ fit to the observed ones.
# 9.1 Визуализация
gallup %>% 
  ggplot(aes(x = date, y = approve)) + 
  geom_line()

# 9.2 Unit-root тесты
adf.test(na.omit(gallup$approve))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(gallup$approve)
## Dickey-Fuller = -3.7492, Lag order = 8, p-value = 0.02154
## alternative hypothesis: stationary
pp.test(na.omit(gallup$approve))
## Warning in pp.test(na.omit(gallup$approve)): p-value smaller than printed
## p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  na.omit(gallup$approve)
## Dickey-Fuller Z(alpha) = -52.203, Truncation lag parameter = 6, p-value
## = 0.01
## alternative hypothesis: stationary
kpss.test(na.omit(gallup$approve))
## Warning in kpss.test(na.omit(gallup$approve)): p-value smaller than printed
## p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(gallup$approve)
## KPSS Level = 1.2296, Truncation lag parameter = 6, p-value = 0.01
# 9.3 ACF/PACF и ARIMA
par(mfrow = c(1, 2))
acf(na.omit(gallup$approve))
pacf(na.omit(gallup$approve))

manual_arima <- arima(gallup$approve, c(1,1,1))

# 9.4 Auto-ARIMA
auto_model <- auto.arima(gallup$approve)
AIC(manual_arima)
## [1] 3326.202
AIC(auto_model)
## [1] 3324.202

  1. (1 point) Now join two data sets in one called ts_approval. Visualise it. Do you still believe it is stationary? In fact, there are plenty of missings in the resulting data set. It might cause some problems with data wrangling and model fitting. For example, you have a missing for 2001-01-01. Try to insert 65% presidential approval to the time point. Let’s try to find all missings in the time series:

    1. Order your dataset by date.
    2. Use missingMonths lapply algorithm to find the missed months. Try to figure out what’s going on under the hood of the function.
    3. Create a full_vector vector by running the following sequence: seq(start, by='1 month', length=<number of months in a period of interest>). This should result in sequential list of months from February 1953 to June 2006.
    4. Create a data.frame df_nans. Adapt the code: data.frame(date=full_vector, approve=with(<your merged data set>, approve[match(full_vector, date)])). Describe the patterns in missing data.
    5. Make the approve column a timeseries (ts) object with an appropriate time span.
    6. Apply a function ggplot_na_distribution(<your ts>) from imputeTS package to your time series to inspect the missigness pattern. Could you believe the missings are at random?
    7. Impute the missed observations with na_mean algorithm from the same imputeTS package. Are you satisfied with a quality of imputation? Find the better algorithm from the cheat sheet provided in the Google Drive class folder. Proceed with the best imputation technique.
# 10. Объединение датасетов
ts_approval <- full_join(data, gallup, by = "date", suffix = c("_1", "_2"))
ts_approval$approve <- coalesce(ts_approval$approve_1, ts_approval$approve_2)
ts_approval <- ts_approval %>% select(date, approve) %>% arrange(date)

# Вставка значения для 2001-01-01
ts_approval[ts_approval$date == "2001-01-01", "approve"] <- 65

# Поиск пропущенных месяцев
start_date <- min(ts_approval$date, na.rm = TRUE)
end_date <- max(ts_approval$date, na.rm = TRUE)
full_vector <- seq(start_date, end_date, by = "1 month")

df_nans <- data.frame(
  date = full_vector,
  approve = with(ts_approval, approve[match(full_vector, date)])
)

# Создание ts объекта и импутация
ts_obj <- ts(df_nans$approve, start = c(1953, 2), frequency = 12)
ggplot_na_distribution(ts_obj)

# Импутация
ts_imputed <- na_interpolation(ts_obj, option = "spline")
df_nans$approve_imputed <- as.numeric(ts_imputed)

URL: https://cran.rstudio.com/web/packages/imputeTS/vignettes/Cheat_Sheet_imputeTS.pdf

  1. (1 point) Run the ARIMA model on the imputed data set (set the seasonal argument to FALSE). Choose the \(p,i,q\), based on the plots and tests you used above.
arima_imputed <- arima(ts_imputed, c(1,1,1))
  1. (0.5 points) Compare the model of your choice and the model provided by auto-arima by residual diagnostics. Describe the behaviour of residuals. What’s unusual in the distribution of the residuals?
auto_imputed <- auto.arima(ts_imputed, seasonal = FALSE)

# Диагностика остатков
checkresiduals(arima_imputed)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)
## Q* = 28.205, df = 22, p-value = 0.169
## 
## Model df: 2.   Total lags used: 24
checkresiduals(auto_imputed)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)
## Q* = 27.971, df = 23, p-value = 0.2169
## 
## Model df: 1.   Total lags used: 24
  1. (1 point) Load econ_join.csv or time.csv dataset. Check missings and impute if necessary. Come up with any idea on the temporal relationship between two variables. First, build the distributed-lag model (DLM) and interpret the coefficients: 1) impact effect, 2) cumulative dynamic multipliers and 3) long-run cumulative dynamic multiplier. Then, run autoregressive distributed-lag (ADL) model and interpret the coefficients.
econ <- read.csv("/Users/andreii/Documents/R_files/ds_hw1/econ_join.csv", sep = ";")

skim(econ)
Data summary
Name econ
Number of rows 148
Number of columns 12
_______________________
Column type frequency:
character 1
numeric 11
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
Date 0 1 10 10 0 148 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
X 0 1.00 74.50 42.87 1.0 37.75 74.50 111.25 148.0 ▇▇▇▇▇
econ_join 0 1.00 20.03 5.01 8.0 17.00 20.00 23.25 33.0 ▂▆▇▅▁
econ_join_lag 7 0.95 20.08 5.06 8.0 17.00 20.00 24.00 33.0 ▂▆▇▅▁
net_way 10 0.93 -1.59 25.96 -71.0 -12.00 1.50 15.00 42.0 ▁▁▃▇▃
net_way_lag 14 0.91 -1.70 26.30 -74.0 -14.00 2.00 16.00 44.0 ▁▁▃▇▃
unemployment 10 0.93 7.55 2.29 4.0 5.60 7.15 8.60 14.1 ▇▇▅▂▁
unemployment_lag 9 0.94 7.63 2.31 4.1 5.60 7.30 9.00 14.6 ▇▇▅▃▁
cpi 10 0.93 1.25 3.35 -0.3 0.40 0.70 1.30 38.4 ▇▁▁▁▁
cpi_lag 9 0.94 1.01 1.22 -0.5 0.40 0.80 1.20 11.6 ▇▁▁▁▁
misery 10 0.93 8.80 4.60 4.1 6.20 8.00 10.28 50.3 ▇▁▁▁▁
misery_lag 9 0.94 8.65 3.09 4.4 6.30 7.80 10.40 24.8 ▇▅▁▁▁ 
# DLM модель
dlm <- lm(econ_join ~ lag(unemployment, 1) + lag(unemployment, 2) + lag(unemployment, 3), data = econ)
summary(dlm)
## 
## Call:
## lm(formula = econ_join ~ lag(unemployment, 1) + lag(unemployment, 
##     2) + lag(unemployment, 3), data = econ)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.1123  -3.1799  -0.4367   2.5088  12.9853 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           12.4588     1.3531   9.208 6.61e-16 ***
## lag(unemployment, 1)   1.2719     0.6084   2.091   0.0385 *  
## lag(unemployment, 2)   0.6919     0.8132   0.851   0.3964    
## lag(unemployment, 3)  -0.9572     0.6098  -1.570   0.1189    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.476 on 132 degrees of freedom
##   (12 observations deleted due to missingness)
## Multiple R-squared:  0.2406, Adjusted R-squared:  0.2234 
## F-statistic: 13.94 on 3 and 132 DF,  p-value: 5.954e-08
# ADL модель  
adl <- lm(econ_join ~ lag(unemployment, 1) + lag(unemployment, 1) + lag(unemployment, 2), data = econ)
summary(adl)
## 
## Call:
## lm(formula = econ_join ~ lag(unemployment, 1) + lag(unemployment, 
##     1) + lag(unemployment, 2), data = econ)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.7697 -3.2884 -0.5057  2.8459 12.9976 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           12.1668     1.3402   9.078 1.24e-15 ***
## lag(unemployment, 1)   1.1716     0.6076   1.928   0.0559 .  
## lag(unemployment, 2)  -0.1342     0.6101  -0.220   0.8262    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.492 on 134 degrees of freedom
##   (11 observations deleted due to missingness)
## Multiple R-squared:  0.2237, Adjusted R-squared:  0.2121 
## F-statistic: 19.31 on 2 and 134 DF,  p-value: 4.274e-08
  1. (1 point) Use time.csv data set. With dynlm package and some others from seminar on orange juice, run the dynamic linear model, explaining a share of those who think the country goes wrong way by unemployment. Choose an appropriate number of lags and plot dynamic and cumulative effects.
library(dynlm)

# Подготовка данных
time_data <- read.csv("/Users/andreii/Documents/R_files/ds_hw1/time.csv", sep=";")
ts_wrongway <- ts(time_data$wrongway, start = c(2000, 1), frequency = 12)
ts_unemp <- ts(time_data$unemployment, start = c(2000, 1), frequency = 12)

# Динамическая модель
dyn_model <- dynlm(ts_wrongway ~ L(ts_unemp, 0:4))

# Динамические эффекты
dynamic_effects <- coef(dyn_model)[-1]
cumulative_effects <- cumsum(dynamic_effects)

# Визуализация
par(mfrow = c(1, 2))
plot(0:4, dynamic_effects, type = "b", main = "Dynamic Effects")
plot(0:4, cumulative_effects, type = "b", main = "Cumulative Effects")