knitr::opts_chunk$set(echo = TRUE)
After loading the data we need to clean the data in order to have consistent frequencies and data types. This will allow us to join the datasets.
Korea.int <- subset(Korea.int, select = -c(DATE, Interest.Rate))
Korea.int$Year <- as.numeric(Korea.int$Year)
Korea.int <- na.omit(Korea.int)
colnames(Korea.int)[2] <- "KR.Int.Year.Average"
US.int <- subset(US.int, select = -c(DATE, FEDFUNDS))
US.int$Year <- as.numeric(US.int$Year)
US.int <- na.omit(US.int)
colnames(US.int)[2] <- "US.Int.Year.Average"
Interest.Dif <- inner_join(Korea.int, US.int, by = "Year")
KR.unem$Year <- as.numeric(KR.unem$Year)
colnames(KR.unem)[2] <- "KR.P.Unem"
US.unem$Year <- as.numeric(US.unem$Year)
colnames(US.unem)[2] <- "US.P.Unem"
Unem.Dif <- inner_join(KR.unem,US.unem, by = "Year")
KR.US.inflation$Year <- as.numeric(KR.US.inflation$Year)
US.KR.bus.conf$Year <- as.numeric(US.KR.bus.conf$Year)
colnames(US.KR.bus.conf)[which(names(US.KR.bus.conf) == "USA")] <- "US.bus.Conf"
colnames(US.KR.bus.conf)[which(names(US.KR.bus.conf) == "KR")] <- "KR.bus.Conf"
US.KR.CAPGDP$Year <- as.numeric(US.KR.CAPGDP$Year)
colnames(US.KR.CAPGDP)[2] <- "USA.CAPGDP"
colnames(US.KR.CAPGDP)[3] <- "KR.CAPGDP"
US.KR.comp.ind <- subset(US.KR.comp.ind, select = -c(Leading.Indicator, Indicator))
US.KR.comp.ind$Year <- as.numeric(US.KR.comp.ind$Year)
US.KR.comp.ind <- na.omit(US.KR.comp.ind)
colnames(US.KR.comp.ind)[2] <- "US.Comp.Ind"
colnames(US.KR.comp.ind)[3] <- "KR.Comp.Ind"
US.KR.cons.conf <- subset(US.KR.cons.conf, select = -c(Year.1, US, KR.1))
US.KR.cons.conf $Year <- as.numeric(US.KR.cons.conf$Year)
US.KR.cons.conf <- na.omit(US.KR.cons.conf)
colnames(US.KR.cons.conf)[which(names(US.KR.cons.conf) == "USA")] <- "US.Cons.Conf"
colnames(US.KR.cons.conf)[which(names(US.KR.cons.conf) == "KR")] <- "KR.Cons.Conf"
US.KR.exch <- subset(US.KR.exch, select = -c(Date, Price, Change.Percent, Factor))
US.KR.exch$Year <- as.numeric(US.KR.exch$Year)
US.KR.exch <- na.omit(US.KR.exch)
US.KR.GDP$Year <- as.numeric(US.KR.GDP$Year)
US.KR.GDP$US.GDP.Millions.of.dollars. <- as.numeric(US.KR.GDP$US.GDP.Millions.of.dollars.)
US.KR.GDP$KR.GDP.Millions.of.dollars. <- as.numeric(US.KR.GDP$KR.GDP.Millions.of.dollars.)
After cleaning and formatting data, to make coding easier, I will create one large dataset that contains the inner join of all the datasets by Year. This way I contain a single dataset from which I can pull values from and I do not have to worry about missing data.
Ind.Variables <- US.KR.cons.conf %>%
inner_join(US.KR.GDP, by = 'Year') %>%
inner_join(US.KR.CAPGDP, by = 'Year') %>%
inner_join(US.KR.comp.ind, by = 'Year') %>%
inner_join(US.KR.bus.conf, by = 'Year') %>%
inner_join(KR.US.inflation, by = 'Year') %>%
inner_join(Unem.Dif, by = 'Year') %>%
inner_join(Interest.Dif, by = 'Year')
US.KR.exch <- subset(US.KR.exch, Year >= 1999 & Year <= 2021)
Now, I need to perform the OLS regression. As one final step in the preparation, I will transform the dataset into matrices so I can perform the linear algebra.
# First add the one column into my initial matrix. This way I do not have to manually input it into the matrix.
Ind.Variables$Ones <- 1
Ind.Variables <- Ind.Variables[, c(1,20,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19)]
#Now I create the Matrix X
X <- data.matrix(subset(Ind.Variables, select = -c(Year)))
Y <- data.matrix(subset(US.KR.exch, select = -c(Year, YPCFactor)))
Finally, we can begin the OLS Regression Analysis.
Beta.Hat.1 <- solve(t(X) %*% X, tol = 1e-30) %*% t(X) %*% Y
Y.hat.1 <- X %*% Beta.Hat.1
e.hat.1 <- Y-Y.hat.1
sig.2.hat.1 <- (t(e.hat.1) %*% e.hat.1)/(2)
Var.beta.hat.1 <- as.vector(sig.2.hat.1) * solve(t(X) %*% X, tol = 1e-30)
SE.beta.hat.1 <- sqrt(diag(Var.beta.hat.1))
Z.values.1 <- Beta.Hat.1 / SE.beta.hat.1
P.values.1 <- pnorm(Z.values.1)
R.squared.1 <- 1 - sum((Y-Y.hat.1)^2)/sum((Y-mean(Y))^2)
R.squared.1
## [1] 0.8635023
P.values.1
## YPC
## Ones 0.3773681
## US.Cons.Conf 0.5465204
## KR.Cons.Conf 0.6810842
## US.GDP.Millions.of.dollars. 0.2641259
## KR.GDP.Millions.of.dollars. 0.7543938
## US_Growth 0.4746029
## KR_Growth 0.2641480
## USA.CAPGDP 0.2918098
## KR.CAPGDP 0.5581806
## US.Comp.Ind 0.6544395
## KR.Comp.Ind 0.2211706
## US.bus.Conf 0.1196342
## KOREA 0.7483717
## US.Inflation.Rate 0.5818821
## KR.Inflation.Rate 0.2886928
## KR.P.Unem 0.7892924
## US.P.Unem 0.4746560
## KR.Int.Year.Average 0.7835719
## US.Int.Year.Average 0.2298402
Using an alpha value of 5%, there are no significantly significant predictors in this model. Now lets compare it to our benchmark.
Y.hat <- mean(US.KR.exch$YPC)
e.hat <- Y-Y.hat
Se.1 <- sd(e.hat)/sqrt(length(e.hat))
Z.values.benchmark <- Y.hat.1 / Se.1
P.values.benchmark <- pnorm(Z.values.benchmark)
Average.variance <- sqrt(sum(e.hat^2) / nrow(e.hat))
P.values.benchmark
## YPC
## 1 1.430219e-01
## 2 9.927053e-01
## 3 9.997999e-01
## 4 3.797832e-04
## 5 4.523471e-03
## 6 6.273108e-04
## 7 3.444996e-02
## 8 7.403964e-03
## 9 9.977576e-01
## 10 1.000000e+00
## 11 2.847737e-10
## 12 2.362554e-04
## 13 9.828649e-01
## 14 9.588718e-01
## 15 9.180788e-03
## 16 6.908679e-01
## 17 9.967971e-01
## 18 6.163560e-01
## 19 1.248921e-04
## 20 5.467794e-01
## 21 9.999628e-01
## 22 2.835067e-03
## 23 9.995696e-01
Both Models are not significant and both are not helpful in answering our question. What about some other models?
For a few alternative models, I repeated the steps I did with the first model, except I took out the variables with the least amount of significance (highest P-values) and re-ran the analysis.
X.2 <- data.matrix(subset(Ind.Variables, select = -c(Year, KOREA, KR.GDP.Millions.of.dollars., KR.Int.Year.Average)))
Beta.Hat.2 <- solve(t(X.2) %*% X.2, tol = 1e-30) %*% t(X.2) %*% Y
Y.hat.2 <- X.2 %*% Beta.Hat.2
e.hat.2 <- Y-Y.hat.2
sig.2.hat.2 <- (t(e.hat.2) %*% e.hat.2)/(nrow(X.2)-ncol(X.2))
Var.beta.hat.2 <- as.vector(sig.2.hat.2) * solve(t(X.2) %*% X.2, tol = 1e-30)
SE.beta.hat.2 <- sqrt(diag(Var.beta.hat.2))
Z.values.2 <- Beta.Hat.2 / SE.beta.hat.2
P.values.2 <- pnorm(Z.values.2)
R.squared.2 <- 1 - sum((Y-Y.hat.2)^2)/sum((Y-mean(Y))^2)
R.squared.2
## [1] 0.8194854
P.values.2
## YPC
## Ones 0.94866492
## US.Cons.Conf 0.18177674
## KR.Cons.Conf 0.32170145
## US.GDP.Millions.of.dollars. 0.50072692
## US_Growth 0.65095466
## KR_Growth 0.74140524
## USA.CAPGDP 0.71204918
## KR.CAPGDP 0.40000826
## US.Comp.Ind 0.60651446
## KR.Comp.Ind 0.38292970
## US.bus.Conf 0.06222657
## US.Inflation.Rate 0.51155264
## KR.Inflation.Rate 0.50226018
## KR.P.Unem 0.73564704
## US.P.Unem 0.19634749
## US.Int.Year.Average 0.33558157
X.3 <- data.matrix(subset(Ind.Variables, select = -c(Year, KOREA, KR.GDP.Millions.of.dollars., KR.Int.Year.Average, KR_Growth, US.GDP.Millions.of.dollars., USA.CAPGDP)))
Beta.Hat.3 <- solve(t(X.3) %*% X.3, tol = 1e-30) %*% t(X.3) %*% Y
Y.hat.3 <- X.3 %*% Beta.Hat.3
e.hat.3 <- Y-Y.hat.3
sig.2.hat.3 <- (t(e.hat.3) %*% e.hat.3)/(nrow(X.3)-ncol(X.3))
Var.beta.hat.3 <- as.vector(sig.2.hat.3) * solve(t(X.3) %*% X.3, tol = 1e-30)
SE.beta.hat.3 <- sqrt(diag(Var.beta.hat.3))
Z.values.3 <- Beta.Hat.3 / SE.beta.hat.3
P.values.3 <- pnorm(Z.values.3)
R.squared.3 <- 1 - sum((Y-Y.hat.3)^2)/sum((Y-mean(Y))^2)
R.squared.3
## [1] 0.8068098
P.values.3
## YPC
## Ones 0.9933949062
## US.Cons.Conf 0.1929830866
## KR.Cons.Conf 0.5310456671
## US_Growth 0.4415626736
## KR.CAPGDP 0.5581081556
## US.Comp.Ind 0.6876542543
## KR.Comp.Ind 0.2024768639
## US.bus.Conf 0.0005495937
## US.Inflation.Rate 0.8072833776
## KR.Inflation.Rate 0.4688601769
## KR.P.Unem 0.9433211432
## US.P.Unem 0.2040291985
## US.Int.Year.Average 0.1997281172
X.4 <- matrix(cbind(Ind.Variables$Ones, Ind.Variables$US.bus.Conf), ncol = 2)
Beta.hat.4 <- solve(t(X.4) %*% X.4, tol = 1e-40) %*% t(X.4) %*% Y
Y.hat.4 <- X.4 %*% Beta.hat.4
e.hat.4 <- Y-Y.hat.4
sig.2.hat.4 <- (t(e.hat.4)%*% e.hat.4)/(nrow(X.4) - ncol(X.4))
Var.beta.hat.4 <- as.vector(sig.2.hat.4) * solve(t(X.4) %*% X.4, tol = 1e-40)
SE.beta.hat.4 <- sqrt(diag(Var.beta.hat.4))
Z.values.4 <- Beta.hat.4 / SE.beta.hat.4
P.values.4 <- pnorm(Z.values.4)
R.squared.4 <- 1 - sum((Y-Y.hat.4)^2)/sum((Y-mean(Y))^2)
R.squared.4
## [1] 0.550242
P.values.4
## YPC
## [1,] 9.999998e-01
## [2,] 2.002646e-07
Now lets do some graphical analysis of our regressions. Because most of our models are multivariate in nature, we will only be examining the residual plots of each model, including the benchmark.
e.hat <- as.data.frame(e.hat)
e.hat$Number <- seq(1,23)
e.hat.1 <- as.data.frame(e.hat.1)
e.hat.1$Number <- seq(1,23)
e.hat.2 <- as.data.frame(e.hat.2)
e.hat.2$Number <- seq(1,23)
e.hat.3 <- as.data.frame(e.hat.3)
e.hat.3$Number <- seq(1,23)
e.hat.4 <- as.data.frame(e.hat.4)
e.hat.4$Number <- seq(1,23)
ggplot(e.hat) + geom_point(mapping = aes(Number, YPC), size = 3, color = "firebrick2") + ylab("Residual") + ggtitle("Residual of Benchmark") + geom_hline(aes(yintercept = 0), size = 1)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggplot(e.hat.1) + geom_point(mapping = aes(Number, YPC), size = 3, color = "firebrick2") + ylab("Residual") + ggtitle("Residual of Model 1") + geom_hline(aes(yintercept = 0), size = 1)
ggplot(e.hat.2) + geom_point(mapping = aes(Number, YPC), size = 3, color = "firebrick2") + ylab("Residual") + ggtitle("Residual of Model 2") + geom_hline(aes(yintercept = 0), size = 1)
ggplot(e.hat.3) + geom_point(mapping = aes(Number, YPC), size = 3, color = "firebrick2") + ylab("Residual") + ggtitle("Residual of Model 3") + geom_hline(aes(yintercept = 0), size = 1)
ggplot(e.hat.4) + geom_point(mapping = aes(Number, YPC), size = 3, color = "firebrick2") + ylab("Residual") + ggtitle("Residual of Model 4") + geom_hline(aes(yintercept = 0), size = 1)
Let’s remove the outliers in the plot.
e.hat.1$Model <- "Model 1"
e.hat.2$Model <- "Model 2"
e.hat.3$Model <- "Model 3"
e.hat.4$Model <- "Model 4"
e.hat$Model <- "Benchmark"
Residuals <- rbind(e.hat, e.hat.1, e.hat.2, e.hat.3, e.hat.4)
#Lets plot Box and Whisker Plots to identify the outliers
ggplot(Residuals) + geom_boxplot((mapping = aes(x = Model, y = YPC)), fill = "firebrick2") + theme(axis.ticks.x = element_blank()) + ylab ("Residuals") + xlab("")
#Without the outliers
ggplot(Residuals) + geom_boxplot((mapping = aes(x = Model, y = YPC)), fill = "firebrick2", outlier.shape = NA) + theme(axis.ticks.x = element_blank()) + ylab ("Residuals") + xlab("") + coord_cartesian(ylim = c(-20, 20))
Now, lets recreate the residual plots without the outliers.
ggplot(Residuals) + geom_point(mapping = aes(Number, YPC), size = 1, color = "firebrick2") + facet_wrap(~Model) + ylab("Residual") + xlab("Number") + geom_hline(aes(yintercept = 0), size = 1) + coord_cartesian(ylim = c(-20, 20))
Now lets create QQplots for all the models, excluding outliers.
Residuals <- Residuals[!(Residuals$Number == "10" & Residuals$Model == "Benchmark"),]
Residuals <- Residuals[!(Residuals$Number == "5" & Residuals$Model == "Model 2"),]
Residuals <- Residuals[!(Residuals$Number == "5" & Residuals$Model == "Model 3"),]
Residuals <- Residuals[!(Residuals$Number == "10" & Residuals$Model == "Model 4"),]
ggplot(Residuals, mapping = aes(sample = YPC)) + stat_qq(shape = 1, color = "firebrick2") + stat_qq_line(color = "royalblue2") + facet_wrap(~Model) + ylab("Sample Quantities") + xlab("Theoretical Quantities") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + theme_gray() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())